Conservation of Retrotransposon Fragments in Animal Genomes

1. Introduction

Population genetics of livestock species addresses several practical problems, including the identification of markers associated with adaptation, productivity, reproduction, and the conservation of biodiversity. Methods used to address these problems range from phenotypic characteristics, including A.S. Serebrovsky’s “signalia” [1], to modern applications of artificial intelligence [2]. Major challenges include identifying markers associated with phenotypic variability and understanding the complex interplay among adaptation, productivity, and reproduction, including interactions among the metabolic pathways underlying these traits.

The search for such genomic markers has remained an important problem for decades. Advances in genomics and epigenomics have made regulatory sequences a promising source of markers. Transposable elements (TEs) and their recombination products contribute substantially to these regulatory sequences [3]. The relationship between TEs and their hosts is complex: TEs contain nucleotide motifs recognized by host regulatory factors that support TE activity, while host genomes contain mechanisms that limit the potentially detrimental genomic effects of TEs. These interactions involve genetic exchange and rapid sequence evolution, with evolutionarily conserved sequences reflecting successful long-term interactions between TEs and their hosts [4].

Mammalian genomes, almost half of which consist of TEs, have been described as TE-rich ecosystems in which transposable elements contribute to gene regulatory networks and genome evolution [5]. Genome size, the abundance of noncoding sequences, and TE content are closely associated with the evolution of multicellular complexity [6], while the emergence of phenotypic innovations has also been associated with TE expansion [7,8]. The distribution of different TE families in livestock genomes is species-specific [9]. The evolutionary age of TEs contributing to regulatory elements such as promoters and enhancers is also related to the functions of their target protein-coding genes: regulatory TEs associated with housekeeping genes tend to be older than those associated with tissue-specific genes [10]. TE-associated noncoding sequences evolve rapidly. Estimates of noncoding DNA gain and loss over approximately 90 million years of human–mouse divergence suggest that only about 40% of the present-day human genome shares common ancestry with the mouse genome; thus, at least 60% of noncoding sequences have changed in each lineage [11]. Moreover, even laboratory mouse strains of common origin differ in their TE composition [12].

Despite their extensive diversity and widespread horizontal transfer [13,14,15], individual TEs can contribute directly to complex phenotypic traits. For example, specific TEs have been shown to play a critical role in regulating neural crest cell migration and, consequently, craniofacial development [16]. Of particular interest are TEs that retain conserved domains across different endogenous retroviruses [17] and TEs involved in the same biological processes, such as the formation of topologically associating domains (TADs) and similar regulatory motifs in paralogous and orthologous genomic regions [18,19]. TEs involved in the evolution of the placenta provide one example. Four major types of mammalian placenta are distinguished according to the degree of syncytialization and the distribution of chorionic villi: epitheliochorial (diffuse) in pigs and horses, synepitheliochorial (cotyledonary) in cattle and sheep, endotheliochorial (zonary) in dogs and cats, and hemochorial (discoid) in humans, rabbits, and mice. These placental types are thought to have evolved through independent integrations of exogenous retroviruses (endogenous retroviruses after integration, ERVs) into host genomes [20,21]. Mechanisms behind placental innovations include the reuse of proteins already expressed in a tissue (co-option), induction of genes normally expressed elsewhere in the organism (recruitment), and introduction of new genes through gene duplication or retroviral insertion (horizontal transfer) [20]. It has been proposed that mammals evolved lineage-specific placental morphologies using ERVs independently acquired in each lineage. Peg10 and Peg11/Rtl1, which have fundamental roles in placental function, were already present in ancestral mammals. In cattle, 7,624 ERV-derived genes have been identified, including 1,542 env-derived genes. Of these, 284 are expressed in the placenta, including BERV-K1, BERV-K2, BERV-K3, Syncytin-Rum1, BERV-P, bERVE-A, and bERVE-B. Epigenetic regulation, maternal–fetal conflict, and gene replacement have been proposed to contribute to this expression and, collectively, to the development of ruminant-specific placental morphology [20].

In our previous studies, we used LTR fragments of BERV-K1 and BERV beta3 (BERV beta3 according to [22]), ERV2-1_BT and ERV2-2_BT according to [23] as IRAP primers to genotype genomic regions flanked by inverted copies of these sequences [24,25]. These analyses revealed substantial interspecific and interbreed differences in IRAP profiles of cattle, sheep, and horses [26,27], as well as in domestic rabbits [28] and even honey bees (unpublished data). Because BERV-K1 and BERV beta3 LTR fragments are evolutionarily conserved, we tried to identify DNA motifs recognized by transcription factors that could explain the functional significance of this conservation. The presence of regulatory motifs in the 3′ LTRs of BERV-K1 and BERV-K2 that may facilitate nuclear export of retroviral mRNAs, as observed in other retroviruses [29], was discussed in early studies of these placenta-expressed retroviruses [30]. A number of host transcription factors that bind different regions of TEs and support TE propagation and transcription have been identified [3], while TEs, in turn, may increase the adaptive potential of the host [31]. Comparing the presence and polymorphism of conserved BERV-K1 and BERV beta3 fragments between cattle and rabbits, which have different placental types, may help identify regulatory mechanisms associated with these differences.

2. Materials and Methods

The analysis included the Bos taurus genome assembly ARS-UCD2.0 (2024; GCF_002263795.3, GenBank GCA_002263795.4), the Oryctolagus cuniculus genome assembly mOryCun1.1 (2024; GCF_964237555.1, GenBank GCA_964237555.1), and the Apis cerana genome assembly AcerK_1 (2026; GCA_029169275.1, GCF_029169275.2). GenBank, RepBase, and miRBase resources were used.

Nucleotide sequences from different regions of endogenous retroviruses were compared using the NCBI Global Align tool (Needleman–Wunsch algorithm), while searches for homologous sequences were performed using BLASTn. The presence of recombination products from different viruses within endogenous betaretroviral sequences of cattle and domestic rabbit was assessed using the GIRI Repbase database and CENSOR [32].Transcription factor-binding sites were identified using MEME Suite with the JASPAR CORE database [33].

3. Results

The BERVK1 and BERV-beta3 sequences selected for genotyping different animal species are fragments of the long terminal repeats (LTRs) of these betaretroviruses.

In cattle, BERVK1 is located on chromosome 7 within an intron of FAT2, an atypical cadherin belonging to the protocadherin superfamily and involved in cell adhesion and epithelial–mesenchymal transition [34]. BERVK1 is 10,185 bp long, with the 5′ LTR at positions 1–1285, the 3′ LTR at positions 8898–10185, and the fusogenic env gene at positions 6790–8931 (GenBank: AB751366.1). The fragment selected as a primer (5′-TATCAGGCCTCTCCGCATG-3′ or its reverse complement, 5′-CATGCGGAGAGGCCTGATA-3′) is located at positions 42–60 of the 5′ LTR and 8939–8957 of the 3′ LTR. Sequence identity between the 5′ and 3′ LTRs reaches 98%. The 5′ and 3′ LTRs are 1,285 and 1,288 bp long, respectively. Relative to the 3′ LTR, the 5′ LTR contains three nucleotide deletions, corresponding to positions 537, 697, and 1034 of the 3′ LTR, whereas the 3′ LTR lacks the terminal five nucleotides, AGTGG, present in the 5′ LTR. The two LTRs also differ by eight transversions, four of which occur within positions 542–555 of the 3′ LTR, and 11 transitions. Sequences homologous to the 19-nt primer fragment were detected on all Bos taurus chromosomes, including X and Y, on both DNA strands, with 100% query coverage and 89–100% sequence identity.

Protein-coding transcripts. Within several protein-coding genes (on either DNA strand) we found sequences fully matching the BERVK1 LTR fragment previously used as a genotyping primer. For example, on chromosome 1, the sequence is located on the minus strand within an intron of NLGN1 (Chr1: 93,783,791–93,783,809). NLGN1 encodes neuroligin-1, a postsynaptic cell-adhesion protein located at excitatory glutamatergic synapses and involved in synapse formation, stabilization, and plasticity [35]. On chromosome 2, the fragment is located on the plus strand within an intron of HYCC2 (Chr2: 89,643,832–89,644,002), which encodes hyccin PI4KA lipid kinase complex subunit 2 [36]. HYCC2 is a component of the complex required for localization of phosphatidylinositol 4-kinase (PI4K) to the plasma membrane. The fragment was also found in transcript isoforms of a myotubularin-related gene (ChrX: 95,781,091–96,066,580); myotubularin-related proteins participate in lipid metabolism, actin cytoskeleton remodeling, and vesicular trafficking [37]. Additional matches were found in transcripts of the transcriptional regulator ZFP37 (Chr8: 102,210,192–102,305,553), which is highly expressed in testis and is involved in chromatin organization [38], and in ABO (Chr11: 104,176,830–104,214,825), which encodes a glycosyltransferase involved in blood-group antigen synthesis and cell-adhesion-related processes [39]. Several copies of this fragment were also detected on chromosome Y and in a number of long noncoding RNAs (lncRNAs). For example, a 179,598-bp chromosome Y clone (GenBank: AC216981.4) contains a BERVK1-derived region at positions 14,015–15,026. This region shows 94% identity to positions 1–988 of the BERVK1 LTR and contains the sequence 5′-CATGCGGAGAGGCCTGATA-3′. The env gene product has not only fusogenic activity but is also associated with extracellular vesicles: the BERV-K1 syncytin protein has been detected in bovine placental extracellular vesicles [40].

BERV-beta3 is located on chromosome 5 [22]. Its complete sequence is 7,666 bp long. Its LTRs are considerably shorter than those of BERVK1: the 5′ LTR is 406 bp and the 3′ LTR is 407 bp, and their sequence identity (93%) is lower than that between the BERVK1 LTRs (98%). The two BERV-beta3 LTRs differ by two nucleotide deletions in the 3′ LTR, five transversions, and 11 transitions. Unlike BERVK1, BERV-beta3 contains, in addition to env (positions 5427–7238), the viral protein-coding genes gag (478–2001), pro (1950–2765), and pol (2765–5404). The fragment selected as a primer (5′-GGACCGTCTCCTTCAAGGC-3′ or its reverse complement, 5′-GCCTTGAAGGAGACGGTCC-3′) occupies positions 25–43 of the 5′ LTR and 7285–7303 of the 3′ LTR. Homologous sequences were detected on all B. taurus chromosomes except chromosome Y, on both DNA strands, with 100% query coverage and 81–100% sequence identity. Complete matches were found in several lncRNAs and in transcripts of protein-coding genes. These included the TRA/TRD locus (T-cell receptor genes), including TRDV1-1; an uncharacterized protein-coding locus, LOC112441588, located within an olfactory receptor gene cluster on chromosome 15; HSD17B4 (Chr7: 34,291,146–34,425,550), involved in long-chain fatty-acid and steroid metabolism; and the calcium-dependent secretory phospholipase A2 gene PLA2G2C (Chr2: 132,459,318–132,490,763), with one transition at position 3834 of the transcript.

A number of transcripts matched 15 of the 19 nucleotide positions. These included COX4I1 (Chr18: 11,758,912–11,767,079), encoding a cytochrome c oxidase subunit; CFAP57 (Chr3: 102,759,678–102,839,317), a gene involved in cilia and flagella function; inverted repeats in the mRNA of CATSPERG (Chr18: 48,146,347–48,177,540), which regulates the CatSper calcium channel and sperm hyperactivation during fertilization; RNF31 (Chr10: 20,992,427–21,003,728), encoding an E3 ubiquitin ligase that is the major catalytic component of the linear ubiquitin chain assembly complex (LUBAC) and contributes to NF-κB activation [41]; the RNA polymerase II-associated gene LEO1 (Chr10: 58,440,313–58,458,755); STRA6 (Chr21: 34,478,073–34,507,705), encoding a retinol transporter and signaling receptor; and RYR3 (Chr10: 28,727,820–29,286,475), encoding ryanodine receptor 3, a calcium-release channel.

In the domestic rabbit (Oryctolagus cuniculus), an endogenous retroviral env gene with fusogenic activity and placenta-specific expression has been described within a provirus with the typical LTR–gag–pol–env–LTR organization [42]. In the rabbit genome, syncytin-Ory1 is located on chromosome 13. Homologous sequences were detected on many rabbit autosomes and at several sites on chromosome X. The Repbase collection contains rabbit betaretroviral LTR sequences designated ERV2-2_OC-LTR; these sequences did not contain regions homologous to either the BERVK1 primer (5′-TATCAGGCCTCTCCGCATG-3′ or its reverse complement) or the BERV-beta3 primer (5′-GGACCGTCTCCTTCAAGGC-3′ or its reverse complement). Nevertheless, homologous fragments were found in several lncRNAs and protein-coding genes. BERVK1-related fragments were identified in transcripts of SH2D4B (Chr15: 71,324,418–71,438,499), which encodes an SH2-domain adaptor protein; RNF208 (Chr1: 237,502,326–237,504,233), encoding an E3 ubiquitin ligase; and RFC5 (Chr21: 14,855,106–14,865,686), encoding replication factor C subunit 5.

BERV-beta3-related fragments were found in transcripts of CFAP57 (Chr7: 129,769,715–129,850,806). As in B. taurus, rabbit CFAP57 retains genetic linkage with genes including the ortholog of bovine C3H1orf210, TMEM125, and EBNA1BP2, a ribosome-biogenesis factor originally identified as an Epstein–Barr virus nuclear antigen 1-binding protein. Additional BERV-beta3-related sequences were identified in transcripts of RNF31 (Chr12: 51,038,560–51,051,473), the E3 ubiquitin ligase component of LUBAC, and COBL (Chr16: 55,540,709–55,725,661), a regulator of actin dynamics.

Sequences homologous to the BERVK1 and BERV-beta3 LTR fragments were also identified in the genome of the Asian honeybee, Apis cerana. A BERVK1-related sequence was found in a transcript of GUK-holder (GUKH; linkage group LG1: 10,419,444–10,521,427), a gene involved in cytoskeletal organization and cell polarity. BERV-beta3-related sequences were found in transcripts of sidekick (SDK; LG3: 10,241,397–10,283,451), a cell-adhesion gene [43]; a zinc-finger protein gene (ZFP568); DOLK, encoding dolichol kinase, an enzyme involved in protein glycosylation; the focal-adhesion protein gene BY; the nuclear hormone receptor Hr39/FTZ-F1 beta; the Krüppel-like factor luna (LOC108003801); an ortholog of the intracellular cholesterol transporter NPC; and an α1-tubulin gene (LOC107994074).

Overall, between the two betaretroviral LTR fragments, the BERVK1 fragment was substantially more species-specific than the BERV-beta3 fragment. In all three species, sequences homologous to these fragments were found in transcripts of genes whose products are primarily involved in transport, cell adhesion, cytoskeletal remodeling, protein ubiquitination, and transcriptional regulation, namely in fundamental cellular processes shared among multicellular organisms. These observations suggest that regulatory factors controlling their activity may also share features, universal across different species.

Noncanonical PU.1 binding motif. One such evolutionarily ancient group of transcription factors (TFs) is the ETS-domain family [44], whose core DNA-binding motif contains the GGAA/T sequence [45]. The noncanonical GGAGA motif generally shows weak ETS binding but, in certain sequence contexts, can bind with high affinity to the ETS-family transcription factor PU.1 [46]. The GGAGA sequence is present in both 5′-CATGCGGAGAGGCCTGATA-3′ (BERVK1 LTR, minus strand) and 5′-GCCTTGAAGGAGACGGTCC-3′ (BERV-beta3 LTR, minus strand). Because these fragments, which contain potential transcription factor-binding sites, occur within protein-coding sequences, their conservation may also be influenced by structural constraints on the encoded proteins, particularly when the fragments contribute to functionally important protein domains.

The idea that protein-coding sequences can have dual functions—encoding both amino acids and regulatory motifs, or “dual layers” of the genetic code—was proposed more than a decade ago [47] and continues to be supported by recent studies [48]. To examine whether such dual coding could contribute to the conservation of BERVK1 and BERV-beta3 LTR fragments across species, we determined whether these sequences contribute to conserved protein domains. For several genes, this was indeed the case.

In Bos taurus CFAP57 (XP_024845953, 1,305 aa), the BERV-beta3 LTR fragment encodes Lys-Glu-Lys-Val (KEKV), a sequence consistent with a coiled-coil region containing alternating hydrophilic and hydrophobic residues. The fragment lies within an evolutionarily conserved DNA-interacting domain annotated as an SMC chromosome-segregation protein domain (aa 820–1166; “SMC_prok_B,” “chromosome segregation protein SMC, common bacterial type”; TIGR02168; CDD: 274008). In Oryctolagus cuniculus CFAP57 (XP_002715219.2, 1,249 aa), the same fragment also encodes KEKV and occurs at the same position as in B. taurus.

In the E3 ubiquitin-protein ligase RNF31, the BERV-beta3 LTR fragment in B. taurus (NP_001422102, 1,070 aa) encodes Ala-Leu-Lys-Glu-Lys (ALKEK) and lies within the linker region of RNF31/HOIP. This region is important for autoinhibition and its release during formation of linear polyubiquitin chains required for activation of NF-κB, a major regulator of inflammation and cell survival [49,50]. The corresponding sequence in O. cuniculus RNF31 (XP_051678665, 1,067 aa) shows the same pattern as in B. taurus.

In Apis cerana Sidekick (SDK; XP_016912836, 2,220 aa), the BERV-beta3 LTR fragment encodes Glu-Gly-Glu-Gly (EGEG) within the extracellular domain responsible for binding to the plasma membrane of a neighboring cell and thus contributing to the architecture of cell–cell adhesion [43]. In Apis cerana Rfp 568 (LOC107994489, ID: XM_017051451.3, XP_016906940), the BERV-beta3 LTR fragment also encodes Glu-Gly-Glu-Gly (EGEG) and is located in the linker between a zinc-finger-associated domain (ZAD) and the first zinc finger (ZnF). ZAD-containing zinc-finger proteins can participate in chromatin organization and transposon regulation; for example, the ZAD protein Kipferl directs the heterochromatin protein Rhino to piRNA clusters in Drosophila [51].

Together, these examples suggest that the evolutionary importance of the proteins containing these sequences may contribute substantially to the conservation of LTR fragments derived from B. taurus-specific endogenous retroviruses, despite the different functions of these proteins. Conservation may also be influenced by the presence of GGAGA, a noncanonical PU.1-binding motif, within these fragments. PU.1 can compete with nucleosomes promoting local chromatin accessibility [46].

Host mechanisms suppressing TE expression. Another potential source of LTR fragments conservation is the host defense systems that suppress endogenous retrovirus expression. These include small noncoding RNAs, such as PIWI-interacting RNAs (piRNAs; 24–31 nt) and microRNAs (miRNAs; approximately 22 nt), as well as proteins containing a KRAB (Krüppel-associated box) domain. The KRAB domain recruits proteins involved in chromatin modification and transcriptional repression. It binds TRIM28 (tripartite motif-containing protein 28), which serves as a scaffold for additional repressive chromatin-modifying factors, including the H3K9 methyltransferase SETDB1, the histone deacetylase-containing NuRD complex, and heterochromatin protein 1 (HP1).

To examine this possibility, we searched these host regulatory systems for sequences homologous to the LTR fragments of the endogenous betaretroviruses expressed in the Bos taurus placenta.

The PIWI system consists of PIWI-family proteins and their associated small noncoding RNAs, piRNAs, which guide PIWI proteins to complementary RNA or DNA sequences and mediate transcriptional and post-transcriptional silencing, including repression of TE expression and mobilization. These mechanisms are particularly important in germ cells [52]. Databases of piRNA sequences are now available, and piRNAs have been characterized in B. taurus ovaries and testes [53]. Among these, two variants overlapping the BERVK1 LTR fragment were identified: 5′-TATCAGGCCTCTCCGCATGACTT-3′ in ovarian samples (GSM1584497) and 5′-ACAAAGGACATGCGGAGAGGCCTGATATGC-3′ in testes (GSM1584507), in which nucleotides 9–27 are fully complementary to the BERVK1 fragment.

Individual piRNAs containing several nucleotide substitutions (2–5 substitutions) relative to these sequences were found in piRNA-cluster databases not only for cattle but also for rabbit and bee [54]. The same was observed for the BERV-beta3 LTR fragment. Related piRNAs were detected in bovine oocytes at the germinal-vesicle stage (GV oocyte, GSM1584496; 5′-TAGCCTTGAAGGAGAAGGTCCTCGGA-3′, with one C/A substitution at position 13) and at metaphase II (MII oocyte, GSM1584501; 5′-TAGCCTTGAAGGAGAAGGTCCTTGGG-3′) [53]. These results suggest that incorporation of the LTR fragments into piRNAs capable of restricting TE transcription and mobilization may contribute to their conservation across species.

Other small noncoding RNAs, miRNAs, regulate gene expression post-transcriptionally. miRNA sequences also occur within the LTRs of some retroviruses; for example, bovine leukemia virus contains a cluster of five miRNAs [55]. Extensive complementarity between a miRNA and its target can promote target RNA cleavage, whereas partial complementarity commonly results in translational repression and/or mRNA destabilization. miRNA orthologs are often evolutionarily conserved among vertebrate species [56], and miRNAs have been characterized in detail in five domesticated mammalian species [57].

Several B. taurus miRNA sequences showed complete or partial homology to the LTR fragments examined here: three to the direct BERVK1 LTR sequence and four to its complementary sequence. Regions homologous to 5′-UAUCAGGCCUCUCCGCAUG-3′ were found in the bta-mir-2285cr-1 precursor on chromosome 10 (Chr10: 37,152,687–37,152,748), which has been associated with lipid metabolism [58]; the bta-mir-7691 precursor on chromosome 7 (Chr7: 81,526,133–81,526,199), located within an intron of SSBP2 (single-stranded DNA-binding protein 2), a gene involved in maintaining genome stability [59], with a direct repeat at positions 52–65 and a complementary sequence at positions 3–20; and the bta-mir-324 precursor on chromosome 19 (Chr19: 26,938,003–26,938,091), which has been associated with lipid-related signaling [60].

Four miRNAs showed homology to the complementary BERVK1 sequence, 5′-CAUGCGGAGAGGCCUGAUA-3′. These included the bta-mir-2375 precursor on chromosome 23 (Chr23: 3,280,212–3,280,285), which has been associated with immune responses and antiviral defense [61,62]; bta-mir-2887-1 on chromosome 2 (Chr2: 121,575,106–121,575,162), located within an intron of DST (dystonin), a member of the plakin protein family that links different components of the cytoskeleton and also contributes indirectly to vesicular trafficking [63], with additional copies on chromosomes 18 (bta-mir-2887-1; Chr18: 59,227,371–59,227,427) and 27 (bta-mir-2887-2; Chr27: 6,217,555–6,217,611); and bta-mir-2285aq-1 on chromosome 22 (Chr22: 17,212,965–17,213,019), which has been associated with lipid metabolism [64].

Sequences homologous to the direct BERV-beta3 LTR fragment, 5′-GGACCGUCUCCUUCAAGGC-3′, were identified in two miRNA precursors. The bta-mir-2285ah precursor (Chr7: 65,993,836–65,993,914), which has been implicated in regulatory networks controlling the corpus luteum [65], is located within an intron of LARP1, encoding La-related protein 1, a regulator of translation that interacts with ribosomes [66]. The second was the bta-mir-155 precursor on chromosome 1 (Chr1: 10,848,650–10,848,712). miR-155 is evolutionarily conserved and is one of the best-characterized miRNAs involved in immune regulation [67]. In humans, miR-155 has numerous experimentally supported immune-related target genes [68,69].

The complementary sequence, 5′-GCCUUGAAGGAGACGGUCC-3′, was also identified in two miRNA precursors. One was bta-mir-494 (Chr21: 65,920,313–65,920,397), which has numerous orthologs/paralogs and participates in PI3K–Akt signaling associated with coordinated development of bovine cumulus cells and oocytes [70]. It forms a cluster on chromosome 21 with four other miRNAs (MIR329B, MIR329A, MIR1193, and MIR543). The second, bta-mir-2474 (Chr8: 62,142,389–62,142,447), has been associated with adipogenesis and forms a cluster with MIR2472 and MIR2473 [71]. This region overlaps SHB, encoding an SH2-domain-containing adaptor protein involved in regulation of signaling pathways [72]. Together, these results indicate that conserved LTR fragments of the B. taurus placenta-expressed betaretroviruses contribute to miRNA precursor sequences.

In Oryctolagus cuniculus, despite the large number of miRNAs and their orthologs represented in current databases, we did not identify miRNA precursors containing regions homologous to the bovine LTR fragments. In contrast, two such precursors were identified in Apis melifera genome, sister species for Apis cerana, containing sequences homologous to the direct and complementary BERVK1 fragments: miR-3722 [73] and miR-6066 (miRBase: MI0020338; LG: 7,148,947–7,149,024) [74].

Thus, the association between small noncoding RNAs involved in TE silencing and the LTR fragments of B. taurus placenta-expressed betaretroviruses differs among species: homologous sequences were found among piRNAs in all three species examined, whereas no corresponding miRNAs have been identified to date in O. cuniculus.

Substantial differences were also observed between the nucleotide sequences of fusogenic env genes in mammals with different placental types (B. taurus, cotyledonary; O. cuniculus, hemochorial, as in primates). CENSOR analysis using Repbase identified no fragments of other transposable elements within the BERVK1 env sequence deposited in GenBank (AB751366.1) when compared with the corresponding consensus sequence in Repbase [23]. BERVK1 env and the nonfusogenic BERV-beta3 env differed substantially. In particular, BERV-beta3 env contained a short fragment at positions 120–169 showing 80% identity to mouse mammary tumor virus (MMTV), as well as a fragment at positions 557–607 showing 86% identity to the Gypsy-22_CQ-I LTR described in the genome of the southern house mosquito, Culex quinquefasciatus, a vector of several mammalian viruses [75].

The same analysis of the fusogenic rabbit syncytin-Ory1 sequence (GenBank: GU196371.1) identified fragments homologous to transposable elements from other animal species. Most notably, a 1,157-bp region (positions 523–1679) showed 68% identity to the proviral DNA of the betaretrovirus MacERV4_int from the common marmoset. Thus, the rabbit fusogenic gene was substantially more similar to a primate env sequence than to bovine env sequences: ORY1 showed 48% identity to BERVK1 env. It was also more similar to the primate sequence than the two bovine placenta-expressed env genes, BERVK1 and BERV-beta3, were to each other. A short 35-bp region of ORY1 (positions 383–418) showed 86% identity to a polinton DNA transposon described in Copidosoma floridanum [76]. C. floridanum is a polyembryonic parasitoid wasp of the family Encyrtidae (Hymenoptera), in which a single egg can give rise to more than 2,000 genetically identical offspring. Compared with BERVK1 env, BERV-beta3 env and ORY1 are less conserved and illustrate their independent evolutionary histories.

Of particular interest is the localization of the BERVK1 LTR fragment within the coding sequence of a KRAB-domain-containing protein. KRAB zinc-finger proteins are major regulators of endogenous retrovirus transcription in mammals [77]. The KRAB domain itself does not bind DNA; DNA recognition is mediated by adjacent zinc-finger domains. KRAB interacts with TRIM28, which recruits chromatin-modifying proteins that promote transcriptional repression [78]. In B. taurus, all isoforms of the KRAB-containing transcriptional regulator ZFP37 (NC_037335.1; Chr8: 102,210,192–102,305,553) contain the complete 5′-CATGCGGAGAGGCCTGATA-3′ sequence corresponding to the BERVK1 LTR fragment within the coding strand. ZFP37 belongs to the C2H2 (Cys2-His2) zinc-finger protein family and is characterized by an N-terminal KRAB domain and 12 tandem zinc fingers at its C terminus. Between these regions is a basic histone-H1-like linker capable of binding double-stranded DNA [38]. The presence of the BERVK1 LTR fragment within the ZFP37 coding sequence suggests that PU.1, an ETS-family transcription factor capable of recognizing the noncanonical GGAGA motif, may participate in regulation of a transcriptional repressor of TEs, suggesting potential interplay between transcriptional regulators with opposing effects.

4. Discussion

The conservation of short nucleotide sequences can reflect several different genomic factors, including the functional importance of the protein domains in which they occur and the contribution of regulatory motifs to the balance between transcriptional activation and repression, both within individual genomes and across species. The coexistence of these factors complicates interpretation of the functional significance of conserved sequences but also makes them useful markers for identifying different components of regulatory networks underlying phenotypic variation.

Identifying the different factors that may contribute to sequence conservation is also important because the same conserved sequence may be incorporated into different regulatory contexts in different species. Sakurai [79] described ERV contributions to placental development in terms of “software” and “hardware.” The “software” comprises regulatory elements involved in maternal–fetal immune and metabolic interactions, including interferon signaling, whereas the “hardware” comprises structural components, particularly fusogenic genes involved in placental tissue architecture.

The conserved BERVK1 and BERV-beta3 LTR sequences examined here illustrate that these components need not evolve in parallel. In cattle and rabbit, identical sequences occur within the coding region of CFAP57, a potentially “hardware” gene involved in the organization of cilia and flagella. At the same time, these species differ substantially in the nucleotide sequences of the fusogenic BERVK1 env and syncytin-Ory1 genes and in placental architecture: cattle have a synepitheliochorial placenta, whereas rabbits have a hemochorial placenta.

Importantly, the high sequence homology between env ORY1 and env genes of primates with hemochorial placentas suggests that the markedly different fusogenic activity of env genes in ruminants and primates may reflect differences in the relative contributions of cell fusion and other forms of intercellular interaction to placental development. The latter may include the formation of extracellular vesicles by trophoblast cells, as described in Bos taurus. [40].

Sequences homologous to the conserved BERVK1 and BERV-beta3 LTR fragments are distributed across multiple chromosomes in both cattle and rabbit. However, rabbit miRNAs currently represented in available databases do not contain these conserved sequences, unlike cattle. Given the evolutionary relationships between TEs and miRNAs, miRNA sequences may retain signatures of previous TE-derived sequence incorporation. One possible explanation for the occurrence of homologous BERVK1 and BERV-beta3 LTR fragments in distantly related species, including the Asian honey bee, is therefore the exaptation of fragments derived from different TEs. An alternative, not mutually exclusive mechanism is the acquisition by TEs of host regulatory motifs, such as the noncanonical PU.1-binding motif of the evolutionarily ancient ETS transcription factor family identified within both BERVK1 and BERV-beta3 LTR fragments.

Interestingly, across the species examined here, transcripts containing sequences homologous to the BERVK1 and BERV-beta3 LTR fragments encoded proteins involved predominantly in cytoskeletal remodeling, cell adhesion, and intercellular interactions. Homologous sequences were also identified in a number of lncRNAs and genes encoding transcriptional regulators. This functional distribution suggests a potential strategy for increasing the specificity of population-genetic analyses: a conserved ERV LTR sequence could be combined with fragments of genetically linked protein-coding genes, especially interesting for animal genotyping.