Widespread HCD-tRNA derived SINEs in bivalves rely on multiple LINE partners and accumulate in genic regions

Abstract

Background

Short interspersed nuclear elements (SINEs) are non-autonomous non-LTR retrotransposons widespread across eukaryotes. They exist both as lineage-specific, fast-evolving elements and as ubiquitous superfamilies characterized by highly conserved domains (HCD). Several of these superfamilies have been described in bivalves, however their overall distribution and impact on host genome evolution are still unknown due to the extreme scarcity of transposon libraries for the clade. In this study, we examined more than 40 bivalve genomes to uncover the distribution of HCD-tRNA-related SINEs, discover novel SINE-LINE partnerships, and understand their possible role in shaping bivalve genome evolution.

Results

We found that bivalve HCD SINEs have an ancient origin, and they can rely on at least four different LINE clades. According to a “mosaic” evolutionary scenario, multiple LINE partner can promote the amplification of the same HCD SINE superfamilies while homologues LINE-derived tails are present between different superfamilies. Multiple SINEs were found to be highly similar between phylogenetically related species but separated by extremely long evolutionary timescales, up to ~ 400 million years. Studying their genomic distribution in a subset of five species, we observed different patterns of SINE enrichment in various genomic compartments as well as differences in the tendency of SINEs to form tandem-like and palindromic structures also within intronic sequences. Despite these differences, we observed that SINEs, especially older ones, tend to accumulate preferentially within genes, or in their close proximity, consistently with a model of survival bias for less harmful, short non-coding transposons in euchromatic genomic regions.

Conclusion

Here we conducted a wide characterization of tRNA-related SINEs in bivalves revealing their taxonomic distribution and LINE partnerships across the clade. Moreover, through the study of their genomic distribution in five species, we highlighted commonalities and differences with other previously studied eukaryotes, thus extending our understanding of SINE evolution across the tree of life.

Background

Bivalves (Class Bivalvia) are a rich and widespread clade of aquatic-only molluscs that diversified back in early Cambrian, more than 500 Mya [1]. This class include multiple economically and ecologically important species. For example, they have colonized freshwater environments [2] and deep-sea vents multiple times during their evolutionary history [3]. They can be useful bioindicators for marine pollutants [4] and they can represent biological models for the study of adaptations to climate change [5], innate immunity [6], sex determination [7], longevity [8, 9] and mitochondrial biology [10]. Moreover, they are characterized by peculiar genomic features that have been hypothesized to be linked to transposable elements (TEs) activity, such as transmissible cancers [11], high levels of hemizygosity [12] and gene presence-absence variation [13].

Their important role as promising model system for addressing both general biology and human health questions, together with the increased cost-efficient accessibility of third-generation sequencing technologies, has led to a major increase in their genomic resources in recent years [14]. This has opened the possibility to explore, in a broader context, also usually neglected genomic components and their evolutionary dynamics, such as TEs and other repetitive sequences which constitute a high proportion of bivalve genomes [15].

Repetitive DNA elements usually replicate in a selfish manner, independently from host genome replication, with variety of effects on the host fitness ranging from neutral to deleterious. However, multiple cases of co-option in novel functions have been described in literature [16]. Furthermore, their evolutionary trajectory can be influenced by the dynamic of the host population, which in turn may be affected by the changes in TE activity [17]. Therefore, our understanding of TE distribution and evolution across the tree of life represent an important step in a broader understanding of evolution of living forms.

Short Interspersed Nuclear Elements (SINEs) are a sub-class of non-autonomous, non-LTR retrotransposons that depend on the protein machinery of their autonomous counterpart LINEs (Long Interspersed Nuclear Elements) to reintegrate into the genome after their transcription by RNA polymerase III (Pol III) [18,19,20]. Moreover, while many non-autonomous elements usually originate from their autonomous counterparts though sequence decay or internal deletion, such as Miniature Inverted-repeat Transposable Elements (MITEs [21],) or Short Internally Deleted Elements (SIDEs [22]), SINEs emergence is only partially dependent from their LINE partners [19]. Their canonical structure comprises a head, a body, and a tail region [20]. The head can originate from one of the three RNA type synthetized by the RNA Pol III, namely tRNAs, 5S rRNAs, or 7SL RNAs, and contain its promoter region. Even if elements originated from all three RNA types have been observed across a wide range of eukaryotes, tRNA-derived SINEs appear the most common [19]. The body, when present, contains a domain of unknown origin and function, which appears to be element-specific [20]. However, in some instances, SINEs may carry bodies with highly conserved domains (HCD) across distinct SINE lineages and hosted by distantly related species. Although the role of HCDs is still unclear, they have been useful for classifying SINEs at the superfamily level [23,24,25]. Finally, the 3’ tail region serves as recognition for the LINE-derived reverse transcriptase (RT) and it may terminate with tandem repeats or an A-rich segment [20]. The SINE-LINE partnership can be specific if a LINE derived segment—usually originated from the 3 ‘UTR LINE region—is required for RT recognition, or aspecific when homology is not necessary [19]. The modular structure of SINEs suggested a characteristic evolutionary model called “mosaic evolution” under which different SINE lineages can exchange their modules through recombination [26]. This feature could allow their long-term persistence under a strict vertical inheritance evolutionary scenario in different genomic context, for example after the extinction of the original LINE partner [24, 25].

A few analyses in bivalves already identified HCD SINEs belonging to the superfamilies Core [23, 27], V [28, 29], Meta [23], Deu [30], and MD [23], where the latter are composed by a dimerization of Meta and Deu domains. Despite covering a low percentage of bivalve genome [15], their emergence was traced back to their most recent common ancestor [23] similarly to what was hypothesized for the diversity of LINE clades [15]. Moreover, for some of these elements Nishihara et al. [23] and Matetovici et al. [31] identified also their putative autonomous partners: CR1, L2, and Nimb. However, these studies were limited by the reduced number of whole genome assemblies available at that time and by the lack of a comprehensive LINE reference library for bivalves.

Here, we leveraged the recent increase in bivalve genomic resources to comprehensively characterize HCD SINE diversity and richness across bivalve evolutionary history. The newly generated SINE library was used to screen for putative and previously unknown SINE-LINE partnerships, revealing that at least 4 different LINE lineages could act as RT donors in seven different SINE-LINE partnerships. Moreover, since SINEs can be important contributors to gene and genome evolution, we conducted a case study analysis on a subset of five species to investigate the possible impact of SINE in genome evolution. Our findings showed that gene-related genomic regions are enriched in SINEs—particularly in old copies—and that they can be organized in tandem-like and palindromic structures, potentially affecting epigenetic regulation. With respect to the recent TE survey conducted in bivalves [15], this analysis further increases and refine the SINEs dataset, along with that of partner LINEs, also providing a deeper analysis of their genomic distribution and organization in bivalve genomes.

Methods

Genomic dataset for de novo SINE prediction and manual curation

We selected 11 genomes from NCBI [32] and GigaDB (Supplementary Table 1) for the de novo mining and manual curation of SINE elements. SINE candidates were mined from each assembly using RepeatModeler2 [33] and SINE_Scan v1.1.1 [34, 35]. For each species, we merged SINE_Scan representative sequences with all TE consensus resulting from RepeatModeler2 and annotated as SINE by RepeatClassifier or by deepTE [36], which was run on “Unknown” elements to increase the chance of include as many SINEs candidate as possible. Candidates elements were then subject to a “Blast-Extend-Extract" process [37] by blasting back each element against its source genome (Blastn v2.6.0: qcov_hsp_perc 70, perc_identity 70; [38]), extracting the top 50 hits + 300 bp at both ends with bedtools v2.26.0 [39] and aligned with MAFFT v7.475 [40]. From each alignment, we built a novel consensus sequence using the online Advance Consensus Maker tool (https://www.hiv.lanl.gov/content/sequence/CONSENSUS/consensus.html). Boundaries of the elements were manually identified looking for the characteristic decay of the alignment towards terminal regions and the consensus sequences was curated implementing a majority rule approach, following the guidelines of Goubert et al. [37] and Peona et al., [41]. To confirm a candidate as a tRNA-related SINE we required: (a) the presence of a microsatellite or a poly-A region at the 3’ end; (b) the presence of a tRNA- related region on the 5’ end predicted by tRNAscan-SE [42], through homology searches on the GtRNAdb 2.0 [43], http://gtrnadb.ucsc.edu) or manually looking for RNA Pol III A and B boxes, and (c) a length between 200 and 700 nucleotides. The presence of characteristic TSDs between 6 and 18 bps was manually checked, although it was not required to confirm a candidate as a SINE.

SINE-LINE partnerships

To identify partnerships of SINEs with their autonomous LINE counterparts, we queried all confirmed SINEs (blastn: word size = 7, gap opening penality = 2, gap extension penality = 2, Match score = 2, Mismatch score = -3, evalue = 0.01) against a bivalve-specific library of LINE elements [15]. All positive hits were manually checked to confirm that the homolog region would fall at the 3’ tail of the SINEs and within the 3’ UTR of the LINE element.

Co-evolutionary dynamics of SINEs and their LINE counterparts were studied in the genomes in which we found evidence of SINE-LINE partnerships. For this purpose, we first selected all assemblies for which we identified homology between the tail region of a confirmed SINE mined from the same assembly and any LINE 3’ UTR region. We then attempted to build a species-specific representative sequence of the LINE counterpart by blasting the original LINE element against the genome with decreasing thresholds in terms of identity and required alignment length. In this way, we obtained a novel consensus sequence following a Blast-Extend-Extract process, as previously described (See “Genomic dataset for de novo SINE prediction and manual curation”). When no homology was identified with a blastn search, we performed more sensitive tblastn searches (E-value 1E-05) of the amino acid translation of ORF2. Species-specific LINE consensus sequences were then checked for conservation of the homologous region with the species-specific SINE tail region. We assessed the completeness of the LINE consensus sequence with TE-aid [37] (min ORF length = 300 aa). For confirmed partnerships, we used all species-specific SINE-LINE partner pairs as custom libraries for RepeatMasker in sensitive mode against the source assembly. TE landscapes, describing the divergence of each TE copy from its consensus sequence in terms of percentage of Kimura distance after CpG corrections, were calculated using the calcdivFromAlign.perl script provided with the RepeatMasker installation. Concurrent activity between SINE-LINE partners was further tested for each species with Spearman's rank correlation tests between accumulation profiles (i.e., number of base pairs occupied in each bin of CpG corrected Kimura divergence) of the two elements.

Superfamily and family level classification of confirmed SINEs

For HCD SINEs classification we follow the superfamily classification scheme of Nishihara et al., [23] based on the presence of characteristic central domains previously identified in bivalve genomes (Meta, V, Deu, Core). We started provisionally annotating each element using the RepeatClassifier utility from the RepeatModeler package. Elements that should share the same central domain were aligned using MAFFT and we then manually checked for the presence of the characteristic domain.

To obtain species-specific SINE families, we clustered all HCD SINEs mined from the same source genome following the 80–80 rule [18] using cd-hit-est v4.7 (-G 0 -c 0.8 -aS 0.8 -t 1; [44]). Clusters were further refined into families following the definition from SINE Base [20], where a SINE family is described as “a set of elements sharing the same modules in the same order, excluding the tail region”, where the tail represent the putative LINE-derived region + the poly-A/microsatellite. To achieve this, we ensured that each cluster contains only elements with the same modules; when this criterion was not met the original cluster was split into different families. BLASTn was used to identify homologous regions within SINE tails between different families, with a permissive e-value of 0.05 but considering only alignments >  = 50 bp. Homologous relationships were visualized as a network with Cytoscape v3.10.2 [45].

Finally, all families were merged into a multi-species SINE library together with 19 elements previously described and deposited in RepBase (Supplementary Table 2) and re-clustered with the same methods describe above to identify families shared by multiple species (cross-species families).

Copy number estimation of tRNA-related SINEs across bivalve diversity

To estimate the distribution of the four SINE superfamilies across bivalve diversity we downloaded additional 20 genomes from NCBI (Supplementary Table 1) and performed homologous searches with blastn of all species-specific SINE families as queries (E-value 1E-05). To avoid crossmatch with tRNA donors and LINE homologous regions, we excluded hits shorter than 150 bp (i.e., approximately shorter than the 50% of the entire SINE length). After this step, we merged overlapping hits resulting from different families of the same superfamily using bedtools merge and counted the number of occurrences of each superfamily in each genome. A maximum of 150 random copies belonging to the superfamilies V, Meta, and Core were extracted from each genome and aligned using MAFFT in auto mode. TrimAl v1.4 [46] was used to remove gap positions (–gappyout mode) and spurious sequences from the alignment (-resoverlap 0.50 -seqoverlap 55). We inferred a Maximum Likelihood tree via FastTree v2.1.10 [47] using a GTR + Gamma model. After adding the gastropod Biomphalaria glabrata (GCF_947242115.1) to the genomic dataset as outgroup, we inferred phylogenetic relationships between the analysed species by extracting and concatenating the 331 complete and single-copy BUSCO genes (Metazoa odb10 dataset; [48]) present in at least 90% of the species. We provided the partitioned supermatrix to IQ-TREE2 [49] with ModelFinder [50] and 1,000 ultrafast bootstrap replicates [51] to respectively find the best-fit partitioning scheme and evolutionary model and assess nodal support of the inferred maximum likelihood tree. Divergence time estimation was performed with LSD2 [52] within IQ-TREE2 using the diversification of analysed bivalves and the crown node of all bivalve orders (except for Adapendonta and Cardiida) as calibration points based on the median divergence time reported in the TimeTree5 online database [53] (last accessed in July 2024).

Additionally, we estimated the taxonomic distribution of species-specific SINE families by blasting each SINE against our extended bivalve genomic dataset. We considered a family present in a genome when we could identify at least 10 hits with a query sequence identity and coverage of at least 80%.

Genomic distribution of SINEs and prediction of tandem-like SINE structures

All cross-species HCD SINE families were used as input library for RepeatMasker v4.1.0 [54] in sensitive mode (-s) to study their genomic occurrence in five species with available gene annotation. Specifically, for Crassostrea gigas (Ostreida), Mytilus californianus (Mytilida), and Pecten maximus (Pectinida) the RefSeq gene annotation was downloaded from NCBI repository, while for Ruditapes philippinarum (Venerida) and Scapharca broughtonii (Arcida) they were retrieved from Xu et al., [55] and Bai et al., [56], respectively. We considered five different features: exons, introns, annotated UTRs, 2,500 bp flanking the genes, and all other intergenic sequences (thus excluding 2,500 bp gene flanking regions). For each feature, we counted the number of intersections with SINE insertions with Bedtools intersect. Over- and under-representation of SINEs in each feature was tested by constructing—with Bedtools shuffle—null distributions from 1,000 random reshuffling iterations of all annotated SINE insertions across the genome (excluding genomic gaps). At each iteration, the number of intersections between each feature and the random intervals were counted. The observed number of intersections of SINEs in each feature was then compared to the null expectation. To directly test the hypothesis of SINE preferential accumulation in 2,500 bp gene flanking regions compared to all other intergenic regions, we split both features into intervals with a window of 500 bp with Bedtools windows and selected 10,000 random intervals for 100 iterations. We then counted the number of overlaps with SINE annotations at each iteration as previously described. Results for intergenic and gene-flanking genomic regions were statistically compared using t-test. Taking advantage of the high-quality repeat annotation of C. gigas, whose repeatome is almost completely characterized [15], we also studied the accumulation patterns of LINEs across the same genomic intervals. Briefly, all LINEs from C. gigas available in RepBase, as well as those identified by Martelossi et al. [15], were combined, redundancy reduced with cd-hit-est following the 80–80 rule and used to annotate the genome with RepeatMasker. Overlaps between LINE insertions and genomic features were counted and statistically tested as previously performed for SINEs. Finally, we additionally hypothesize that gene- genomic regions, here defined as UTRs + exons + introns + 2,500 bp gene flanking regions, are characterized by older SINEs copies compared to intergenic ones. To test this, all gene- genomic regions were merged, and we calculated, for both gene-related and intergenic genomic regions, the percentage of Jukes-Cantor (JC) distance of each SINE copy to its consensus as a proxy for the time of insertions. Distributions were then tested with t-test. The same analyses were also performed for LINE insertions in C. gigas.

The same five genomes were scanned to identify presence of tandem SINE arrays. For this purpose, we only kept high-scoring SINEs, i.e. RepeatMasker annotated insertions with a score higher than 400 and with a length of at least 150 bp. This was necessary to remove possible misannotations such as host tRNA genes. We consider tandem-like SINE structures when multiple elements coming from the same family were detected one after the other.

Collection of seed alignments for DFAM submission

SINE family consensus sequences were used to build up seed alignments for DFAM submission [57]. For this purpose, we used the generateSeedAlignments.pl script provided with RepeatModeler installation with the flags –taxon, specifying the species name as reported in NCBI taxonomy, and –assemblyID followed by the NCBI accession number of the assembly. Resulting Stockholm files were submitted to DFAM. Consensus sequences can be found in Supplementary Data S1.

Results

An improved HCD tRNA-related SINE library for bivalves

By combining RepeatModeler, SINE_Scan, and homology searches, we identified 201 SINEs across the 25 selected bivalve genomes analyzed for the initial screening of SINE candidates (Supplementary Table 2). All confirmed elements exhibited signatures of tRNA-related origin based on tRNA prediction analyses, homology searches against the GtRNAdb and/or manual identification of putative A and B boxes, the typical RNA polIII promoter (Supplementary Table 2, Fig. 1). The tail region of candidate SINEs was also checked for the presence of microsatellites, and we successfully identified characteristic TSDs with sizes ranging between 6–18 bp for 181 (90%) of these elements (Supplementary Table 2). Comparative analyses using the domains described in Nishihara et al., [23] allowed us to subdivide these elements into the five known HCD superfamilies: Meta, V, MD, Core, and Deu. Specifically, we classified 31 elements as Core, 16 as Deu, 34 as MD, 40 as Meta, and 53 as V. Additionally, we found 27 other SINEs without clear homology to the aforementioned domains which we simply classified as tRNA-related SINEs. Within the five known HCD superfamilies, we identified 10 putative different tRNA donors, which are also shared between different superfamilies (Supplementary Table 2).

Fig. 1

Schematic representation of identified HCD tRNA-derived SINEs in bivalves. For each superfamily we reported the tRNA-related heads identified with tRNA-Scan SE and the putative LINE donors

These 201 elements were firstly clustered in species-specific families (Supplementary Table 2) and then, along with 19 publicly available bivalve SINEs into 76 distinct cross-species families based on criteria of reciprocal homology and order of SINE modules (Supplementary Table 2). Among these families, 17 are composed of unknown SINEs, 14 of Core SINEs, eight of Deu SINEs, six of MD SINEs, eight of Meta SINEs, and 23 of V SINEs (Supplementary Table 2). Unknown families were excluded from following analyses. The presence of 22 families shared by multiple species (five Core, four Deu, one MD, five Meta, seven V) belonging to the same bivalve order (Supplementary Table 2) underlies the possible long-term conservation of HCD SINE. Some notable examples are the families Bpla_SINE-1_Meta (tRNA head: Ser) shared between Mytilinae and Bathymodiolinae, Tgra_SINE-7_Meta (tRNA head: Ser) shared between all Arcidae and Oden_SINE-1_CORE (tRNA head: Ser) shared between O. denselamellosa, C. gigas and S. glomerata.

Bivalves HCD SINEs depend on at least 4 different LINE lineages

Using curated LINE libraries previously obtained from mollusc genomes [15] together with all newly generated SINE consensus sequences, we searched for putative SINE-LINE partnerships. Our results highlights that at least four different LINE clades can match any of the SINE tails (Fig. 1; See Supplementary Table 2 for all recognized homologies). Homologies between SINE and LINE 3’ ends can be shared between different superfamilies and span between 35 and 61 bp with an identity ranging from 72% to 95% (Supplementary Fig. 1). Nishihara et al. [23] and Matetovici et al., [31] found similarities between tail regions of V and Core families with CR1 and L2 elements and between Meta SINEs and Nimb LINEs (LINE I superfamily). Here we found that not CR1, but CR1-Zenon elements, a LINE clade closely related to CR1 and widespread in bivalves but apparently poor in other molluscs [15], are likely responsible for the retro-transcription/reintegration of V, Core, Meta, and Deu families, while Nimb LINEs may promote V and Meta replication. It is interesting the case of the M.phylippinarum_91 CR1-Zenon LINE element, which show clear homology to both Meta, Deu, V and Core elements (Supplementary Fig. 1). Additionally we also found only one family coming from the Venerida Archivesica marissinica genome with a tail region highly similar to 3’ ends of CR1 elements (Amar_SINE-2_CORE).

The homology relationships of SINE tail regions strongly support previous results. Indeed, we identified 21 different clusters characterized by SINEs with homologous or partially homologous tail regions (Supplementary Fig. 2). Sixteen of these clusters contain only families without a clear LINE-derived region, whereas the others are characterized by at least one element with one of the previously described LINE homologies. Interestingly, while tails coming from the same SINE superfamily tend to be more closely associated within the similarity network (i.e., more similar), multiple cases of homologous or partially homologous tails can be found between elements coming from different HCD superfamilies. Notably, the two biggest clusters contain almost all SINE tails that we identified as Nimb and CR1-Zenon related.

To study the co-evolutionary dynamics between autonomous and non-autonomous elements, we first reconstructed, where possible, a full-length species-specific LINE counterpart. We managed to reconstruct 10 species-specific SINE-LINE partnerships across 9 different species. All reconstructed LINE partners resulted longer than 4,000 bp with all but one (Myes-2_LINE-I) possessing an ORF2 which encodes for an RT + EN domain. Among these, eight out of nine also present an ORF1 demonstrating the completeness of their structure (Supplementary Fig. 3). Repeat landscapes analyses revealed contrasting patterns between LINE and SINE activity (Supplementary Fig. 4). Indeed, despite a positive and significant correlation between all accumulation profiles (all p-values < 0.05, Supplementary Fig. 5) both visual inspection of repeat landscapes profiles as well as correlation analyses point to possible different evolutionary scenarios. Specifically, the accumulation patterns of BivaV-SINE2_CrGi#V, Medu_SINE-2#Meta, and Sbro_SINE-2#Meta resulted lower and less significantly correlated to their LINE counterparts (0.3 < Spearman’s rho < 0.44; 0.002 < p-values < 0.03; Supplementary Fig. 5) compared to other analysed partnerships. On the contrary the partnerships Cgig_SINE-10#CORE / Cgig-1_LINE#L2, BivaV-SINE1_MiYe#V / Myes-2_LINE#Nimb, and Amar_SINE-2#CORE / Amar-1_LINE#CR1 show both strong correlations (Spearman’s rho > 0.8) and overlapping activity profiles.

HCD tRNA-derived SINEs are widespread in bivalves and maintained activity after bivalve order diversification

To have a broader picture of SINEs HCD family and superfamily distribution across bivalve diversity, we added other 20 assemblies to our starting genomic dataset for a total of 45 analyzed species representative of 10 different bivalve orders (Supplementary Table 1). After inferring phylogenetic relationships and divergence time estimation based on 331 single copy BUSCO genes present in at least the 90% of the species, we used these genomes as database for homology searches using all previously confirmed SINEs as queries. Recovered phylogenetic relationships are largely coherent with published phylogenomic results (Fig. 2A and B; [58, 59]), identifying all analysed bivalve orders as monophyletic with high bootstrap values. The only exception is Gari tellinella (superfamily Psammobiidae, Gtell in Fig. 2A and B) which is recovered in sister relationship with analyzed Adapendonta species and separated from other Cardiida belonging to the Cardioidea superfamily, as already reported [58].

Fig. 2

Taxonomic Distribution of HCD superfamilies and families in bivalves. A Taxonomic distribution of the 5 known HCD SINE superfamilies and B number of shared families between each pair of analyzed species. Species name abbreviations refer to Supplementary Table 1. Phylogenetic relationships and divergence time estimation are represented at the top and bottom of the two panels. All nodes receive a bootstrap value ≥ 95 except those labelled in B

Concerning HCDs distribution, the Core superfamily was found across all members of the orders Pectinida, Ostreida (except for Pinna nobilis), Unionida, of the superfamily Cardioidea as well as in two Adapedonta (Solen grandis and Sinonovacula constricta), two Arcida species (Anadara kagoshimensis and Scapharca broughtonii) and in four Venerida genomes (Archivesica marissinica, Cyclina sinensis, Mactra quadrangularis, Ruditapes philippinarum) (Fig. 2A). No Core element could be identified in Mytilida, Myda and Lucinida representatives.

The Meta superfamily appears widespread across analyzed Arcida, Mytilida, Unionida (except for Margaritifera margaritifera), Adapedonta and in the three Venerida Saxidomus purpuratus, Spisula solida and Mercenaria mercenaria (Fig. 2A).

We identified the Deu superfamily in the majority of the Arcida, Cardioidea, and Venerida as well as in some Mytilida (Mytilus coruscus, Mytilus edulis, Mytilus californianus, and Mytilisepta virgata) and Ostreida (Ostrea denselamellosa, Saccostrea glomerata and P. nobilis) (Fig. 2A). Interesting, the deep-sea symbiotic clam A. marissinica hosts ~ 17 times more Deu elements than the second richest species S. glomerata (118,575 and 6,816, respectively) confirming an in increased of activity of specific TE groups in this lineage, possibly related to its colonization of hydrothermal vents [15, 60].

The MD superfamily appears as the least represented across bivalves, while the V superfamily resulted the most ubiquitous with elements identified across all species except for S. glomerata, Cangeria kusceri, and Fragum whitleyi, confirming what was previously found by Nishihara et al., [23] (Fig. 2A). Interesting both Meta and V superfamilies are present in similar high copy number across four out of the six Unionida species here analyzed (Hyriopsis cumingii, Unio delphinus, Megalonaias nervosa and P. streckersoni,from 27,000 to 101,712 copies). Similarly, also the LINE complement of M. nervosa and P. streckersoni appeared to be different compared to other bivalves in Martelossi et al. [15]. However, the apparent absence of SINE V and META amplification in the P. streckersoni sister species Venustaconcha ellipsiformis points to two possible different evolutionary scenarios: an amplification in their most recent common ancestor with subsequent genomic deletion in V. ellipsiformis or an independent amplification along different Unionida lineages. The fast increase in high quality Unionida genomic resources will possibly shed light on which process is driving the different TE landscape of this group.

Coherently with the results obtained from the cross-species family clustering (see Results: “An improved HCD tRNA-related SINE library for bivalves”), we found that bivalves belonging to the same order and those phylogenetically close are more likely to possess shared SINEs (Fig. 2B). For example, we found additional evidence for shared families between Mytilinae and Bathymodiolinae, which diverged around 324 Mya in our estimation and ~ 400 Mya according to Lee et al. [61] as well as between all Arcida (divergence time: ~ 178 Mya here and ~ 177 Mya in [62]). Interestingly, we also found shared families between Adapendonta, G. tellinella, Myda, and Venerida, which diverged ~ 400 Mya in our study but up to ~ 500 Mya based on mitochondrial markers in Wang et al. [63]. Despite this general trend, we also observed some shared families more complex to explain under a strict vertical evolutionary scenario, such as between Unionida and Venerida, as well as between Pinna nobilis (Ostreida, Pnib in Fig. 2B) and A. marissinica (Venerida; Amar in Fig. 2B).

Phylogenetic analyses of 150 random copies of the V and Meta superfamilies for each species (Fig. 3A-B) indicate that the great majority of elements are specific for a given bivalve order, as for Unionida, Mytilida, Arcida, and Venerida, while a few other elements are shared by different bivalve orders. On the contrary, the phylogenetic pattern of the Core superfamily is less clear as multiple groups of SINEs can be observed from the same bivalve order (Fig. 3C).

Fig. 3

Phylogeny of HCD superfamilies in bivalves. Phylogenetic trees of 150 random copies extracted from each genome for the superfamilies V (A), Meta (B), and Core (C). Colours of the tip labels represent bivalve order and reflect the colouring scheme of Fig. 2

SINEs accumulation in gene-related genomic regions and organization in complex tandem-like structures

To detect potential preferences in the genomic occurrence of SINEs with respect to coding regions, we carried out a case study using five species with available gene annotation, testing the hypothesis of a higher accumulation of older SINEs in gene-related compared to intergenic genomic regions. HCD SINE insertions were found within 0.9%, 10%, 8.8%, and 4.1% of the genes in C. gigas, M. californianus, P. maximus, and S. broughtonii, respectively, but reached 41% in R. philippinarum. Compared to the null expectation, exons and UTRs consistently exhibited significantly fewer insertions, while gene flanking and intergenic genomic regions were generally enriched with SINEs (Table 1). On the other hand, we observed a significant overrepresentation of insertions in the introns of R. philippinarum and S. broughtonii (Table 1), where we even found up to 61 and 145 insertions within a single gene, respectively (Table 1). These two species exhibited different accumulation patterns of SINEs within introns, with the former showing a low number of insertions in a high number of introns, while the latter showed a high number of insertions in a low number of introns (Supplementary Fig. 6).

Table 1 Genomic distribution of observed and simulated SINE insertions with respect to different genomic backgrounds

Moreover, for C. gigas, M. californianus, R. philippinarum, and S. broughtonii, gene flanking regions (2,500 bp flanking the gene) showed an enrichment of SINEs compared to intergenic ones (t-test; p-value < 0.01), whereas for P. maximus, we observed the opposite trend (Fig. 4A). Gene-related genomic regions (defined as exons, introns, UTRs and 2,500 bp gene-flanking regions) appear also characterized by older SINEs compared to intergenic ones, based on the Jukes-Cantor distance from consensus sequences across all species (Fig. 4B; t-test, all p-values < 0.001). Interestingly, we did not observe the same accumulation pattern when analysing the LINEs counterparts in C. gigas. Here intergenic genomic regions resulted significantly more affected by insertions compared to gene-flanking ones (Fig. 4C; t-test, p-value < 0.001; Supplementary Table 3) and particularly by older insertions (Fig. 4D; t-test, p-value < 0.001).

Fig. 4

Genomic occurrence of HCD SINEs. A Number of overlaps between SINE insertions with random gene-flanking regions (2,500 bp upstream and downstream the gene) and random intergenic genomic regions. B Jukes-Cantor (JC) distances of each SINE insertion from its consensus sequence as a proxy of the time of insertion. Gene-related = Insertions founded within genes (exons, introns, and UTRs) or in their 2,500 bp flanking regions. C and D are respectively specular to A and B but refer to LINE insertions in C. gigas. C Number of overlaps between LINE insertions with random gene flanking regions versus random intergenic genomic regions. D JC distance of LINE copies in gene-related versus intergenic genomic intervals. All comparisons are statistically significant (t-test, p-value < 0.01). Cgig = C. gigas, Mcal = M. californianus, Pmax = P. maximus, Rphi = R. philippinarum, Sbro = S. broughtonii

The same five genomes were also scanned for tandem-like HCD SINEs, considering only high-scoring insertions. All tandem arrays consist of two or three elements across all genomes, except for S. broughtonii, where we found 64 elements organized in tandem arrays of 4-15 units. Furthermore, while C. gigas hosts the smallest number of tandem-like SINE structures (three), in S. broughtonii, 3% of the high-scoring SINEs are organized in tandem arrays or palindromic structures, with 137 of them also incorporating one or multiple elements coming from a different family (Fig. 5A). The high number of tandem arrays in the blood clam S. broughtonii, even within intronic sequences (659 tandem arrays), could be an important contributor to the previously observed pattern of few introns impacted by a high number of insertions. We suggest that these SINE-rich introns could also drive the observed enrichment of SINE insertions in intronic sequences despite the low number of affected genes. Direct tandem arrays constitute most of tandemly repeated SINEs, while palindromes account for the 23% (403 structures) of which 126 overlap with gene annotations and particularly within intronic sequences (Fig. 5B-C).

Fig. 5

Tandem-like SINE structures in bivalve genomes: A Number of tandem-like SINE structures identified in each of the five analysed bivalve genomes. Tandem-like + Different SINEs means that together with tandem SINEs structures coming from the same family, we also detected elements coming from different families. Cgig = C. gigas, Mcal = M. californianus, Pmax = P. maximus, Rphi = R. philippinarum, Sbro = S. broughtonii. B and C examples of respectively direct and inverted SINE repeats present in intronic sequences of the S. broughtonii genome. D and E organization and structure of direct tandem-like SINE structures identified in the S. broughtonii genome. The tandem repeat in D is the same rapresented in B

Due to the high number of tandem-like SINEs in S. broughtonii, we looked deeper into the structures and arrangement of direct tandem arrays composed of more than 3 elements belonging to the same SINE family. We identified multimers of tandemly arranged SINEs with two main structures: TSD₁-(SINE-TSD₁)_n and TSD₁-(SINE)_n-TSD₁. The former is composed of highly similar SINE elements truncated at the same position of the 3’ tail and separated by the same TSDs (Fig. 5D), whereas the latter usually consists of a first SINE with a complete tRNA-related head but a truncated 3’ tail, multiple repetitions of a truncated SINE at both the 5’ tRNA head and 3’ tail, and a final SINE with the same tRNA head truncation but a complete tail (Fig. 5E).

Discussion

Transposable elements (TEs) are among the most significant sources of genetic variation across the eukaryotic tree of life. The advancement in the sequencing field are leading to a greater appreciation of TEs in the context of understanding genome evolution, gene regulation, and species diversification. However, while genomic resources are rapidly expanding for most eukaryotic clades, accurate TE identification and annotations are still lacking in non-model species [64], hindering our ability to comprehend their taxonomic distribution and effects on host biology. In this study, we took advantage of the increased number of bivalve genomes to comprehensively characterize tRNA-related HCD-containing SINEs across their diversity and to investigate their genomic distribution patterns. Our examination of 49 assemblies confirms that the five identified HCD-SINEs superfamilies have an ancient origin in bivalves, have been retained for a long evolutionary timescale [23], and can be derived from at least 10 different tRNA genes. Simultaneously, we observed important order-specific activity of the V and Meta superfamilies based on their phylogenetic clustering patterns. Based on analyses of LINE-derived tail regions, we found that at least four different LINE lineages (CR1; CR1-Zenon; L2; Nimb) can act as RT donor to four different SINE superfamilies for a total of seven SINE-LINE relationships of which three were previously unknown (specifically the partnerships between SINE V and LINE Nimb, SINE Core and LINE CR1, and SINE Meta and CR1-Zenon). Therefore, we greatly increase the number of putative SINE-LINE partnerships compared to the previous studies of Nishihara et al. [23] and Matetovici et al. [31]. The presence of multiple LINEs partner to SINE families with the same HCD, and the identified homologies between tail regions coming from families with different central domains suggest that recombination and module shuffling might have played a role in the emergence and subsequent amplification of new SINE families during bivalve diversification, coherently with a “mosaic” evolutionary scenario of SINEs [26]. Because of the strict relationship between SINEs and their LINE counterparts, we might expect almost overlapping landscapes of activity in the case of partnerships between the two elements. However, in multiple instances more complex evolutionary scenarios emerge, possibly due to different competitive dynamics for the LINE-derived enzymatic machineries [65, 66]. Indeed, specific SINE lineages could be particularly efficient in parasitizing their LINE counterparts, preventing them from expanding. The hijacked LINE, in turn, might increase its replication rate only when the SINE partially loses its parasitizing capacity and, consequently, its replication rate. Another limitation in inferring the co-evolutionary dynamics of SINEs and LINEs using repeat landscape profiles is the inability to account for different deletion rates among various transposons. Some TEs might be more susceptible to genome elimination compared to others, resulting in their underrepresentation in older divergence bins. These phenomena could contribute to the different patterns that emerge from our analyses. Indeed, despite a consistently significant and positive correlation between accumulation profiles of SINEs and their LINE counterparts, the strength of the correlation varied significantly between species, and the repeat landscapes showed substantial overlaps only in a few instances. Another possible explanation on the incongruence between SINE-LINE activity is that we might have missed the true LINE partner, annotating a sibling LINE lineage as the autonomous counterpart. While we tried to overcome this issue by attempting to reconstruct putative species-specific full-length partner LINE, we could not discard this possibility which should be further addressed in future studies.

Interestingly, when performing the cross-species family clustering we identified 22 HCD-SINE families characterized by the same highly similar modules in species separated by exceptionally long evolutionary time. For instance, the family Bpla_SINE-1_Meta was found in both Bathymodiolinae and Mytilinae which diverged more than 300 Mya in our study and ~ 400 million years ago in Lee et al., [61]. These results were strongly supported by the presence of multiple families shared between phylogenetically related but anciently diverging lineages when analysing the complete dataset of 45 bivalve genomes. Here we identified shared families even between groups of species that diverged between ~ 400 in our results and ~ 500 Mya in Wang et al. [63]. These represent exceptionally extreme cases of what was previously observed in grasses where SINE families were found to be far more conserved than LTR and TIR elements and retained for at least ~ 60 million years [34, 35]. The apparent long-term retention of TE families could also be explained by horizontal transposon transfer (HTT), which was already observed for SINEs in a few instances [28, 67, 68]. Indeed, we also found patterns of shared families difficult to explain under a strict vertical evolutionary scenario, such as between the Ostreida Pinna nobilis and the Venerida A. marissinica.

The ability to reconstruct consensus sequences shared between distantly related species could be favoured by the persistence of old insertions across the genome. In this context, it is interesting that, despite being underrepresented in exons and UTRs, SINE insertions tend to accumulate in gene-flanking regions (except for the Pectinida P. maximus) and, in the case of R. philippinarum and S. broughtonii, also within intronic sequences. It must be noted that we did not consider different GC contexts in our permutation tests, and different transposable elements have been shown to be enriched in genomic regions with different GC content [69, 70]. Even though we could not distinguish between preferential accumulation in GC-rich genomic regions, which are usually gene-dense [71], or in the genes per-se, the close association of HCD SINEs with gene bodies could increase the probability of their co-option from the host genome. Indeed, SINEs can frequently act as cis-regulatory elements, provide novel exons, or contribute to the mRNA processing, potentially enhancing the plasticity of tissue-specific transcripts, as recently observed in Drosophila [72]. A close association between SINEs and genes was also observed in plants [34, 35, 73, 74], fishes [75], mammals [76], and insects [67]. Open-chromatin genomic regions are known to be enriched in short and fragmented TEs [76, 77]. Furthermore, older SINE families were found to be more represented in euchromatic genomic regions compared to the younger ones both in grasses [34, 35] and in the coelacanth [75], consistently with our results. Indeed, we found that gene-related genomic regions (i.e., intragenic + 2,500 bp gene flanking regions) are enriched in older HCD SINE insertions, in terms of distance from their consensus sequence, compared to intergenic ones. If we assume no difference in insertion preference between old and young HCD SINEs, this pattern may suggest that gene-related regions could serve as safe ecological niches where short, non-coding transposons can survive. Coherently, we did not observe the same accumulation pattern when analysing LINE elements in the model bivalve C. gigas. One explanation is that short transposon insertions, like SINEs, could be favoured in the proximity of genes by a combination of (1) reduced competition with longer, more harmful TEs and (2) lower efficiency of TE-purging processes [34, 35, 78]. Indeed, deletions of transposable elements are mainly caused by ectopic DNA repair mechanisms, such as non-allelic homologous recombination and microhomology-mediated end joining [79, 80]. All these processes promote genome instability and may affect genomic flanking sequences, giving rise to complex and potentially harmful variants if a gene or a gene-interacting region is involved [81].

Methylation of SINE-derived direct repeats has been linked to the epigenetic regulation of downstream genes in Arabidopsis thaliana [82], and double-stranded hairpin structures in the mRNA derived from palindromic structures, resulting from alternating orientations of SINE insertions, might serve as substrates for DICER enzymes [74]. We found that HCD SINEs in bivalves can be organized in such tandem-like and palindromic structures also within genes, with an increased tendency in the Arcida S. broughtonii. In this species, approximately 3% of its HCD SINEs are organized in a similar manner. For comparison, in the potato genome, about 2% of the SINEs are included in tandem-like arrays [74]. Recently, Vassetzky et al. [83] also described common tandem repeats derived from different portions of SINE transposons in squamates [83]. Even though in the four species analysed here it appears to be a much less common phenomenon, it must be noted that we applied stringent cutoffs in terms of hit length (minimum length 150 bp) and RepeatMasker score (minimum 400) to avoid spurious matches. Therefore, we might have missed tandem structures composed of highly fragmented and/or diverging SINEs, and more systematic and targeted screening is needed to more precisely quantify the impact of SINEs in the formation of tandem repeats in bivalves. Finally, the high number of SINE direct repeats identified in S. broughtonii raises interesting hypotheses about their potential origin and genome evolutionary dynamics of this species. Indeed, one possible outcome of unequal homologous recombination between target site duplications (TSDs) is the formation and the expansion of SINE tandem arrays [84]. The high number of such structures in S. broughtonii could therefore imply higher recombination rates in this species compared to other analysed bivalves.

Conclusions

Here we present for the first time a wide characterization of tRNA-related HCD SINEs in bivalves looking at their distribution, LINE partnerships and genomic occurrence. Thanks to a novel, manually curated SINE library we found that bivalves HCD SINEs could derived from at least 10 different tRNAs and depend on at least four different LINE lineages. Different LINEs can promote the amplification of the same HCD superfamily whereases homologues tails are shared between different superfamilies, suggesting a “mosaic” evolutionary scenario of SINE modules. Additionally, some SINE families are apparently shared between distantly related species underlying the possible long-term retention of highly similar HCD SINE linages characterized by the same tRNA-related head, central domain and LINE-derived tail. Genomic occurrence analyses across five different bivalve species highlighted their potential different effects in genome evolution. Indeed, different species show overrepresentation of SINE insertions across different genomic compartments as well as different tendencies to form tandem-like and palindromic structures which could be present in intronic sequences. Despite these differences, we found a consistent trend of accumulation of old SINEs in close proximity of genes, as previously observed in plants and other metazoan. This result suggest that evolutionary dynamics of SINEs might partially follow a common evolutionary route across eukaryotes in which euchromatic genomic regions serve as safe niches for their survival. Overall, this study represents a step forward in a broader understanding of TE evolutionary dynamics in a highly overlooked but important taxonomic group like bivalves, and open interesting questions about the possible role of SINEs in bivalve genome organization, biology and evolution.