Impact of a horizontally transferred Helitron family on genome evolution in Xenopus laevis

Abstract

Background

Within eukaryotes, most horizontal transfer of genetic material involves mobile DNA sequences and such events are called horizontal transposable element transfer (HTT). Although thousands of HTT examples have been reported, the transfer mechanisms and their impacts on host genomes remain elusive.

Results

In this work, we carefully annotated three Helitron families within several Xenopus frog genomes. One of the Helitron family, Heli1Xen1, is recurrently involved in capturing and shuffling Xenopus laevis genes required in early embryonic development. Remarkably, we found that Heli1Xen1 is seemingly expressed in X. laevis and has produced multiple genomic polymorphisms within the X. laevis population. To identify the origin of Heli1Xen1, we searched its consensus sequence against available genome assemblies. We found highly similar copies in the genomes of another 13 vertebrate species from divergent vertebrate lineages, including reptiles, ray-finned fishes and amphibians. Further phylogenetic analysis provides evidence showing that Heli1Xen1 invaded these lineages via HTT quite recently, around 0.58—10.74 million years ago.

Conclusions

The frequently Heli1Xen1-involved HTT events among reptiles, fishes and amphibians could provide insights into possible vectors for transfer, such as shared viruses across lineages. Furthermore, we propose that the Heli1Xen1 sequence could be an ideal candidate for studying the mechanism and genomic impact of Helitron transposition.

Background

The transmission of genetic material from parental organisms to their descendants usually occurs in a so-called vertical mode. Still, horizontal transfers (HT) between organisms are known to occur, providing evidence for DNA movement across mating barriers. Such HT happens frequently in prokaryotes to drive gene acquisition and adaptive evolution. In eukaryotes, HT events are increasingly reported, especially those involving repeated genomic components, essentially transposable elements (TE). Such events are called Horizontal TE transfer (HTT) and are significant actors in eukaryote genome evolution, including metazoans [1].

TEs are DNA sequences that can “jump” from one genomic locus to another; in doing so, they usually duplicate and can finally make up a considerable proportion of genetic materials [2]. TEs can be classified into Class-I (RNA TE) and Class-II (DNA TE) based on their sequence homology and their transposition mechanisms [2]. In many cases, the new arrival of TEs in a genome induces additional mutations and eventually can be co-opted by the host for additional functions [3]. However, in most circumstances, TE-induced mutations are deleterious to their host organisms, and multiple genetic and epigenetic mechanisms involving piwi-interacting small RNAs (piRNAs) and Kruppel-associated box-containing zinc finger proteins (KRAB-ZFPs) pathways have evolved to suppress TE activities at both the transcriptional and post-transcriptional level [4]. Once their activity is suppressed, TEs gradually lose the ability to transpose and eventually become transposon relicts before being extinct due to genetic drift. An explanation for the long-term persistence of TEs across genomes is that HTT can serve as an escape mechanism. This allows TEs to evade “suppressor” hosts that can repress their activity, enabling them to colonize “naive” genomes where their replication is supported [4].

Since the first report of HTT, numerous other cases have been reported, and the surge of genomic sequencing revealed that these events are more frequent than initially believed [5, 6]. A recent large-scale survey has identified nearly a thousand such HTT events, mostly involving ray-finned fishes [6]. Zhang and collaborators showed that Class-II TEs are significantly more engaged in HTT, evolved neutrally within their genome of origin, and under purifying selection after an HTT event.

Among different types of Class-II TEs, Helitrons are especially noticeable because they are the only known eukaryotic TE superfamily replicating via rolling-circle replication [7]. Helitrons can be classified into two groups, Helitron1 and Helitron2, which differ in their terminal features and transposase homology but are widespread throughout the eukaryotic world [8, 9]. Numerous studies showed they can capture host genes and reshuffle host genomes, significantly influencing host evolution [10,11,12,13]. For example, Helitron1s invaded the vesper bat lineage about 30–36 mya, made up ~ 3% of the M. lucifugus genome, and amplified the nucleotide binding protein-like (NUBPL) gene to about 1000 copies [14].

Still, many factors and mechanisms that govern HTT remain to be elucidated, such as viral transduction, parasitism and symbiosis, and endocytosis of extracellular vesicles [15,16,17]. In particular, the encapsulation of LINE1 retrotransposons into extracellular vesicles followed by HTT between cultivated human cells and the encapsulation of various TEs into large double-stranded DNA virus particles infecting arthropod hosts have received experimental support [17]. How much transfer occurs between two species or how much variability there is in the efficiency of HTT between major lineages requires the analysis in selected organisms, including model species for which a large body of genomic data is available. Since it has been shown that amphibian genomes are especially prone to HTT and since Helitrons are understudied class-II TEs [6], we investigated in more detail the landscape of Helitrons in Xenopus, the most studied amphibian model organism.

Xenopus, belonging to the Pipidae family, is a genus of African clawed frogs, including two subgenus: Silurana and Xenopus. Over the past half-century, these frogs, especially X. tropicalis (Silurana subgenus) and X. laevis (Xenopus subgenus), have served as important model organisms in cell and developmental biology. With the advent of the genomic era, they have gained increasing prominence in studies of genome evolution due to their frequent whole-genome duplication (WGD) events, which rapidly produced a range of polyploids, from diploids to dodecaploids. These polyploids are allopolyploids where genomic dominance would happen after WGD, resulting in faster genomic DNA loss on one subgenome (S subgenome) than another one (L subgenome) [18]. Correspondingly, abundant omics datasets are now available from public databases, including the genome assemblies of X. tropicalis (2n = 2x = 20), X. laevis (2n = 4x = 36) and X. borealis (2n = 4x = 36) [18]. X. tropicalis is the only known diploid species of Xenopus frogs, and X. laevis and X. borealis are allotetraploids whose ancestor originated from the hybridization of two diploid Xenopus species. X. tropicalis diverged from X. laevis and X. borealis ~ 48 mya, and X. laevis diverged from X. borealis ~ 17 mya [18].

In this work, we carefully annotated the Helitron families for Xenopus genomes and compared their distribution pattern across three species. We further investigated the impacts of Helitron amplification on their hosts, including possible Helitron-mediated gene capture events, the Helitron expression and insertion polymorphisms within X. laevis population. Our findings revealed that the Heli1Xen1 family exhibited significant activity in all three species. Then, we explored the prevalence of Heli1Xen1 and found its presence in 14 vertebrate species from divergent vertebrate lineages, including reptiles, ray-finned fishes and amphibians. Further phylogenetic analysis confirmed that this Helitron family had recently been horizontally transferred among these lineages. Altogether, our results demonstrate that the Heli1Xen1 is very young and might still be active in X. laevis, underscoring the significant role of Helitrons in genome evolution. Moreover, the frequent Heli1Xen1-mediated HTT events among distinct lineages that share hosts with certain viruses support the models of viruses acting as vectors for HTT [15].

Methods

Annotation of Helitron sequences for Xenopus species

We downloaded the reference genome assemblies of X. tropicalis (v10.0, GCF_000004195.4), X. laevis (v10.1, GCF_017654675.1) and X. borealis (GCA_024363595.1) from NCBI. To annotate the Helitrons for these species, we applied HELIANO (v1.2.1) on each genome with the parameter ‘-s32 0 -is1 0 -is2 0 --nearest -flan_sim 0.7 -sim_tir 0.9’.

Phylogenetic analysis for Helitron sequences from Xenopus

We predicted the RepHel transposase sequences of autonomous Xenopus Helitron insertions with the program getorf (EMBOSS:6.6.0.0) with the parameter ‘-minsize 400’. To analyse their phylogenetic relationship, multiple alignment was carried out for RepHel sequences, including two from Repbase (HeliNoto and Helitron- 9_OS) and three outgroup sequences generated manually. The HeliNoto (from the fish Chionodraco hamatus) and Helitron- 9_OS (from the plant Oryza sativa) were chosen as positive controls because they are known autonomous Helitron2 and Helitron1 sequences from Repbase. The outgroup sequences were generated by concatenation of the Rep sequences from plasmid/phage and the Pif1 Helicase sequences from vertebrates based on the phylogeny of Rep and Hel domains [19]. Outgroup1 is the concatenation of replication protein from the plasmid of Brevibacillus borstelensis (BAA07788.1) and the pif1 protein from X. tropicalis (XP_004914012.2); Outgroup2 is the concatenation of replication protein from plasmid pVT736 - 1 and pif1 protein from X. laevis (NP_001089530.1). Outgroup3 is the concatenation of Replication-associated protein from phage (AAA88392.1) and human pif1 protein (KAI4058258.1). After alignment, we applied FastTree (v2.1.11) with default parameters to build the phylogenetic tree [20].

Identification of tandem arrays that originated from Helitron sequences

To get the tandem array annotation for the X. laevis genome, we ran the program trfBig from the ucsctools package with default parameters. We then extracted the monomer sequence of each tandem array sequence and compared the monomer sequence of tandem arrays with Helitron sequences using blastn. The tandem arrays with monomer sequences sharing > = 80% identity with Helitron sequences and covering > = 85% of Helitron sequences were considered Helitron-derived.

Detection of gene capture events by Heli1Xen1

We collected all coding sequences (CDS) from the X. laevis NCBI gene annotation project. The genes with a ‘LOC’ marker were excluded from the analysis. Then, we compared collected CDS sequences and Helitron sequences using blastn with default parameters. Hits with more than 80% identity and at least 100 bp matched length were kept.

Characterization of Heli1Xen1 expression in X. laevis

To check if the Heli1Xen1 loci in the X. laevis genome could be combined with any regulatory elements, we extracted the Chip-seq profiles of foxn4 and myb in bigwig format from Xenbase [21]. There are 28 Chip-seq transcription factor profiles available from Xenbase representing 13 different transcription factors, among which six of them (myb, foxn4, beta Catenin, ep300, rfx2, foxj1) can be found enriched in the 5’ of the Heli1Xen1 locus on Chr3L. We selected myb and foxn4 as examples because their peaks at this region are the most obvious. To check if the Heli1Xen1 loci are transcribed in X. laevis, we collected the RNA-seq data from a previous study (PRJDB2519) [22]. The RNA-seq reads were extracted from the X. laevis blastema cells in proliferating and non-proliferating stages. Then, we mapped the RNA-seq reads with STAR (v2.7.10b) onto the X. laevis reference genome (v10.1) [23] with parameters ‘–runMode alignReads –winAnchorMultimapNmax 100 –outFilterMultimapNmax 100 –outSAMattributes All –outSAMstrandField intronMotif –outSAMattrIHstart 0’. Finally, we extracted the uniquely mapped reads with samtools (-q 255) and visualized the results with JBrowse [24].

Characterization of Heli1Xen1 insertion polymorphism in X. laevis

Genomic reads from three outbred clutches of embryos from the X. laevis Nasco strain (Xla.pigmentedNasco in Xenbase) were collected from a previous study (PRJNA952800) [25]. These genomic reads were mapped against the X. laevis reference genome (v10.1) with bwa-mem2 (v2.2.1) with default parameters, and the output of sorted bam files obtained with samtools (v 1.16.1) were used as inputs of the program MEGAnE (v1.0.1) to detect the potential Helitron caused absence and insertion [26,27,28]. Candidates in MEGAnE output with the ‘PASS’ flag were selected for further analysis.

Discovery of Heli1Xen1 insertions among various vertebrate species

To explore Heli1Xen1 distribution across eukaryotic genomes, we used the Heli1Xen1 RepHel transposase nucleotide sequence as a query to search the NCBI whole-genome shotgun database with blastn (v 2.13.0 +). We used only the RepHel sequence instead of the whole-length sequence as a query because Helitrons are known to be variable in length, which obscures the interpretation of sequence similarity results. Genome assemblies containing highly similar copies to Heli1Xen1 RepHel (> = 80% identity and > = 80% coverage) were then collected and downloaded from NCBI (Supplementary Table S1). To obtain all Heli1Xen1 homologous sequences for each genome, we ran HELIANO (v1.2.1) de-novo using the parameter ‘-p 1e- 5 -s 32 –nearest -flank_sim 0.6’. We then used the Heli1Xen1 consensus as a query sequence to run a blastn search against all HELIANO predictions. Each species'predicted families who shared > = 80% identity and > = 80% coverage with the Heli1Xen1 query were regarded as members of the Heli1Xen1 family.

Making consensus sequences

We performed multiple alignments of all predicted Heli1Xen1 autonomous copies for each species using mafft (v7.526) with the parameter ‘–auto’. The consensus sequences were obtained using cons (EMBOSS:6.6.0.0) with the default parameter. Only one autonomous Heli1Xen1 insertion was found in genomes of the reptiles Hemicordylus capensis, Lacerta agilis and Lacerta bilineata; and of the frogs Hymenochirus boettgeri and Xenopus borealis, and they were used as representative sequences. All sequences are available in Supplementary Data.

Estimation of Heli1Xen1 conservation score

We obtained consensus sequences from each Heli1Xen1 group (Heli1Xen1-a or Heli1Xen1-b) and each corresponding species and made multiple alignments via mafft (v7.526) with default parameters. We then applied phyloFit with default parameters from the PHAST package (v1.5) to train phylogenetic models. Afterwards, we ran phastCons (–target-coverage 0.25 –expected-length 30) to estimate the phastCons score representing each base's conservation score [29]. The consensus sequences of X. laevis and P. senegalus were selected as reference sequences for Heli1Xen1-a and Heli1Xen1-b, respectively.

Estimation of Kimura substitution for protein-coding genes

Among those vertebrates whose genomes contain full-copy Heli1Xen1, the genome annotation of eight species is available from NCBI, and they are the amphibians Bufo bufo, Bufo gargarizans and Geotrypetes seraphini; the reptiles Hemicordylus capensis, Lacerta agilis and Zootoca vivipara; and the ray-finned fishes Erpetoichthys calabaricus and Polypterus senegalus. To extract the orthologous gene groups among those vertebrates, we used the annotation for X. laevis as the reference. We used its unique gene set as queries to extract the corresponding orthologous genes in other species. Since X. laevis is a tetraploid species, there is a majority of duplicated genes (also called homoeologs or ohnologs), labelled with either the suffix “L” or “S” depending on their location on a long or short chromosomal copy. The unique gene set of X. laevis was made by collecting all singleton genes (either L or S form) together with only the L copy of duplicated genes. We compared each species' gene catalogue to the reference with blastn (v 2.13.0 +) with default parameters and defined the reciprocal best hits as orthologous pairs. Afterwards, we applied MACSE (v2.07) with default parameters to each orthologous gene group for codon alignments [30]. Finally, we estimated each group's Kimura substitution (Ks) with the function kaks in the seqinr package [31].

Phylogenetic analysis

We fetched species trees from the Time Tree database [32]. To obtain the phylogenetic relationship of Heli1Xen1, we carried out the following steps: 1) multiple alignments of consensus sequences using mafft (v7.526) with default parameters; 2) trimming of the alignment using trimal (v1.4) with the parameter ‘-strictplus’; 3) construction of a maximum-likelihood tree using FastTree (v2.1.11) with default parameters [20, 33].

Results

Helitron distribution among Xenopus species

Although Xenopus frogs are important model species for genomic studies, their Helitron repertoire has not been deeply investigated. Indeed, the published annotations of high-quality reference genomes for the diploid and allotetraploid Xenopus species report only three non-autonomous Helitrons in Xenopus tropicalis [18, 34]. In this work, we annotated Helitrons with our recently developed tool HELIANO on three Xenopus frog genomes [8], one (X. tropicalis) from the Silurana subgenus and two (X. laevis and X. borealis) from the Xenopus subgenus. After curation, we annotated 84 insertions (three autonomous) in the X. tropicalis genome, 2,310 insertions (15 autonomous) in the X. laevis genome and 1,861 insertions (two autonomous) in the X. borealis genome (Supplementary Table S2). In X. tropicalis, Helitron insertions covered 0.01% (159 kbp) of the genome, while the coverage reached 0.13% (1,789 kbp) in X. laevis and 0.05% (1,254 kbp) in X. borealis. These insertions were dispersed over all chromosomes (Supplementary Table S2).

We then asked how these Helitron sequences evolved across these Xenopus species. We found that these Helitron sequences belonged to three families according to sequence homology: two Helitron1 families, Heli1Xen1 and Heli1Xen2, and one Helitron2 family, Heli2Xen (Fig. 1A). For Helitron1, we detected the Heli1Xen1 family in the X. borealis and X. laevis genomes and the Heli1Xen2 in the X. laevis genome only (Fig. 1A). For Helitron2, we found the Heli2Xen family in X. tropicalis and X. laevis genomes (Fig. 1). In association with the full-length Helitron sequences found in these Xenopus genomes, we identified nearly three orders of magnitude more non-autonomous Helitron copies (Fig. 1B). We examined the length distribution of these non-autonomous insertions and found abundant miniature copies of a few hundred base pairs for Heli1Xen1 (Helitron1) and Heli2Xen (Helitron2) families in X. laevis (Fig. 1B). We found that the landscape of Heli1Xen1 non-autonomous insertions was particularly remarkable with more than one thousand copies and a minimal size of 152 bp in X. laevis genome. The pattern observed for Heli1Xen2 was markedly different, with a peak at 2,067 bp (Fig. 1B). For the Heli2Xen family, the distribution differed between X. laevis and X. tropicalis genomes: the latter was characterized by a relatively flat distribution and only a few non-autonomous insertions, while there was a clear peak of 480 bp insertions, found with 670 copies in X. laevis genome.

Fig. 1

Helitron distribution in the genome of Xenopus frogs. A Phylogenetic tree of Xenopus Helitron transposases. The HeliNoto (from the fish Chionodraco hamatus) and Helitron- 9_OS (from the plant Oryza sativa) were chosen as positive controls because they are known autonomous Helitron2 and Helitron1 sequences from Repbase. Outgroup1 - 3 are three outgroup sequences. B Distribution of non-autonomous Helitrons in three Xenopus frogs. Frequency along the y-axis of Helitron insertions according to their length (x-axis, in bp). XT: X. tropicalis; XB: X. borealis; XL: X. laevis. The y-axis is square root adjusted. The star in red represents the WGD event that happened ~ 34 mya

Previous studies showed that some Helitron families can generate tandem arrays in genomes due to their rolling circle replication mechanism [35, 36]. Therefore, we asked if some Xenopus Helitron sequences were organized in such tandem arrays and found 43 tandem repeats in the X. laevis genome overlapping Helitron sequences (Supplementary Table S3), indicating their Helitron origin. The monomer size of these tandem arrays varied between 121 and 1,996 nt, and the number of repeats varied between two and seven. Moreover, we found that 39 (~ 91%) of these tandem arrays were homologous to the family Heli2Xen, and the remaining four were homologous to the family Heli1Xen1. A closer inspection of these tandem arrays revealed that all Heli2Xen cases are made of minisatellite repeats of 166 bp. These minisatellite arrays are also found on both ends of the unique Heli2Xen autonomous copy found on chromosome 3L (Supplementary Table S2, Supplementary Figure S1). When we analysed all non-autonomous Heli2Xen insertions, we found that these minisatellite arrays occupy the major part of the insertions and display sequence variation.

In summary, we observed a mosaic distribution of Helitrons across African clawed frogs, suggesting rapid turnover through extinction and acquisition, particularly between the Silurana and Xenopus subgenera. Additionally, we found that Helitron amplification patterns varied among frog species, even within the same Helitron family, exhibiting differences in copy number and non-autonomous element length. Furthermore, we observed a higher prevalence of tandem arrays originating from Helitron2 than Helitron1.

Non-autonomous Heli1Xen1 captured coding sequences in the X. laevis genome

Helitron can capture and disseminate nearby genomic sequences, resulting in gene duplication, genome rearrangement, and epigenetic conflicts [10,11,12, 37, 38]. Therefore, we wanted to characterize such gene capture events by Helitrons in Xenopus genomes. We focused on the X. laevis genome because it contains all three Helitron families with numerous insertions and has a good gene annotation. We collected all X. laevis protein-coding sequences (CDS) and compared them to Helitron sequences. We found that only Helitrons of the Heli1Xen1 family contained fragments of CDS, with 19 fragments distributed across 64 Heli1Xen1 insertions, which shared the same pattern in LTS and RTS region with that of autonomous Heli1Xen1 insertions (Fig. 2A). Moreover, 13 gene CDS fragments were captured and distributed multiple times (Fig. 2B-E, Supplementary Table S4). For example, part of the CDS of armadillo repeat containing 5 gene (armc5.L) has been copied 17 times, and the CDS of H2.0 like homeobox gene (hlx.L) has been copied six times.

Fig. 2

Landscape of captured gene fragments within Helitron insertions. A Sequence logos of terminal sequences of LTS and RTS for nonautonomous Heli1Xen1 involving in gene capture events. Red arrows on top of the graph indicate the pattern of Helitron insertion sites. B Nucleotide identity between the donor gene and captured fragments. C Length of captured CDS fragments by Helitron sequences. D Relative start position of captured fragments within the donor gene, zero corresponds to the 5’ end of the CDS and 0.8 to the 3’ end. E Number of Helitron insertions that contain the gene fragments

We then characterized the features of these captured gene sequences by Heli1Xen1 Helitrons. We found that almost all of these captured CDS fragments belonged to the 5’ end of their genes (17 out of 19, Fig. 2D), which is consistent with previous studies [39, 40]. The two exceptions are armc5.L and transmembrane protein 111 (tmem111.S) genes. The captured armc5.L gene fragments belonged to its 3’ end, while the captured tmem111.S fragments overlapped with the complete transcript of the tmem111.S gene. We found that the Heli1Xen1 family captured almost the full transcript of the tmem111.S gene and amplified it to two copies, distributed on chromosomes 4L and 7L (Supplementary Figure S2). Altogether, our results show that the Heli1Xen1 family has captured and amplified 19 gene fragments in the X. laevis genome.

Heli1Xen1 is probably still active in the X. laevis genome

Since we found that Heli1Xen1 Helitrons captured gene sequences, we wondered if this family is still active in the X. laevis genome. We detected three autonomous insertions, with one of the insertions on Chr3L containing a complete uninterrupted ORF structure, annotated by the NCBI genome annotation pipeline as ‘uncharacterized protein LOC108706566’ (Fig. 3 and Supplementary Figure S3). This ORF encodes a 1401 aa protein with Rep and Helicase domains necessary for transposase function (Supplementary Figure S3). We further looked at the transcription landscape of this locus using publicly available datasets from Xenbase [21]. We found pieces of evidence from Chip-seq profiling of two transcription factors (foxn4 and myb) at embryonic neurula stage NF18, with significant peaks for foxn4 and myb at the 5’ end of this locus corresponding to the presumed promoter region, indicating that this locus could be a target of the two transcription factors which can initiate transcription (Fig. 3). We also mined publicly available transcriptome data [22] and remapped them onto the genome. Some RNA-seq reads corresponding to blastema cells were uniquely mapped on this locus, although the coverage was low (Fig. 3).

Fig. 3

Genomic landscape of the X. laevis Heli1Xen1 locus. This locus is annotated as gene LOC108706566 (the blue part indicates the 5’ UTR region, and the yellow part indicates the CDS region) by the NCBI annotation pipeline and as a Helitron insertion by HELIANO. The Chip-seq profile of animal cap NF18 shows significant foxn4 and myb Chip-seq peaks (TPM) at the 5’ end of this locus. The RNA-seq reads collected from blastema cells in either non-proliferating or proliferating status are uniquely mapped on this locus [22]

Finally, we reasoned that TE polymorphisms between X. laevis frogs of distinct populations could also indicate Heli1Xen1 TE activity. We collected resequencing data of two X. laevis samples of the same outbred Nasco strain generated by a previous study and searched structural polymorphisms of Heli1Xen1 insertions using MEGAnE. As a result, we detected 90 insertions and 56 absences of Heli1Xen1 in the outbred genomes in comparison to the inbred J strain reference genome, indicating that Heli1Xen1 contributed to genomic polymorphisms in X. laevis population (Fig. 4A, Supplementary Figure S4). Insertion lengths varied from 66 to 2,136 bp, and deletion lengths ranged from 131 to 3,017 bp. Remarkably, we observed a prominent peak at approximately 200 bp in the size distribution of these polymorphisms, accounting for about 57.78% of insertion polymorphisms and 37.50% of absence polymorphisms (Fig. 4B).

Fig. 4

Heli1Xen1 polymorphism in X. laevis genomes. A Heli1Xen1 polymorphism distribution across chromosomes of X. laevis. The blue bar represents the density of Heli1Xen1 insertions every 100 kbp window in the X. laevis reference genome. Dots in purple represent the insertion of Heli1Xen1 in the outbred genomes. Green triangles represent the absence of Heli1Xen1 in the outbred genomes. B The histogram shows the frequency of detected Heli1Xen1 polymorphisms according to their length expressed in bp

We conclude that the Heli1Xen1 family might be transcriptionally active and have generated genomic polymorphisms within the X. laevis population.

Heli1Xen1 sequences in other vertebrates

Because Heli1Xen1 showed signs of activity in X. laevis, we investigated the possibility of its distribution in other species. We looked for homologous sequences in the whole-genome shotgun database and found Heli1Xen1 homologs with high levels of similarity in 13 vertebrate genomes from two reptile families: Lacertidae and Cordylidae; three amphibian families: Dermophiidae (caecilian), Pipidae and Bufonidae (frogs), and one ray-finned fish family: Polypteridae. We annotated all Heli1Xen1 insertions in these genomes using HELIANO and detected numerous non-autonomous Heli1Xen1 insertions, i.e., from hundreds to more than twenty thousand (Fig. 5). As expected, we found fewer autonomous Heli1Xen1 insertions (one to 61), fitting the ‘master copy’ model [41] (Fig. 5).

Fig. 5

Heli1Xen1 sequences in vertebrates. The heatmap (left) shows the pairwise identity of Heli1Xen1 RepHel consensus sequences from the different vertebrates where it was found. Each cell indicates the percentage of nucleotide identity. The bar plot (right) shows the number of autonomous and non-autonomous Heli1Xen1 insertions in each vertebrate genome. The phylograms on the left and the top of the heatmap are based on Time Tree database

To examine the homology of Heli1Xen1 across different species, we built consensus sequences for autonomous insertions from each species and performed pairwise alignments. Because Helitrons can capture non-Helitron sequences, which would dramatically increase the sequence divergence, we extracted the RepHel region from each consensus sequence and estimated their pairwise identity. Overall, we observed that the pairwise identity for each pairwise alignment group is high, from 81 to 100%, which convinced us of their homology (Fig. 5). In addition, we could classify these consensus sequences into at least two groups, Heli1Xen1-a and Heli1Xen1-b, whose pairwise identity within the group is much higher (~ 95%) than outside group (Fig. 5). The Heli1Xen1-a group contained sequences from Xenopus, Lacerta and Zootoca genomes, and the Heli1Xen1-b group contained sequences from Bufo, Polypterus and Erpetoichthys genomes.

Further phylogenetic analysis of these RepHel sequences confirmed the two groups, suggesting different phylogenetic relationships (Fig. 6). Then, we analyzed the conservation across the consensus sequence for each group. The result showed that the RepHel domain of the Heli1Xen1 sequence is much more conserved than terminal regions, and the LTS region is much more variable than the RTS (Fig. 6B, C). Furthermore, we found that RepHel of Heli1Xen1-b conservation is high all over the sequence, while there are conservation drops in HeliXen1-a, possibly indicating significant mutations in the ORF for this group (Fig. 6B, C).

Fig. 6

Classification of homologous HeliXen1 sequences across vertebrates. A Phylogenetic relationships of HeliXen1 RepHel sequences from different vertebrates. Heli1Xen1-a and Heli1Xen1-b are highlighted in different colours. Heli1Xen2 and Helitron- 9_OS are used as the outgroup. B-C Plots show the PhastCons score across Heli1Xen1-a (B) and Heli1Xen1-b (C) consensus sequences from X. laevis and Polypterus senegalus (a ray-finned fish), respectively. The position of RepHel ORF is highlighted in red

Evidence for Horizontal transfer of Heli1Xen1

Such a high level of nucleotide similarity between sequences of such divergent species is unexpected and questions the scenario of transmission via vertical inheritance of these Heli1Xen1 Helitrons. An alternative explanation would be that HTT events enabled the spread of Heli1Xen1 across these species or their close ancestors. HTT events can be inferred when the nucleotide divergence between TE copies from two host lineages is much lower than the expectation estimated from the divergence of orthologous genes.

To estimate the nucleotide divergence of homologous Heli1Xen1 sequences, we used the Heli1Xen1 consensus of X. laevis as the reference. We estimated the Kimura substitutions between the reference and those from other species. We collected a unique gene set from X. laevis to assess the divergence between orthologous genes, which were then used as the reference gene set. We identified each unique gene's orthologous gene and estimated the Kimura substitution (Ks) for species whose gene annotation is available. Finally, we plotted the distribution of Ks values for orthologous gene pairs. The result showed that the median Ks value for genes varied from ~ 1.27 to 2.17 among species (Fig. 7). In contrast, we observed that the Ks value for Heli1Xen1 was significantly lower than that for genes, from 0.05 to 0.17 (Fig. 7). This result indicated that Heli1Xen1 was likely transmitted horizontally rather than vertically among these vertebrates.

Fig. 7

Histogram of the Kimura substitution distribution between orthologous genes of X. laevis and other vertebrates. The red dotted line indicates the Kimura substitutions between Heli1Xen1 from X. laevis and each vertebrate genome

Species-specific amplification of Heli1Xen1

To understand the amplification pattern of this Helitron family across the species phylogeny, we analyzed the length distribution and terminal features of Heli1Xen1 sequences for each species. We detected some huge non-autonomous insertions in a few species, but they occurred rarely and were mainly caused by nested TEs, which inflated the Helitron length. The non-autonomous insertions were always significantly shorter than autonomous ones by an order of magnitude (Fig. 8). We highlighted the peaks of the non-autonomous length distribution, which represented the most frequent Helitron sizes. We found variations ranging from 152 bp (X. laevis) to 1,747 bp (B. bufo) (Fig. 8). The size distribution also varied significantly, most notably for the G. seraphini genome and to a lesser extent for Bufo genomes, with a wide distribution of a few tens of non-autonomous elements of up to 20 kbp, in correlation with the total number of non-autonomous Helitrons.

Fig. 8

Length distribution of Heli1Xen1 insertions in different species genomes. Non-autonomous insertions are shown in red, and autonomous insertions in blue. Species are ordered based on their phylogeny relationship. The X-axis represents the insertions' length (bp), and the Y-axis represents the number of insertions log10-adjusted. The dotted line represents the peaks of size distribution

Further analysis of their terminal sequences revealed that Heli1Xen1 for each species started with ‘TC’ and terminated with ‘CTAG’, which fulfilled terminal features of canonical Helitrons. Besides, we observed conserved patterns for both LTS and RTS: e.g., the TA-rich region (TATATAAA) for LTS and the hairpin pattern (‘CCCGTTA…TAACGGG’) in RTS, which are the other characteristic features of canonical Helitrons (Supplementary Figure S5). The conservation of these terminal sequence features argues for a well-supported delineation of Heli1Xen1 insertions that follow the definitions of bona fide Helitrons. However, we found that the LTS and RTS were not the same among these species, especially for the LTS, which was more variable than the RTS (Supplementary Figure S5). Altogether, these results suggest that the amplification pattern of Heli1Xen1 in these species is species-specific.

The recent invasion of Heli1Xen1 in vertebrate genomes

Because TEs evolve mostly neutrally since their insertions, their invasion time can be approximated by measuring the number of substitutions among copies and the neutral substitution rates characteristic of the species lineage [42]. To estimate the neutral substitution rates, we collected the divergence time between X. laevis and the other five classes from the Time Tree database [32] and their synonymous substitutions that were estimated before. For Pipidae frogs, we adopted the substitution rate of 3.20 × 10^–9 substitutions per year based on the estimation from a previous study [18]. Finally, based on the number of substitutions between Heli1Xen1 copies, we inferred that the full-length Heli1Xen1-a invaded Xenopus about four mya and invaded Zootoca less than one mya (Table 1). Similarly, the full-length Heli1Xen1-b invaded Bufo about five mya and invaded fish about seven mya (Table 1). This result indicated that the Heli1Xen1 family is relatively young and that several independent transfer events occurred even between closely related species. Based on the more recent estimation of the invasion time (Table 1) than the Time Tree inference of divergence time, independent Heli1Xen1 transfer happened in X. laevis and X. borealis, and in B. bufo and B. gargarizans.

Table 1 Table for Heli1Xen1 Ks value and invasion time in each species

Discussion

Helitrons can be difficult to identify because they lack the obvious features of some class I and class II TE and because only the RepHel transposase is a well conserved feature in autonomous elements. The detailed analysis of Helitrons in Xenopus genomes leads to empirical evidence on the rapid extinction and acquisition of Helitrons, the high frequency of coding exon capture and amplification by Helitrons, the recent and ongoing transposition activity of the Heli1Xen1 elements in the X. laevis genome and the recurrent HTT of Heli1Xen1 to some amphibians, squamates and actinopterygians species.

The first genome-scale analyses of the X. tropicalis genome sequence led to the identification of the non-autonomous Helitrons: Helitron-N1_XT, Helitron-N1 A_XT and Helitron-N2_XT sequences [43,44,45,46]. We did not identify Helitron-N1_XT and Helitron-N1A_XT sequences using HELIANO because they are unrelated to a complete autonomous insertion, a limitation underscored in our previous work, and the Helitron-N2_XT was further proved belonging to the family Heli2Xen [8]. The short description of Helitron-N2_XT and Helitron-N1_DR (a zebrafish Helitron) notes that these two transposons are all composed of minisatellite repeats and are likely derived from an ancestral HTT [46, 47]. Given an estimated Kimura substitution of 0.24 for the two sequences and the substitution rate 2.52e- 9, we estimated the approximate divergence time between the two sequences to be ~ 47.62 mya, which is indeed much more recent than the speciation time between X. laevis and Danio rerio (430 mya from Time Tree database) and therefore corresponds to another HTT. The observation that some families of Helitrons can be composed in large part of minisatellite repeats has been made on elements found in various species, such as DINE and HINE in Drosophila, in bivalves and bats. Our findings that Heli2Xen autonomous element is bounded by minisatellite arrays further documents these previous observations and provides a nice example linking an autonomous Helitron with the birth and dissemination of minisatellite repeats [48].

Helitron-mediated capture of host genes is well known and corresponds to a transduplication event in which the captured gene sequence is located between the TE terminal structures [38, 49]. We found that Xenopus Helitron activity also generates protogenes by capturing and duplicating genome-resident gene sequences with 64 Heli1Xen1 insertions containing protein-coding sequences transduplicated from genome-resident genes. This estimation is conservative since it relies on a strict definition of Helitrons and excludes truncated insertions. Moreover, we recorded only well-annotated protein-coding genes from X. laevis and excluded non-coding genes. This number of Helitrons with gene capture events is much lower than what has been reported in the maize (~ 20,000), in the rice (~ 508), in the silkworm (~ 3,546) or the bat (~ 12,382) genomes but only marginally less than in Arabidopsis (~ 216) [39, 50, 51] [52]. Thus, our finding provides another example of the extent of real quantitative differences in the frequency of genome-resident gene capture events by Helitrons between species [53]. Across these different species, the vast majority of Helitron insertions do not contain the whole transcribed unit or all exons from a given gene, and there is a skew toward the capture of the 5’ exons [39, 40]. In addition, the orientation of the gene fragments is strongly biased toward conservation, compared to the orientation of the Helitron [50, 54]. Our findings agree with these features, which are likely the outcome of the underlying molecular mechanisms and the evolution of these insertions.

Gene fragment capture may lead to epigenetic conflicts between TEs and host genes by increasing the sequence promiscuity between the host genes and TE [55, 56]. This promiscuity occurs through the combination in cis of transposon and host sequences, which increases the difficulty of effectively silencing the transposons. According to this model of genomic conflict, acquiring a host gene fragment by a TE would either lead to the silencing of the donor host gene or to a titration of a small RNA-mediated response by the host against the TE. In addition to epigenetic conflicts, gene fragment capture impacts the transcriptome since it can lead to chimeric transcript expression and the birth of new genes. Yet, most Helitron amplified protogenes are evolutionarily neutral, and only some are under purifying selection [39, 57].

We found a stark difference in the Helitron coverage between the three Xenopus genomes studied, with a ten-fold higher abundance in X. laevis compared to X. tropicalis. We are more cautious regarding an interpretation of the result of the X. borealis genome since the assembly contains too many gaps. One line of interpretation of the difference between X. laevis and X. tropicalis is ploidy since X. laevis genome is allotetraploid (2 N = 4X = 36) while X. tropicalis is diploid (2 N = 2X = 20) with an estimated divergence of 48 mya [43]. Two ancestral Xenopus species fused their genomes 17 mya to form the modern X. laevis allotetraploid genome [18]. This genomic shock likely triggered a TE activation wave that could explain the pattern of TE abundance and diversity in X. laevis versus X. tropicalis [58]. It has been hypothesized that the surge of DNA transposon activity has played a role in the subgenome-biased gene loss that occurred after the allopolyploidization event in X. laevis. Still, we did not record a significant bias of Helitrons between the X. laevis L and the S subgenomes [59].

We observed a striking difference in the mean size distribution of non-autonomous insertions of the Heli1Xen2 family (2,067 bp) in comparison with that of the Heli1Xen1 (152 bp) or Heli2Xen (480 bp) families (Fig. 1). Such a difference between non-autonomous Helitrons has been recorded previously with a mean size of 4,616 bp for maize, 441 bp for rice and 950 bp for Arabidopsis elements in association with differences in subfamily structure in these angiosperm species [53]. While the determinants of such differences are elusive, they may reflect different stages of population dynamics or chance events, such as the capture of other TE or tandem repeats.

The recurrent HTT events mediated by Heli1Xen1 in vertebrate genomes highlighted their strong capacity for invasion into organisms via horizontal transfer, as previously reported [14, 36, 38, 60, 61]. Moreover, their activity within Xenopus genomes following invasion underscores their role in driving genome evolution in plants and animals, including vertebrates [38]. Horizontal transfer involving Helitrons has been previously documented in various animals, including vertebrates [62]. Chiropterans are specific hotspots for HTT of Helitrons and display patterns of long-term activity [63]. Our findings are consistent with multiple and independent events of Heli1Xen1 HTT between closely related species (X. laevis and X. borealis; B. bufo and B. gargarizans) but we can not determine whether they are donors or receivers. In another previous large-scale study of HTT in vertebrates, Helitrons were not reported as being involved in HTT between vertebrates, most likely due to the difficulties associated with Helitron annotation [6]. Yet, the pattern of HTT between amphibians, squamates, and ray-finned fishes was frequently found, with DNA transposons contributing to most such events. One possible explanation for this trend could be that viruses and other eukaryotic parasites targeting these sympatric vertebrate species living in the aquatic biomes are major drivers of Helitron’s HTT. Like Baculoviruses acting as cargos to shuttle TE between lepidopterans [64], we hypothesize that Ranaviruses could be agents shuttling TEs as they are large dsDNA viruses capable of infecting these three ectothermic vertebrate classes. In addition, Ranaviruses have broad host specificity and, thus, potential for interspecies and interclass transmission [65]. Since there is a lack of current evidence for this cargo role of Ranaviruses, future experiments could explore this hypothesis as different Frog Virus 3 isolates are used in transmission experiments [66]. A second possible explanation could be that extracellular vesicles would efficiently drive Helitron’s HTT in sympatric aquatic vertebrates. Most amphibians and ray-finned fishes share the life-history traits of external fertilization and development, and the involvement of these traits was suggested in the horizontal gene transfer of a fish antifreeze protein [67]. Early embryogenesis in these species, and especially gastrulation, is associated with cellular movements that trap extracellular material from the perivitelline space into the embryo at a time when the embryo is made of relatively few cells (e.g. 10 000–30 000 in Xenopus) and where germ cell precursors are more easily reachable. In addition, cells undertaking gastrulation movements are known to have increased endocytosis activity [68]. While we did not find published reports of Helitrons sequences in DNA extracted from extracellular vesicles, the findings of encapsulated LINE1 elements in human cells [17] and the involvement of exosomes and other extracellular vesicles in parasite infections [69] warrant future experiments to provide these much wanted experimental pieces of evidence.

Our results open an exciting avenue for the study of active Helitron in a vertebrate animal model amenable to a wide array of experimental analyses at different scales to explore the mechanism and genomic impact of Helitron transposition, similar to what was done on the resurrected Helraiser [70]. We acknowledge that our study's limitation is providing indirect evidence of Heli1Xen1 activity. Future work should provide additional diagnostic indicators of Heli1Xen1 activity, such as full-length RNA transcripts and extrachromosomal circular intermediates [71]. In addition, further experimental work in Xenopus could help tackle questions such as how Helitrons can shape genome structure and induce epigenetic conflict [38]. Already, we provided evidence that Helitron captured several exons. Remarkably, most of the captured gene fragments correspond to the 5’ end of their transcripts, which raises questions about the mechanisms involved. It has previously been reported that such gene captures relied on the activity of LINEs to reverse transcribe and retrotranspose cellular mRNAs into Helitron [39]. Yet, such LINE-mediated reverse transcription products are typically truncated in their 5’ end [72]. Thus, the mechanisms that could explain this 5’ enrichment remain to be deciphered and may involve a different pathway. We also noticed that fragments of seven regulatory genes (dact3, hlx, stcl1, tmem65, pou3f2, hes4 and rgmb) involved in early embryonic development have been captured and shuffled in the X. laevis genome by Heli1Xen1 Helitrons, underscoring their impact on fundamental biological processes. Yet, the physiological consequences need to be explored in more detail.

Conclusion

We concluded that the Heli1Xen1 family is young and might still be active within the Xenopus population. Our data paint a picture of Xenopus Helitrons in which gene capture is rare and with individual capture events distributed marginally by transposition across the genome. This reinforces the notion that Helitrons can significantly affect vertebrate genome evolution. Moreover, the frequent Heli1Xen1-mediated HTT events among amphibians, fishes, and reptiles, which are shared hosts for certain viruses or parasites, should enable us to investigate their roles in acting as vectors for HTT. To advance our understanding of such HTT events, we must complete the evidence presented here by experimental validation.