Recent studies have demonstrated that transposon-derived regulatory network elements are associated with phenotypic diversity (Fueyo et al., 2022; Du et al., 2024; Zhao et al., 2023). These studies are particularly important in the era of genome editing, as the structural and functional organization of potential gene editing targets requires careful study to assess possible negative side effects (Kosovsky et al., 2024).

Transposable elements (TEs) contain sequences necessary for their expression and transposition. These sequences can also regulate the host’s transcription, even when the TE itself has lost transposition function because of accumulated mutations. The mechanisms of host transcriptional regulation by TEs include providing enhancer and promoter sequences, targets of three dimensional chromatin organization, transcription factor (TF)-binding sites, and new regulatory elements including non-coding RNAs, the modification of nearby genes’ methylation states, and chromosomal structural modifications (Du et al., 2024). In mammalian genomes, there are gene blocks that are evolutionarily conserved and gene blocks that are frequently subjected to chromosomal rearrangements (homologous synteny blocks (HSBs) and evolutionary breakpoint regions (EBRs)). The differences between HSBs and EBRs are that the latter have high frequencies of TEs and microsatellite loci, more genes with a low number of base pairs, and they are frequently located in regions without topologically associated domains (TAD). HSBs are enriched with genes involved in anatomical morphogenesis and the development of the central nervous system (Damas et al., 2022).

Accumulated data suggest that when searching for genomic elements associated with phenotypic diversity, non-coding sequences with TEs make the greatest contribution, since they can regulate the expression profiles of adjacent genes (Du et al., 2024). TEs of different origin with the same TF-binding sequences are frequently located in syntenic regions. If one of the TEs carrying the corresponding TF-binding nucleotide motif is damaged, then its function can be performed by another TE variant carrying the same motif (Du et al., 2024). The density of TEs carrying TF-binding sites varies significantly between different genomic regions and depends on the localization of genes involved in different pathways (the lowest density is in the regions of genes associated with the early stages of embryonic development) (Damas et al., 2022; Nikitin et al., 2018; Nikitin et al., 2019).

Many studies are focusing on influencing phenotypic characteristics by editing TE elements (e.g., Buttler et al., 2023; Simpson and Chuong, 2023). However, the high evolutionary rate of TEs, as well as their diversity and multifunctionality, make this research avenue difficult.

Mammalian genomes are dominated by TEs and their derivatives that can comprise approximately half of the genome (Farmiloe et al., 2023). They are divided into two groups based on their mobilization intermediates. Class I elements, also known as retrotransposons, mobilize as an RNA intermediate and include autonomous retrotransposons (long interspersed nuclear elements (LINEs)and endogenous retroviruses (ERVs)), as well as non-autonomous short interspersed nuclear elements (SINEs). Class II elements, also known as DNA transposons, mobilize as a DNA intermediate. The diversity and abundance of TEs can be highly different, even for closely related species. One explanation for this pattern is that TEs with lower genomic abundance more strongly influence closely linked host genes and host activity to block TE transcription and insertion at other positions in the genome (Platt et al., 2018).

The difficulties in understanding TE effects on phenotypic diversity are not exclusively related with their extreme diversity and differences in copy numbers; there could also be species-specific TEs influencing the same phenotypic characteristic. For example, exhaustive studies support the hypothesis that different ERVs could be evolutionary driving forces for trophoblast cell fusion and species-specific placental structures (Sakurai et al., 2023). Thus, the same phenotypic characteristic can be formed with convergent participation of different TEs.

In order to elucidate the potential relationship between the function of protein-coding genes and the distribution of various TEs, we here perform a comparative analysis of TEs and the structural and functional features of the genes encoding melanophilin (MLPH, a transporter of melanosomes, in particular), myostatin (MSTN, an intracellular signaling system; a ligand of the transforming growth factor beta), leptin receptor (LEPR, an extracellular signaling system), and alpha-thalassemia (ATRX, involved in cell division mechanisms) in three mammalian species belonging to different orders: Oryctolagus cuniculus (domestic rabbit), Bos taurus (cattle), and Homo sapiens (human). Considering the associations between the functions of protein-coding genes and the distribution of TEs (Damas et al., 2022; Nikitin et al., 2018; Nikitin et al., 2019), one would expect that the lowest TE density would be in the ATRX gene, which is always located on chromosome X in mammals.

The analysis of TE distribution in genes and their flanking regions (1 million base pairs long, 1 Mb) includes LINEs, SINEs, endogenous retroviruses (ERVs, long terminal repeats; LTR in what follows), DNA transposons (DNA in what follows) and microsatellites. LINEs/L2 and SINEs/MIRs were considered separately because they are classified as evolutionarily “ancient” repeats compared to LINEs/L1 and SINEs/Alu, tRNA (Smalheiser and Torvik, 2005; Jurka et al., 2007). The subdivision of SINEs and LINEs into “ancient” and “young” transposons is associated with different frequencies of LINEs/L2 and SINEs/MIRs in mammals (monotremes, marsupials, and placentals). Estimates of the monotreme–theria divergence time range between 160 and 210 M years ago; platypus, placed with the echidnas into the Monotremata taxon, has mostly ancient LINEs/L2 and SINEs/MIRs (Warren et al., 2008). The 2.3 Gb platypus genome contains 1.9 and 2.75 million copies of LINE2 and MIR, respectively, compared to significantly lower numbers in other mammals (Warren et al., 2008). It has been shown that ancient LINE/L2s and SINE/MIRs are gradually replaced by young LINE and SINE repeats in mammalian genomes (Buc et al., 2018).

All four genes have distinct functions and belong to gene families of varying sizes and syntenic conservation. Only the ATRX gene is consistently located on the X chromosome across mammalian species. Comparative analysis reveals significant differences in the distribution of LINE, SINE, LTR, and DNA transposons across genes and species regardless of the localization and function of the genes considered.

For comparative analyses, O. cuniculus, GCF_009806435.1, B. taurus, GCF_002263795.3, and H. sapiens, GCF_000001405.40 (https://www.ncbi.nlm.nih.gov/, GenBank) genome assemblies were used. MLPH, MSTN, LEPR, and ATRX coordinates were uploaded from GenBank (https://www.ncbi.nlm.nih.gov/). H. sapiens genome assembly: GRCh38.p14 (GCF_000001405.40). B. taurus genome assembly: ARS-UCD1.3 (GCF_002263795.3). O. cuniculus genome assembly: mOryCun1.1 (GCF_964237555.1).

The gene neighborhoods were obtained with Gene Advanced Search Builder (https://www.ncbi.nlm.nih.gov/) using the “sort by chromosome” option and “chromosome & organism” query. Evolutionary conservation of genetic linkage analysis for genes in the neighborhoods of MLPH, MSTN, LEPR, and ATRX was performed in genomic segments approximately 2 × 10⁶ bp long, with the reference gene located in the middle of the segment. Non-annotated loci were excluded. The frequencies and distribution of dispersed and microsatellite repeats were found with RepeatMasker software tool (RepeatMasker.org, https://repeatmasker.org/cgi-bin/AnnotationRequest). Genomic TE annotations were obtained from RepeatMasker (https://repeatmasker.org/; H. sapiens—RepeatMasker version 406; B taurus—RepeatMasker version 405; O cuniculus—RepeatMasker version: 405), and correlation coefficients were estimated with https://www.statskingdom.com/correlation-calculator.html.

To estimate the evolutionary conservation of genetic linkage for the reference genes, the distribution of neighborhood genes was analyzed in 2 Mb segments for H. sapiens, B. taurus, and O. cuniculus genomes. Neogale vison (GCF_020171115.1) was also included in this analysis but was excluded from the analysis of dispersed repeats distribution because this species is not represented in the RepeatMasker database.

To estimate the evolutionary conservation of melanophilin’s genetic linkage, the distribution of neighborhood genes was analyzed in 2 Mb segments for the genomes of H. sapiens, B. taurus and O. cuniculus. The length (1 Mb) for 5′ and 3′ flanking regions was selected because, on average, the length of each chromatin loop is approximately 1 Mb long. Generally, the length of topologically associated regions varies from 100 kbp to 1 Mb (Cao et al., 2023; Glaser and Mundlos, 2022). We have previously examined the evolutionary conservation of genetic linkage for MSTN and LEPR genes in various species, including reptiles (Kosovsky et al., 2024).

The results of genetic linkage for MLPH are presented in Table 1. There is consistent genetic linkage between MLPH and COPS8, COL6A3, PRLH, and RAB17 genes. PRLH (enables neuropeptide hormone activity and prolactin-releasing peptide receptor binding activity) and RAB17 (a key gene regulating melanocytic filopodia formation in melanocytes in mast cells (Beaumont et al., 2011; Babina et al., 2022)) are genetically linked with MLPH even in birds (Bed’hom et al., 2012). All three genes participate in exosome formation and, eventually, in intercellular transport. Two other genes, genetically linked with MLPH, COL6A3 (encodes the alpha 3 chain of type VI collagen, a protein that plays a crucial role in the extracellular matrix), and COPS8 (one of the eight signalosome complex subunits that interact with ubiquitin ligase and is involved in histone modifications) are presumably involved in fundamental processes of intercellular interactions and the differential regulation of gene expression. It is known that E3 ligase triggers extensive gene expression reprogramming by changing global levels of H3K27me3 (a histone modification associated with gene repression) in plants (Wang et al., 2024), bees (Lowe et al., 2022), and in some cellular populations in tumorigenesis (Marine, 2012).

The MLPH, MSTN, LEPR, and ATRX genes of all three species have variable length, number of exons, and dispersed repeat frequencies (non-autonomous retrotransposons SINE, autonomous retrotransposons LINE, endogenous retroviruses LTR-ERV, DNA transposons, and tandem microsatellite repeats, STR) inside the genes and their flanking regions. The number of TEs inside the genes was significantly different: MSTN (extensively used for successful gene editing in livestock) had the least amount of TEs and LEPR and ATRX had the most. There were also significant differences in the ratio of “young” and “ancient” SINE and LINE repeats. For humans and cattle, the balance between SINE and LINE inside LEPR gene and in its flanking regions was shifted toward ancient repeats compared to MLPH and MSTN genes (Figures 1–6; Tables 4, 5).

In this context, it should be noted that only LEPR and ATRX have gene insertions inside: LEPR has the gene LEPROT (its sequence overlaps with the LEPR gene itself) and human ATRX has FABP5P15 pseudogene (fatty acid binding protein 5 pseudogene 15, GRCh38.p14, NC_000023.11), as well as long non-coding RNA (lncRNA). In addition, the high frequencies of young SINE and LINE in LEPR and ATRX suggest they are less insulated from DNA insertions than MLPH and MSTN (specifically MSTN, Figures 1–3). Moreover, LEPR belongs to an evolutionarily conserved gene cluster (Kosovsky et al., 2024), and ATRX in mammals is always located on the X chromosome. Thus, TEs can also target evolutionarily conserved gene blocks, not only the evolutionarily “fragile” chromosomal regions. However, the relatively high number of ancient (SINE/MIR + LINE/L2) versus young TEs (SINE/Alu and SINE/tRNA + LINE/L1) in the flanking regions of LEPR, compared to the MLPH and MSTN flanking regions in humans and cattle (Figure 6), could be associated with LEPR being in a genetically-linked evolutionarily conserved gene cluster.

Questions about the origin, evolutionary mode, and functions of dispersed repeats in the genome are still largely open, despite more than two decades of active research (Jordan et al., 2003). The origin of dispersed repeats is complex; for example, the SINE/Alu repeat originated from the 7SL RNA sequences, a non-coding RNA that is a key component of the signal recognition particle and is essential for protein secretion and translocation into the endoplasmic reticulum. Another example is SINE/tRNA, which originated from tRNA genes; some other repeats are recombination products and overlap, so their origin is untraceable (Sharma et al., 2024). Domestic rabbits have two families of dispersed repeats—OcuSINEA and OcuSINEB—with their central parts originating from LINE and not tRNA, demonstrating that SINE can coopt LINE’s endonucleases (Yang et al., 2021). However, based on the similar frequencies of SINE/Alu in human and SINE/tRNA in cattle and rabbit genomes and the presence of similar regulatory elements in both repeats, it is reasonable to assume that these repeats could have some exchangeable functionalities. For example, SINE/Alu, other SINEs, and some LTR and DNA transposons (e.g., DNA/hATCharlie) contain the same motif for transcriptional regulatory factor CTCF (responsible for chromatin higher order structure). Even without the CTCF binding motif, the clusters of SINE/Alu and SINE/MIR influence local chromatin organization (Gunsalus et al., 2023). In addition, SINE/MIR, LTR, and SINE-VNTR-Alu form a “transposon code”, providing particular patterns of transcription factor binding sites for regulatory networks (Testori et al., 2012; Trizzino et al., 2017). It should be noted that transcription factors that bind heavily depend on the DNA methylation of their binding sites (Rimoldi et al., 2024), and the methylation itself is an epigenetic silencing mechanism of TEs transcription that has both species-specific and common characteristics in vertebrates (Carotti et al., 2023; Ilık et al., 2024). In addition, SINE, LINE, ERV, and DNA transposons are the major sources of regulatory network dynamics, regulating transcription factors, providing promoter, enhancer, and silencer sequences, histone and microsatellite codes, methylation and imprinting, and different non-coding RNAs (Mandal, 2024). For example, in different mammalian species, not only do all TE groups contain transcription factor-binding sites but also many genes encoding these factors contain different TE groups (Du et al., 2024).

In summary, the accumulated evidence suggests that all transposons, even from different classes, despite their structural differences and evolutionary origin, have common characteristics that are related to the mechanisms of their replication, transposition, genomic integration, involvement in 3D genome organization, and interactions with host genes. The interaction could be at the level of neighboring genes and could also be at varying distances. Thus, the common trait of all different TEs is their ability to participate in the host’s regulatory networks by providing a multitude of different regulatory mechanisms to them. Their integration and distribution are defined not by external factors but by instant need inside the genome organization that can be solved by any available source of regulatory element, at any evolutionary moment, thus optimizing survival and successful reproduction. In other words, TE repeat integration is defined not by group membership but by the availability of a TE with necessary regulatory elements. For example, some TEs play a supporting function in the interphase nucleus architecture (Jordan et al., 2003; Falk et al., 2019). This “internal” genomic adaptation was first defined by the architectural term “spandrel” by Gould and Lewontin (1979) and was later adopted by many others (Koonin, 2011). It can be expected that the accumulation of TEs and their turnover is a consequence of the non-adaptive accumulation and exchange of various “junk” genomic elements, functionally convergent in the presence of regulatory motifs (Koonin, 2016).

In the area of finding the most effective markers associated with phenotypic diversity, this approach could shift attention from finding SNPs in causative loci or finding specific genomic elements toward estimating polymorphism for the specific members of regulatory networks.