Identification of telomere-specialized retrotransposons

We searched for non-LTR retrotransposons from the Jockey superfamily in long-read genome assemblies from a total of 106 Drosophila species and three outgroup species of Drosophilidae (Supplemental Table S1). We used an automated pipeline to detect ORFs with homology to known Jockey superfamily ORF1 and ORF2 peptides in each assembly, followed by a phylogenetic approach to assign ORFs to telomeric versus nontelomeric clades (see Methods) (Supplemental Fig. S1). We then used manual curation to generate full-length consensus sequences for each TE family and to confirm that candidate telomeric TEs form head-to-tail arrays at contig ends, consistent with telomere specialization (see Methods) (Supplemental Figs. S1, S2).

We identified a total of 396 putatively telomere-specialized TE families across 109 species (Supplemental Table S2). Of the 396 telomeric TE families, we detected 188 families in head-to-tail arrays at the extreme ends/termini of gene-rich, megabase-long contigs in their host species’ genome assembly. One hundred eighty-nine TE families were found in multiple head-to-tail copies but only on short contigs lacking genes that would allow their assignment to chromosome arms. The remaining 19 TE families were present in a single copy in their host species genome assembly but were supported as telomeric based only on their position in our ORF1 and/or ORF2 gene trees (Supplemental Fig. S3). We were unable to identify Jockey telomeric TEs in a total of 15 species, including the previously reported case in D. biarmipes (Saint-Leandre et al. 2019). In three of these species (Drosophila kurseongensis, Drosophila orena, and an undescribed species from the funebris group, labeled as Drosophila sp. St01m by Kim et al. 2021), we identified fragments of Jockey telomeric TEs lacking intact ORFs, but no full-length TEs (Supplemental Table S3). To confirm the absence of active Jockey telomeric TEs in all 15 species, we analyzed the raw genome sequencing reads and recovered only the same TE fragments present in the genome assemblies of the three species above, confirming the absence of Jockey telomeric TEs in these species. To determine whether other non-Jockey TEs could have become telomere specialized in these species, we manually searched telomeric (based on synteny) scaffold ends in these genome assemblies but were unable to find any TE arrays. Thus, in total, we infer that telomere-specialized TEs were independently lost at least 10 times across the genus (Fig. 1).

Figure 1.

Evolutionary relationships of the species included in this study. Colors correspond to major Drosophila clades, as delineated by Suvorov et al. (2022). The dashed gray branches indicate the 15 species for which no telomeric TEs were identified. The bar plot shows the number of telomeric TE families identified in each species.

Open in new tabDownload PowerPointLink to figure

In the remaining species, we identified about four telomere-specialized families per species, on average (Supplemental Fig. S4). On average, full-length elements were 8.00 kb in size and present in approximately 11 copies per genome, based on read depth (see Methods) (Supplemental Fig. S4). The vast majority of species with telomeric TEs have at least one telomeric TE family that carries both ORF1 and ORF2; however, we identified full-length TEs carrying only ORF1 in the genomes of Zaprionus ghesquierei and Drosophila rufa but only fragments of telomeric TEs with ORF2. Overall, ∼60% of telomeric TE families contained both ORF1 and ORF2, whereas ∼40% contained ORF1 only, having lost ORF2 (Supplemental Fig. S4). We also found 18 TE families in which ORF2 was present but ORF1 was missing, a structure that has not been previously reported. In five of these families, ORF1 was interrupted by a premature stop codon whereas the remainder completely lacked ORF1 sequence (Supplemental Table S2). These TEs have significantly reduced copy numbers compared with TE families with either ORF1 only or both ORF1 and ORF2 (Wilcoxon test P = 0.0004 and P = 0.0022, respectively) (Supplemental Fig. S4). We therefore conclude that these likely represent older fragments of previously active TEs.

To assess the accuracy of our classification of Jockey elements as telomeric versus nontelomeric, we focused on a subset of 30 Drosophila species with chromosome-level genome assemblies. We identified 363 insertions of telomeric clade Jockey TEs across these species, 337 (93%) of which were located within 500 kb of the telomeric end of a chromosome-length scaffold (for discussion of telomeric Jockey elements with locations outside the telomere, see Supplemental Materials). Insertions of TE families belonging to the nontelomeric clade of Jockey elements showed the opposite pattern: Of the 57,760 insertions of nontelomeric Jockey elements across all 30 species, only 321 (0.5%) were located within 500 kb of the telomeric scaffold end (Supplemental Fig. S5). The nontelomeric Jockey families are instead highly enriched in pericentromeric heterochromatin (Supplemental Fig. S5). These results strongly suggest that our phylogenetic classification of telomeric versus nontelomeric Jockey elements is accurate.

Taken together, our results show that the genomes of the vast majority of Drosophila species contain multiple families of telomere-specialized retrotransposons and that the dependence of ORF1-only telomeric TEs on ORF2 proteins encoded by other telomeric TE families is widespread across the genus. The two species for which we were unable to identify any intact telomeric TEs with ORF2 (in either the genome assembly or the raw reads) raise the possibility that some telomeric TEs could be using ORF2 from other nontelomeric Jockey clade TEs for replication.

Evolutionary diversification of telomere-specialized retrotransposons

We used ORF1 and ORF2 peptide sequences to create phylogenetic trees showing the evolutionary relationships among all 396 telomeric TE families. We identify six major clades of telomeric retrotransposons across 109 Drosophila species. These clades are present in both our ORF1 and ORF2 trees and include all original telomeric TE clades described by Villasante et al. (2007): TR1, TR2, TR3, TR4, and HTT, although we note that our more comprehensive trees suggest the TR4 clade is actually part of the HTT clade (Fig. 2A; Supplemental Fig. S6). We therefore use HTT/TR4 here to encompass both groups (see Methods) (Fig. 2A; Supplemental Fig. S6). Our results confirm the monophyly of the Villasante et al. clades and also provide evidence of two novel clades, which we name TR5 and TR6. Many species contain diverse assemblages of telomeric TEs: 94 species have telomeric TEs from more than one TR clade, with 18 species harboring three or more TR clades (Fig. 2A; Supplemental Table S2).

Figure 2.

Evolutionary history of telomeric TE families. (A) Phylogeny of telomeric TE families inferred from ORF2 peptide sequences. Six major telomeric retrotransposon (TR) clades were identified, each with 100% bootstrap support (see Supplemental Fig. S6). The naming scheme is based on TR clades originally described by Villasante et al. (2007), with the addition of the novel TR5 and TR6 clades identified here. The tree branches are colored by TR clade, and the outer circle colors correspond to the host species clades shown in B. The tree is rooted using ORF2 sequences from nontelomeric Jockey clade TEs (data not shown) (see Supplemental Fig. S6). (B) Dated species tree showing major events in the diversification of Drosophila telomeric TEs based on species tree/TE tree reconciliation (see Methods). Rectangles are colored based on the TR clades shown in A. TR clade origins are indicated by filled rectangles, and losses are indicated by an “X.” Horizontal transfer events are summarized by circular (transfers within subgenera) and vertical (transfers between subgenera) arrows. Both types of arrows are colored based on TR clade. Two additional transfer events are shown by curved arrows: one giving rise to the origin of the TR2 clade (black arrow) and the other showing the transfer of a TR2 clade TE to the outgroup species Leucophenga varia (blue arrow). For additional descriptions of TR clade evolution, see Supplemental Text in the Supplemental Materials. Parentheses next to the major Drosophila species clades indicate the number of species harboring the corresponding TR clade over the total number of species within each species clade. (C) Rates of horizontal transfer for nontelomeric versus telomeric Jockey clade TEs.

Open in new tabDownload PowerPointLink to figure

We used our ORF1 tree to investigate the evolutionary relationships among the abundant ORF1-only TEs found across the Drosophila genus. The location of these TEs within our ORF1 tree suggests that ORF2 has been repeatedly lost across dozens of independent TE lineages (Supplemental Fig. S7). However, we also find monophyletic clades containing between two and six ORF1-only TE families, which suggests that not all ORF1-only TEs arise from independent loss of ORF2 (Supplemental Fig. S7). Instead, ancestral ORF1-only TE lineages can birth new ORF1-only TE lineages. The repeated loss of ORF2 from telomeric TEs is particularly striking given that this phenomenon does not seem to occur in nontelomeric Jockey elements (Arkhipova 2012).

We next sought to investigate the evolutionary history of these telomeric TEs across the Drosophila genus. In a mutualistic relationship, in which the symbiont is transmitted vertically from parent to offspring, its evolution will be tightly coupled to that of its host, resulting in a pattern of cospeciation. Such a pattern has been observed in a variety of obligate endosymbionts (Moran and Baumann 1994). Thus, if telomeric TEs are acting as mutualists, the major telomeric TE clades should mirror the major clades found in the Drosophila species tree. We would also expect that the most ancient TE clade should be found broadly, across all extant Drosophila species, because its origin would coincide with or predate the most recent common ancestor of the genus. Instead, we find that many distantly related species harbor closely related telomeric TEs. For example, TR1, TR2, and TR3 clade elements are found in species from both the Sophophora and Drosophila subgenera. Furthermore, the early diverging TE clade, TR1, is found in only six of 10 major species clades, whereas TR2, one of the more derived clades, is found across all 10 species clades (Fig. 2).

There are two evolutionary scenarios that could explain these observations: (1) The divergence of the majority of telomeric TE clades predated the common ancestor of the Drosophila genus, and different TE clades were subsequently lost from different Drosophila lineages, or (2) divergence of the TE clades accompanied the divergence of the genus, followed by frequent horizontal transfer (HT) of TE clades among Drosophila species groups. To assess the evidence for these two scenarios, we performed gene tree/species tree reconciliation (see Methods) to infer duplication (in the case of TEs, this would be divergence of an ancestral TE family into two related families within a single host lineage), transfer, and loss events for the telomeric TEs. To more easily visualize the reconciliation results, we focused on nine major species clades (Suvorov et al. 2022) and pruned the species and TE trees (for the pruned TE tree, see Supplemental Fig. S6) to include the minimum number of species necessary to account for the presence of each TR clade within each species clade (Fig. 2B). We found that the inferred rate of HT exceeded the inferred rates for both duplications and losses (HT, 22; duplications, six; losses, 13). We confirmed that these putatively horizontally transferred TEs showed patterns of sequence similarity consistent with HT (Supplemental Figs. S8, S9) and repeated this analysis on the unpruned species and gene trees, which infers a total of 82 HT events, 12 duplication events, and 70 losses. For comparison, we ran the same analysis pipeline on a set of 465 nontelomeric Jockey TEs recently identified in Drosophila (Tambones et al. 2019). The rate of HT for telomeric TEs is ∼63% higher than that of nontelomeric TEs, a difference that is highly significant (Wilcoxon test P < 2.2 × 10⁻¹⁶) (Fig. 2C). Thus, our reconciliation analysis suggests that the six major telomeric TE clades originated within different Drosophila lineages at different time points during the diversification of the genus and subsequently expanded their host range via HT.

Convergent evolution of telomere specialization

Telomeres can act as a “safe harbor” in which TEs can insert without deleterious consequences. Indeed, transposon insertions have been identified within the canonical telomere repeats in species of fungi, silk moths, and plants (Higashiyama et al. 1997; Fujiwara et al. 2005; Gladyshev and Arkhipova 2007; Rahnama et al. 2020). One might therefore expect other transposons to target their insertions to Drosophila telomeres. To assess this possibility, we asked whether there were Jockey family TEs from the nontelomeric clade that form head to tail arrays at the telomeric regions of any species included in our analysis. We identified Jockey retrotransposons from the nontelomeric clade at the telomeres of six species: Drosophila cardini, Drosophila funebris, Drosophila littoralis, Drosophila virilis, Scaptomyza hsui, and Scaptomyza pallida (Fig. 3). These TEs form a monophyletic clade within the larger clade of nontelomeric Jockey elements, despite their appearance in species that are only distantly related. For example, D. cardini, D. funebris, and S. hsui/pallida diverged from each other more than 20 million years ago (Kumar et al. 2022). The presence of these closely related TEs in distantly related species is indicative of HT (Fig. 3; Supplemental Fig. S10). Hereafter, we refer to this as the NTT clade (nontelomeric clade TEs at telomeres).

Figure 3.

Convergence and neofunctionalization in telomeric TE families. (A) Species tree showing the six species harboring nontelomeric clade TEs at their telomeres (NTTs, red branches) and the four species carrying potentially neofunctionalized ORF2 (blue arrows). The outer circle colors delineate major species clades. (B) Pruned TE tree showing that the NTT families (red branches) form a monophyletic subclade within the larger clade of nontelomeric Jockey elements, inferred from ORF2 peptide sequences. (C) Schematic examples of two species (D. virilis and S. hsui) showing telomeric regions composed of head-to-tail arrays of telomeric TEs and NTTs. In S. hsui, we found a fragmented TE sequence (red arrow) at the telomeric end of the contig whose partial ORF2 sequence placed it within the telomeric clade. (D) Seven telomeric TEs from members of the virilis species group (see A) contain an ORF2 (brown box) that shows no significant homology with the ORF2 (purple box) of other telomeric TE families from this same group via peptide BLAST. Lengths correspond to the full-length TE in base pairs and the ORF2-encoded peptide in amino acids (aa). (E) A profile HMM created from the seven unusual ORF2 sequences shows homology with several regions of the ORF2 from a D. littoralis telomeric TE that also contains the canonical endonuclease and RT domains usually found in Jockey clade ORF2 peptides, plus a domain of unknown function termed X. The unusual ORF2 sequences in E may represent a case of neofunctionalization, in which the ancestral function provided by the endonuclease and RT domains has been lost and replaced by a currently unknown function.

Open in new tabDownload PowerPointLink to figure

These elements are not found elsewhere in the genomes of any of the six species listed above, consistent with telomere specialization; however, in all six species, their arrays are found in between and/or adjacent to arrays of telomeric-clade Jockey TEs (Fig. 3C; Supplemental Fig. S10). Their genomic location suggests that these elements have acquired the ability to target their insertions to the telomeric regions of the genome. Nevertheless, it remains unclear if they are able to act as bona fide telomeres (i.e., by mobilizing directly to chromosome ends) or if they are instead parasitizing the telomeric niche by homing in on pre-existing telomeric TE arrays.

It is also conceivable that other non-Jockey TEs have evolved telomere specialization. To assess this possibility, we calculated the median distance to the telomere for all TE families from 30 Drosophila species with chromosome-level genome assemblies. With the exception of the above NTT elements and telomeric Jockey TEs, no other TE families showed a telomere-biased distribution in their genomic locations, suggesting that NTT and Jockey elements are the only telomere-specialized TEs across these 30 species (Supplemental Fig. S11).

Potential neofunctionalization of ORF2

The ORF2 of all LINE elements, including mammalian LINE-1, contains, at a minimum, an endonuclease and a RT domain (Eickbush and Malik 2002). The domain organization of ORF2 within the Jockey clade of LINEs, which includes all Drosophila telomeric TEs, is similar to that of LINE-1 elements, with the endonuclease domain preceding the RT domain (Eickbush and Malik 2002; Wicker et al. 2007). The presence of a conserved endonuclease domain within the ORF2 carried by Drosophila telomeric TEs is surprising, considering that transposition to chromosome ends should not require DNA nicking (Casacuberta 2017). To better characterize the extent of ORF2 domain conservation across Drosophila telomeric TEs, we used InterPro to identify conserved domains in the ORF2 sequences from our library of telomeric TEs. Seven and nine ORF2 proteins were found to be missing the RT and endonuclease domains, respectively (Supplemental Table S4). Two of these cases were because of a premature stop codon truncating the ORF2 reading frame, whereas the seven remaining ORF2 proteins showed no signs of truncation but were missing both RT and endonuclease domains. All seven ORF2 proteins are from species in the virilis group, and according to InterPro, these unusual ORFs contain no conserved domains at all aside from a single disordered region. The first ORF in these TEs is clearly homologous to the ORF1 protein found in other Drosophila telomeric TEs, which is how these elements were classified as telomeric TEs in the first place. However, a peptide BLAST search was unable to detect any significant amino acid similarity between these unusual ORF2 peptides and the canonical ORF2 peptides from other telomeric TEs (Fig. 3D).

In total, we identified seven telomeric TE families carrying this unusual ORF from the following species: D. virilis, Drosophila americana, Drosophila novamexicana, and D. littoralis. The peptides encoded by these ORFs range in size from 793 to 1219 amino acids, and they range from 30% to 96% identical among the seven TE families (Fig. 3D; Supplemental Fig. S12). We were unable to identify homologs of these peptides outside of the above species, despite employing sensitive hidden Markov model (HMM)–based search strategies (see Methods). We calculated the ratio of nonsynonymous (dN) to synonymous (dS) substitutions for all pairs of sequences and found that, in 20 out of 21 pairwise comparisons, there is statistically significant evidence of purifying selection, which suggests that the ORFs encode functional proteins (Supplemental Table S5). The pair that did not reach significance is highly similar to one another, resulting in a lack of power to detect purifying selection (Supplemental Table S5).

In addition to the conserved endonuclease and RT domains, the ORF2 proteins of some Drosophila telomeric TEs carry an additional glutamine-rich region consisting of approximately 400 amino acids at the C terminus, referred to as the “X domain” (Casacuberta and Pardue 2003). We hypothesized that perhaps the unusual ORF that we discovered was derived from this X domain, having lost both the endonuclease and RT domains. To test this hypothesis, we used a HMM approach (see Methods) to compare the seven unusual ORF2 peptides identified here to all telomeric TE ORF2 peptides containing the endonuclease, RT, and X domains (hereafter referred to as EN/RT/X ORF2 peptides). These searches revealed several regions of homology between a D. littoralis EN/RT/X ORF2 peptide and our HMM (Fig. 3E). Small portions of both the endonuclease and RT domains in the D. littoralis EN/RT/X peptide show weak evidence for homology with our HMM (HMMER E-value = 1.8 × 10⁻⁴ and 0.2, respectively). There is much stronger evidence for homology between the C terminus of the HMM and the X domain of this peptide (HMMER E-value = 2.2 × 10⁻¹⁴). These results suggest that the unusual ORF2 peptides we discovered are in fact derived from the canonical ORF2. However, the endonuclease and RT domain regions have evolved so rapidly that they are no longer recognizable as conserved domains and may no longer retain their ancestral function, raising the possibility that this ORF has acquired a novel function in these species. This represents, to our knowledge, the only known example of a LINE ORF2 peptide that lacks both RT and EN domains.

Gene capture by telomeric transposons

In plants, gene capture may allow TEs to evade host silencing by forcing the host to reduce the efficiency of TE suppression mechanisms in order to avoid self-silencing (Lisch 2009; Muyle et al. 2021). Similarly, our prior work showed that the D. melanogaster telomeric TART element has captured a fragment of the host piRNA pathway gene nxf2 and that TART-derived piRNAs likely target nxf2 for suppression (Ellison et al. 2020). This antisilencing strategy imposes a cost to the host genome while benefitting the TE, consistent with TART acting as a parasite rather than a mutualist. We therefore sought to determine if gene capture occurred in the telomeric transposons of other Drosophila species by searching the full-length consensus of all Drosophila telomeric TEs identified here for homology with D. melanogaster peptides. In Drosophila tristis, we identified multiple fragments of the host gene pangolin embedded within the telomeric TE array of the dot chromosome; however upon further inspection, we concluded that the presence of pangolin coding sequence within the telomere of this species is likely owing to unequal crossing over rather than gene capture by telomeric TEs (Supplemental Fig. S13). After manual curation of all remaining BLAST hits (see Methods), we identified 20 telomeric TE families that show evidence of gene capture (Supplemental Table S6), which represents at least nine independent gene capture events across the genus, including the previously described capture of nxf2 (Fig. 4A). Eight of the nine events involve capture of known piRNA pathway genes: the capture of nxf2 by D. melanogaster plus seven additional capture events involving either piwi or aubergine (aub), both of which are members of the Piwi subfamily of Argonaute proteins that bind piRNAs and are required for transposon silencing (Saito et al. 2006; Vagin et al. 2006; Brennecke et al. 2007; Gunawardane et al. 2007).

Figure 4.

Repeated capture of piRNA pathway genes by telomeric TEs. (A) Species tree showing independent gene capture events (red thick lines) at the nodes. Labels indicate the gene that was captured, and the outer circle is colored by major species clade. (B) Diagrams showing the location of the captured sequence within the full-length TE consensus sequence and within the host gene. Green and purple boxes indicate the locations of ORF1 and ORF2, respectively. TEs lacking a purple box have lost ORF2. Note that the 5′ UTR of TART-A is copied from a portion of the 3′ UTR during replication; thus, the 3′ UTR of TART-A also carries nxf2 sequences (Ellison et al. 2020), but they are not shown here for clarity.

Open in new tabDownload PowerPointLink to figure

We find that a fragment of piwi was captured by a telomeric TE family present in both Drosophila auraria and Drosophila triauraria, and fragments of aub were captured by a telomeric TE family in Drosophila sucinea and nine different species within the vittiger subgroup of the Zaprionus clade. The other novel capture event involves another telomeric TE family from D. auraria and D. triauraria, which carries fragments from a homolog of the F-box domain containing D. melanogaster gene CG12520. F-box proteins are known to provide substrate specificity to the SCF ubiquitin ligase complex (Kipreos and Pagano 2000). CG12520 is expressed in the female germline (Tiwari et al. 2019), and knockdown of CG12520 causes a very modest upregulation of two telomeric TEs (Czech et al. 2013). This gene is largely uncharacterized; for example, its binding substrate remains unknown (Dui et al. 2012). Across all gene capture TE families, the length of the captured gene fragments ranges from 140 bp to 235 bp, and their similarity to their host gene ranges from a percentage identity of 76.4% to 94.82% (Fig. 4B; Supplemental Table S6).

Five of the nine independent gene capture events occurred in the Zaprionus clade. Out of 17 Zaprionus species with sequenced genomes, we find 10 species harboring one or more telomeric TEs that contain aub gene fragments (Fig. 5A). In total, four different regions of the aub gene were captured by different TE lineages within this group (Fig. 5A). Some of these telomeric TEs, such as those from Zaprionus africanus, Zaprionus gabonicus, and Zaprionus indianus, captured two different regions of the aub gene (Fig. 5A). The captured aub fragments were subsequently amplified within these TEs, up to seven copies in the Z. africanus TE and eight copies in Z. gabonicus and Z. indianus TEs (Fig. 5A).

Figure 5.

Multiple captures of aubergine (aub) by telomeric TEs in Zaprionus. (A) Zaprionus species tree showing independent gene capture events. The rectangles are colored based on the aub region that was captured. Diagrams show the location of the captured region within the full-length TE. Note that the region B of aub was captured in the ancestor of Z. vittiger and Z. capensis (solid blue rectangle) and subsequently lost in the Z. davidi/Z. taronus/Z. capensis clade (open blue rectangle), whose TEs now carry region C of aub. (B) Normalized piRNA abundance from the captured host gene (red indicates aub; green, piwi) calculated in 80 bp nonoverlapping windows across the gene transcript. aub is used as a control (blue) in species in which piwi was captured, and piwi is used as a control in species in which aub was captured. The control gene produces significantly fewer piRNAs in all within-species comparisons. Wilcoxon test P < 0.05 in all cases: D. auraria, P = 0.026; D. triauraria, P = 1.3 × 10⁻⁴; Z. davidi, P = 3.0 × 10⁻⁴; Z. gabonicus, P = 6.3 × 10⁻⁴; and Z. indianus, P = 5.2 × 10⁻⁸. (C) Gene tree inferred using aub coding sequence. The clades analyzed in D are indicated by vertical bars. (D) Comparison of rate of evolution of the captured (red) portion of aub versus the uncaptured portion. In the species in which telomeric TEs have captured part of aub (capture clade), the homologous sequence in the aub gene is evolving approximately twice as fast as the uncaptured portion of the gene (Fisher's exact test P = 2.7 × 10⁻⁵). In species in which this region has not been captured, it evolves at a similar rate as the uncaptured portion of the gene (Fisher's exact test P = 0.64).

Open in new tabDownload PowerPointLink to figure

It is possible that these seemingly independent gene capture events actually originated via a single capture event that included a large portion of the aub gene, followed by independent deletions of aub sequence in different TE lineages. Indeed, the Zaprionus TE families that carry aub sequence form a monophyletic clade in the ORF2 TE tree (Supplemental Fig. S14). We tested this prediction by reconstructing the evolutionary history of each capture event (Supplemental Fig. S14). We created gene trees using alignments containing the captured aub-like TE sequences along with their homologous sequences from the aub gene of each Zaprionus species. These trees confirmed our initial inference that five independent capture events occurred within Zaprionus, and allowed us to assign each event to a node in the Zaprionus species tree (Fig. 5A). Unexpectedly, we find that region B of aub was captured in the ancestor of Zaprionus vittiger and Zaprionus capensis and subsequently lost in the Zaprionus davidi/Zaprionus taronus/Z. capensis clade, whose TEs now carry region C of aub (Fig. 5A; Supplemental Fig. S14). We also note that region B of aub was captured independently by two different Zaprionus clades; however, the captured regions overlap but do not share the same boundaries.

The repeated capture of piRNA pathway genes by Drosophila telomeric TEs suggests that acquisition of these gene fragments may increase TE fitness. We previously found evidence that antisense piRNAs produced from the nxf2-like region of TART are capable of targeting the nxf2 host gene for suppression, consistent with an antisilencing strategy to escape host suppression (Ellison et al. 2020). To determine if a similar scenario is occurring in these other species, we performed small RNA sequencing from ovaries for five species in which telomeric TEs had captured either aub or piwi, as well as four related species in which gene capture had not occurred to use as controls. The median identity between captured sequence and host gene is ∼90% across all gene capture TEs (Supplemental Table S6). This identity is noticeably smaller for two gene capture cases: CG12520 in D. auraria/triauraria (∼79% identity) and aub in Zaprionus sp. sepsoides (∼82% identity) (Supplemental Table S6). We find almost a complete lack of piRNAs mapping to the host genes in these latter two cases of elevated divergence (Supplemental Fig. S15), suggesting that the similarity between captured sequence and host gene is too low to direct host gene targeting by TE-derived piRNAs.

In the species with high sequence identity between the captured sequence and the host gene, we found small RNAs from both the gene capture telomeric TEs and the aub (or piwi) host genes whose lengths and 5′ U bias are consistent with piRNAs (Supplemental Figs. S16–S18). Furthermore, in species in which aub (but not piwi) was captured, the abundance of piRNAs from aub is significantly larger than from piwi, whereas the reverse is true for species in which piwi (but not aub) was captured (Fig. 5B). Overall, piRNA abundance from the captured host gene is significantly larger in the species harboring gene capture TEs compared to species whose TEs have not captured any host gene (Wilcoxon test P = 0.0036) (Supplemental Fig. S19), consistent with the host gene being targeted by piRNAs derived from the gene capture region of the telomeric TE. Across all species harboring gene capture TEs, the telomeric TE families that have captured a host gene fragment have significantly higher copy numbers compared with the other noncapture telomeric TEs in the same genome (paired Wilcoxon test, one-sided P = 0.02) (Supplemental Fig. S20), suggesting that gene capture increases TE fitness beyond its role in suppressing the host piRNA pathway (see Discussion).

There are two ways that piRNA biogenesis from the piwi/aub host genes could be initiated: (1) via ping-pong amplification or (2) via a “trigger” piRNA, as has been described for the D. melanogaster row gene, in which phased piRNA production is triggered by a piRNA from the 1360 TE that shares sequence similarity with row (Mohn et al. 2015). piRNA biogenesis from D. melanogaster nxf2 mRNA is more consistent with the latter phenomenon in which a trigger piRNA from the TART telomeric TE induces phased piRNA biogenesis downstream from the region of nxf2 that was captured by TART (Ellison et al. 2020).

To better assess ping-pong and phasing signatures in piwi- and aub-derived piRNAs, we performed additional deep sequencing of piRNAs (see Methods) from four species: D. auraria, D. triauraria, Z. gabonicus, and Z. indianus. For both piwi and aub, we observe sense-strand piRNA production both upstream of and downstream from the region that was captured by a telomeric TE; however, only the downstream piRNAs show a strong signature of phasing (Supplemental Fig. S18). In all four species, antisense piRNAs derived from telomeric TEs align to the captured region of the host gene and overlap with host gene–derived sense piRNAs (Supplemental Fig. S17). In three species (D. auraria, D. triauraria, and Z. gabonicus), the overlap between these sense and antisense piRNAs shows the signature of ping-pong amplification (Supplemental Fig. S17). These results are consistent with ping-pong amplification and/or trigger piRNAs resulting in the biogenesis of phased piRNAs downstream from the captured region of the host genes aub and piwi.

If aub is being targeted for suppression by TE-derived piRNAs, we expect to see accelerated evolution in the corresponding region of aub, which would reduce piRNA targeting by decreasing the similarity between aub and the aub-like sequences within the telomeric TE. We focused on the clade of aub with sequences from the following species: Zaprionus nigranus, Z. taronus, Z. davidi, and Z. capensis, because the telomeric TEs in all of these species have captured the same single region of aub, and the captured sequences have not been amplified within the TE (Fig. 5A). We find that the region of aub that was captured by TEs is evolving approximately twice as fast as the uncaptured portion of the gene, based on a comparison of trees constructed from the captured region of aub versus aub sequence that was not captured by any Zaprionus TEs (Fisher's exact test P = 2.7 × 10⁻⁵; see Methods). As a control, we repeated this analysis on the aub clade from Zaprionus camerounensis, Zaprionus lachaisei, and Z. vittiger, whose TEs have not captured this region of aub. In this case, there was no difference in the rate of evolution (Fisher's exact test P = 0.64) (Fig. 5D).

We also find that aub has undergone two independent duplication events in Zaprionus, once in the Z. capensis lineage and once in Z. africanus (Fig. 5C). This finding alone is remarkable given that a previous survey of 39 Dipteran species (including 22 species of Drosophila) identified only a single duplication event involving aub since its origin more than 150 million years ago, which occurred in stalk-eyed flies (Lewis et al. 2016). Notably, both Zaprionus aub duplications occurred in species with telomeric TEs that have captured an aub fragment (Fig. 5A), raising the possibility that the duplications were selected for as a means to increase aub dosage in response to it being targeted by TE-derived piRNAs.

In addition to duplications, we also find strong support for introgression of aub from Z. taronus to Z. nigranus, likely owing to recent gene flow between these two species (Fig. 5C; Supplemental Figs. S21, S22). This transfer of aub from taronus to nigranus potentially sets the stage for a subsequent invasion of the nigranus genome by taronus telomeric TEs, whose aub-like sequences are 92.26% identical to those of the taronus aub gene. Indeed, two telomeric TEs, both of which carry an aub gene fragment, show evidence of HT from Z. taronus to Z. nigranus (Supplemental Fig. S21). We infer that, based on sequence similarity, the two telomeric TEs initially diverged in Z. taronus, with one family subsequently losing the ORF2 gene (Supplemental Fig. S21). Synonymous divergence between the taronus and nigranus copies of the aub gene (K_s = 0.034) is larger than that of the ORF1 genes for the two TE families (TE-1, K_s = 0; TE-2, K_s = 0.023), suggesting that the transfer of aub from taronus to nigranus occurred first, followed by the transfer of TE-2 and, most recently, the transfer of TE-1 (Supplemental Fig. S21). Notably, these recently transferred TEs (which captured region C of aub) (Fig. 5A) have replaced the ancestral telomeric TEs in Z. nigranus, which carried region B of aub. Thus, different lineages of gene capture TEs may compete with one another to maintain their residence in a given host genome.

Together, these results suggest that capture of piRNA pathway genes by telomeric TEs has occurred repeatedly across the Drosophila genus, likely as a form of counter-defense, in which TEs selfishly target host suppression systems for silencing to increase their own fitness. The accelerated evolution of the region of the host gene that was captured, which was previously described for nxf2 (Ellison et al. 2020) and shown here for aub, suggests there is selection to avoid this targeting by “erasing” the similarity between the captured sequence carried by the TE and the host gene.

Species

Drosophila species used in this study were provided by Matute laboratory (Z. indianus, RCR 17.04; Z. gabonicus, JD'18; Z. davidi, JD'18; Z. nigranus, 18 CAR07 Z03) and ordered from the National Drosophila Species Stock Center at Cornell (D. auraria, strain 14028-0471.00; D. triauraria, strain 14028-0651.00; Z. ghesquierei, strain 50000-2743.00; Zaprionus kolodkinae, strain 50000-2748.00). Z. taronus (strain 50001-1020.00) and Z. camerounensis (strain 50001-1010.01) were ordered from the National Drosophila Species Stock Center at Cornell, but our sequencing data from these species suggested that they were misidentified and that the two stocks are likely to be from unsequenced species related to Zaprionus megalorchis and Z. sepsoides, respectively. We therefore refer to these stocks as Z. sp. megalorchis and Z. sp. sepsoides.

Data acquisition

Out of the 109 long-read genome assemblies used in this study, 106 assemblies were acquired from the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject): Drosophila albomicans (PRJNA630751) (Mai et al. 2020); Drosophila bifasciata (PRJNA565796) (Bracewell et al. 2020); Drosophila innubila (PRJNA524688) (Hill et al. 2019); Drosophila miranda (PRJNA474939) (Mahajan et al. 2018); Drosophila pseudoobscura (PRJNA596268) (Liao et al. 2021); Drosophila serrata (PRJNA355616) (Allen et al. 2017); Drosophila suzukii (PRJNA594550) (Paris et al. 2020); D. triauraria (PRJNA627893) (Torosin et al. 2020); Drosophila mauritiana, Drosophila sechellia, and Drosophila simulans (PRJNA383250) (Chakraborty et al. 2021); Drosophila athabasca, Drosophila lowei, and Drosophila subobscura (PRJNA545704) (Bracewell et al. 2019); and Drosophila ananassae, Drosophila azteca, Drosophila erecta, Drosophila hydei, D. novamexicana, D. orena, Drosophila persimilis, D. virilis, and Scaptodrosophila lebanonensis (PRJNA475270). The remaining 83 NCBI assemblies were obtained from PRJNA675888 (Kim et al. 2021). Three new assemblies were generated in our lab: D. auraria, Z. sp. megalorchis, and Z. sp. sepsoides (see below).

Sequencing and genome assemblies

Using the Monarch HMW DNA extraction kit for tissue, we extracted DNA from approximately 20 females of D. auraria (NEB T3010) and approximately 30 males of Z. sp. megalorchis and Z. sp. sepsoides (NEB T3060S). We used Oxford Nanopore Technologies (ONT) SQK-LSK109 library preparation kit for D. auraria and the ONT SQK-LSK114 library preparation kit for the two Zaprionus species to construct libraries following the ONT ligation sequencing kit protocol. The library of D. auraria was sequenced on a MinION R9.4 flow cell. Each library of two Zaprionus species was sequenced on a MinION R10.4.1 flow cell. Raw signal data were basecalled using the ONT Guppy software package version 4.0.15 for D. auraria and version 6.4.8 for the two Zaprionus species with the default parameters.

We used Flye (version 2.8.1) (Kolmogorov et al. 2019) to assemble the D. auraria genome, followed by polishing with Medaka (version 1.1.3; https://github.com/nanoporetech/medaka) and Nextpolish (version 1.3.1) (Hu et al. 2020). We removed allelic contigs using Purge Haplotigs (version 1.1.2) (Roach et al. 2018). We then performed Hi-C scaffolding of the remaining contigs using the 3D de novo assembly (3D-DNA) pipeline (version 201008). Hi-C data were generated as described by Torosin et al. (2020). The two Zaprionus species genomes were assembled using Flye (version 2.9.2) with the option ‐‐no-alt-contigs and polished using Medaka (version 1.6.0).

Identification of telomeric retrotransposons

We used telomeric retrotransposons of the genus Drosophila available on the Repbase database (TART_DVIR, HeT-A_DYAK, TART-A_DMEL, TART-B_DMEL, TART-C_DMEL, TAHRE_DMEL, and HeT-A_DMEL) as the query to run RepeatProteinMask (RepeatMasker version 4.1.2; http://www.repeatmasker.org) on each species genome with the parameters -noLowSimple -engine ncbi. We clustered similar DNA sequences within species for ORF1 and ORF2 genes separately using CD-HIT-EST (version 4.8.1) (Li and Godzik 2006) with the parameters -d 0 -r 0 -c 0.9 -n 8 -g 1 -T 4 -M 32000. We aligned the DNA sequences of each cluster using MUSCLE (version 3.8.31) (Edgar 2004) with the parameters -maxiters 1 -diags and then ran PILER (version 1.0) (Edgar and Myers 2005) to construct consensus sequences for each cluster, removing consensi <600 bp in length. We used the MACSE pipeline (OMM_MACSE version 10.01) (Ranwez et al. 2018) with the parameters ‐‐no_prefiltering ‐‐no_postfiltering to generate frameshift-aware amino acid translations for ORF1 and ORF2 from each consensus sequence.

We concatenated the ORF1 amino acid sequences generated by the MACSE pipeline from all species, additionally adding known telomeric ORF1 peptides—TART_DVIR, HeT-A_DYAK, TART-A_DMEL, TART-B_DMEL, TART-C_DMEL, TAHRE_DMEL, and HeT-A_DMEL—and the nontelomeric ORF1 sequences—Juan_DMEL, Jockey_DMEL, Doc_DMEL, and F-element_DMEL. We did the same thing for ORF2 (minus the HeT-A peptides because this TE family lacks ORF2) and ran MAFFT (version 7.471) (Katoh and Standley 2013) with ‐‐auto to generate ORF1 and ORF2 amino acid multiple sequence alignments. We trimmed the aligned sequences using ClipKIT (version 1.1.5) (Steenwyk et al. 2020) with the parameters -m kpic-gappy -g 0.5 for ORF1 and -m kpic-gappy -g 0.6 for ORF2. We removed the trimmed sequences with >20% and >10% gaps, for ORF1 and ORF2, respectively. We constructed ORF1 and ORF2 gene trees using IQ-TREE (version 1.6.12) (Nguyen et al. 2015) with the parameters -m TEST -abayes -bb 1000.

We used an unbiased approach to define the telomeric clade for the ORF1 and ORF2 gene trees by calculating branch lengths between each novel ORF1/ORF2 gene and the known telomeric and nontelomeric ORF1/ORF2 genes included in our unrooted trees (see GitHub for code). The novel ORF1 and ORF2 genes that showed shorter branch lengths to their known telomeric homologs (compared with their nontelomeric homologs) were assigned as putative telomeric clade peptides. This approach resulted in a single monophyletic clade of putative telomeric ORF1 and ORF2 genes (Supplemental Fig. S3).

We next sought to identify the genomic locations of each candidate telomeric ORF1 and ORF2 gene. The telomeric clade ORF1 and ORF2 consensus DNA sequences from each species were used as a custom library when running RepeatMasker (version 4.1.2) with the parameters -no_is -norna -nolow -engine ncbi, retaining matches with ≤3% divergence from their consensus. We used the RepeatMasker hits as starting points to generate manually curated, full-length telomeric TE consensi. We identified arrays of RepeatMasker hits at contig ends (or within short unplaced contigs) and generated a dot plot of the array sequence using YASS, a genomic similarity search web tool (Noé and Kucherov 2005). We manually examined the dot plots to identify the boundaries of the monomers that compose the head-to-tail TE arrays. We then compared the monomers from the same TE family to identify the longest, unfragmented, representative sequence for each TE family. For species with ONT data, we aligned the nanopore genomic sequencing data to the representative TE sequence and performed one round of polishing using Medaka (version 1.6.0). All consensi (except those from species lacking Illumina data) (see Supplemental Table S7) were then further iteratively polished using Pilon (version 1.24) (Walker et al. 2014) with Illumina data until no new changes were made to the consensus. We searched each consensus for intact ORF1 and ORF2 open reading frames and performed an additional round of polishing using Pilon with Illumina data (when available) (see Supplemental Table S7) when the consensus sequence contained fragmented ORFs.

All manually curated telomeric TE consensi include both 5′ and 3′ UTRs. Recently derived ORF2-defective, nonautonomous elements share the same 5′ and 3′ UTRs with their autonomous progenitors. Our approach classifies the autonomous and nonautonomous versions of these elements as different TE families (Supplemental Table S2). In the cases in which there was no autonomous progenitor present in the genome, we defined the nonautonomous TE boundaries by aligning multiple copies of the same TE, requiring that multiple copies support the same 5′ and 3′ boundaries.

To estimate the genomic copy number of each telomeric TE family, we aligned the long-read genomic sequencing data from the host species to the TE consensi. To avoid underestimation of copy number owing to 5′ end truncation of telomeric TEs, we calculated the median per-base read depth across a 1 kb interval of the TE consensus beginning 500 bp from the 3′ end of the sequence, similar to the approach described by McGurk et al. (2021). We then estimated copy number by dividing this value by the median whole-genome read depth (averaged across 100 bp windows).

TR clade definitions

We defined TR clades so that they would be as consistent as possible with the original Villasante designations while also being monophyletic (Supplemental Fig. S5). For TR1, Villasante et al. (2007) identified TR1 elements in species from the obscura, virilis, and repleta groups. In our ORF2 tree, these sequences form a clade that is closely related to two other clades: one composed of TRs from species in the Zaprionus, Scaptomyza, and cardini groups and the other containing species from the obscura group. The virilis, repleta, Zaprionus, Scaptomyza, and cardini group species are all members of the Drosophila subgenus; thus, we expanded the TR1 clade to encompass the related TRs from these groups. For TR2, Villasante et al. (2007) identified TR2 elements in species from the ananassae, obscura, virilis, and Hawaiian Drosophila groups. Our TR2 clade represents the smallest possible monophyletic clade that includes these elements. For TR3, Villasante et al. (2007) identified TR3 elements in species from the obscura and repleta species groups. Our TR3 clade represents the smallest possible monophyletic clade that includes these elements. For TR4 and HTT, Villasante et al. (2007) identified TR4 elements in species from the obscura group. These elements are closely related to TRs that we identified in species from the montium group, however expanding the TR4 clade to include the montium group elements would result in a clade that is not monophyletic. We therefore combined both HTT and TR4 elements into a single monophyletic clade, which we refer to as “HTT/TR4.” TR5 and TR6 are novel monophyletic clades that were not previously identified by Villasante et al.

Identification of the telomeric retrotransposons with host gene capture

We used D. melanogaster amino acid sequences (longest isoform per gene, r6.42 annotation) as queries and the DNA consensus sequences of all species’ telomeric transposons as the database to run TBLASTN (version 2.10.1) (Camacho et al. 2009) with the parameters -outfmt 6 -evalue 1 × 10⁻⁵. For each TBLASTN hit, we determined whether there was significant similarity between the telomeric TE and its host species gene at the DNA level by running BLASTN (version 2.10.1) with the parameters -outfmt 6 -evalue 1 × 10⁻³. To search for older gene capture events, we also ran BLASTN with the parameters -outfmt 6 -word_size 4 to identify host gene/TE sequences with higher levels of sequence divergence.

Phylogenetic reconstruction (species, telomeric retrotransposons, and their reconciliation)

We used full-length ORF1 and ORF2 sequences from our curated telomeric TE consensi to create telomeric TE gene trees. We additionally added the following known nontelomeric Jockey clade ORFs from Repbase TE families to use as outgroup sequences: Jockey_DMEL, Jockey-1_DEr, Jockey-1_DEu, Jockey-1_DF, Jockey-1_DGri, Jockey-1_DK, Jockey-1_DRh, Jockey-1_DT, Jockey-1_DVi, Jockey-1_DWi, Jockey-1_DYa, Jockey-10_DAn, Jockey-10_DRh, Jockey-11_DAn, Jockey-11B_DBp, Jockey-14_DBp, Jockey-2_DEl, Jockey-2_DEu, Jockey-2_DRh, Jockey-2_DT, Jockey-3_DEu, Jockey-3_DF, Jockey-3_DK, Jockey-3_DRh, Jockey-3_DVi, Jockey-4_DEl, Jockey-4_DEu, Jockey-4_DF, Jockey-4_DPer, Jockey-4_DVi, Jockey-5_DAn, Jockey-5_DBp, Jockey-5_DTa, Jockey-6_DAn, Jockey-6_DEl, Jockey-6_DEu, Jockey-6_DF, Jockey-7_DAn, Jockey-7_DF, and Jockey-8_DAn. We aligned the ORF peptides using MAFFT (version 7.471) and trimmed the alignments using ClipKIT (version 1.1.5) with the parameters -m kpic-gappy -g 0.1. We removed sequences with >20% gaps from each ORF alignment and then created the telomeric ORF1 and ORF2 gene trees using IQ-TREE (version 1.6.12) with the parameters -m TEST -abayes -bb 1000 and rooted the trees using the nontelomeric outgroup sequences. All tree figures were created using the iTOL tool (Letunic and Bork 2024).

To reconcile the above TE ORF1 and ORF2 trees with the Drosophila species tree, we first pruned the dated species tree from (Suvorov et al. 2022), keeping only the species studied here. We then ran Ranger-DTL-Dated (version 2.0) (Bansal et al. 2018) with the default parameters to perform the gene tree/species tree reconciliation. We calculated support values for each reconciliation event using the AggregateRanger program on the output of 1000 Ranger-DTL-Dated reconciliation runs (Supplemental Table S8). We used a single consensus sequence for each TE family as input for reconciliation. Note that Ranger-DTL identifies TE “duplications” as cases in which an ancestral TE lineage gives rise to two novel TE lineages, which is akin to “speciation” from the point of view of the TE. Ranger-DTL compares the TE tree to the species tree and identifies points of incongruence between the two. It then finds the most parsimonious combination of duplication (i.e., TE “speciation”), HT, and loss (i.e., TE extinction) events to reconcile these points of incongruence using a reconciliation cost for each of these events. HT is assigned the highest cost: Transfer events are assigned a cost that is 1.5-fold higher than the duplication cost and threefold higher than the loss cost, which makes our approach conservative in terms of identifying HT events. To assess the accuracy of HT events identified by Ranger-DTL, we used BLAST to create local alignments between pairs of TE sequences and recorded the alignment bit score, which is a normalized measure of both the percentage identity and alignment length of a given BLAST match. For all Ranger-DTL events involving HTs between two tips of the species tree, we recorded the best alignment bit score for the BLAST search between the donor telomeric TE lineage and the recipient lineage. As a control, we used bit scores from related, vertically inherited, nontelomeric Jockey elements between the same species pair.

We used the same Ranger-DTL approach to infer HT events from nontelomeric Jockey sequences obtained from Tambones et al. (2019). The total number of transfer events for each group (i.e., telomeric vs. nontelomeric) was then normalized by the total number of nodes in the corresponding tree. The box plots in Figure 2C show the distribution of these estimates across 1000 runs of Ranger-DTL.

We generated two novel de novo genome assemblies for this study: Z. sp. megalorchis and Z. sp. sepsoides. We placed these species within the Zaprionus species tree using the same methods as described by Suvorov et al. (2022). Briefly, we identified single-copy orthologs from each Zaprionus species genome assembly by running BUSCO (version 5.4.7) (Manni et al. 2021) with the Diptera data set of the OrthoDB v10 release (diptera_odb10) (Kriventseva et al. 2019) using the default parameters. We used IQ-TREE (version 1.6.12) to obtain the gene trees for each single-copy ortholog, as described by Suvorov et al. (2022). The gene trees were used to infer the Zaprionus species tree using ASTRAL (version 5.7.8) (Zhang et al. 2018) with the default parameters.

To make the aub gene tree for Zaprionus species, we first built aub gene models in each species using genBlastG (She et al. 2011) with the parameters -e 1e-5 -g T -r 3 -c 0.6 and the D. melanogaster Aub peptide sequence as a query. We then used MAFFT (version 7.5.20) with –auto to align the aub coding sequences and IQ-TREE (version 1.6.12) with the parameters -m TEST -bb 1000 to infer the gene trees.

To test for accelerated evolution of the aub gene in Zaprionus, we first identified all regions of aub coding sequence that were not captured by any Zaprionus species (i.e., uncaptured sequence) and concatenated these regions together for each aub gene in each species. We aligned these sequences using MAFFT (version 7.520) with –auto. We then extracted two sequence subsets from this alignment: those from Z. nigranus, Z. taronus, Z. davidi, and Z. capensis (capture clade) and those from Z. camerounensis, Z. lachaisei, and Z. vittiger (control clade). We added Z. inermis aub sequence to each sequence subset as an outgroup and used ModelFinder from IQ-TREE, which identified “TPM2+G4” as the best-fitting DNA substitution model for both sequence subsets. We then inferred gene trees for both the capture and control sequence subsets from the uncaptured aub sequence using IQ-TREE, specifying Z. inermis as the outgroup and “TPM2+G4” as the substitution model, and summed the branch lengths of the resulting trees.

We then repeated the above procedure using Region C of aub, which was captured by the Z. nigranus, Z. taronus, Z. davidi, and Z. capensis telomeric TEs, providing the corresponding uncaptured tree file as a constraint tree topology. We then compared the summed branch lengths between the captured and uncaptured aub regions for both the capture clade and the control clade. We used the Fisher's exact test to compare the proportion of substitutions among all sites between captured and uncaptured sequence subsets.

To test for introgression between Z. taronus and Z. nigranus, we downloaded gene trees from BUSCO single-copy orthologs generated in a previous study (Suvorov et al. 2022). Across all approximately 2500 trees, we counted the number of times we observed the following species pairs as sister taxa, which all share the same most recent common ancestor, according to the species tree by Suvorov et al. (2022): Z. taronus versus Z. nigranus (n = 71), Z. davidi versus Z. nigranus (n = 9), Z. capensis versus Z. nigranus (n = 12), Z. taronus versus Z. vittiger (n = 11), Z. taronus versus Z. camerounensis (n = 6), and Z. taronus versus Z. lachaisei (n = 8). Overall, 2.77% of gene trees show Z. taronus and Z. nigranus as sister taxa compared with an average of 0.35% support for the other pairs shown above.

Small RNA sequencing

For seven Zaprionus species plus D. auraria and D. triauraria, we extracted small RNAs from 10 to 24 pairs of ovaries using the TraPR small RNA isolation and library prep kit (Lexogen 135.08) (Grentzinger et al. 2020), following the kit protocol. The small RNA libraries were sent to Novogene for NovaSeq SE50 sequencing.

We subsequently performed deep sequencing of sodium periodate treated small RNAs, with 2S rRNA blocking using a terminator oligo (Wickersheim and Blumenstiel 2013), from four species whose telomeric TEs had captured the host genes aub or piwi: D. auraria, D. triauraria, Z. gabonicus, and Z. indianus. We unfortunately lost our Z. davidii stock after the initial small RNA library was generated and were unable to obtain a replacement. For protocol, see the Supplemental Methods.

We used Trim Galore! (version 0.6.10; https://github.com/FelixKrueger/TrimGalore) to trim adapter sequences from raw small RNA reads and SortMeRNA (version 4.3.6) (Kopylova et al. 2012) to remove reads arising from tRNAs and rRNAs. TRNA genes were predicted by tRNAscan-SE (version 2.0.11) (Chan et al. 2021), and rRNA databases were generated by BLASTN searches with D. melanogaster rRNA queries. The filtered small RNA reads were aligned using ShortStack (version 3.8.5) (Axtell 2013) with the parameters ‐‐dicermin 20 ‐‐dicermax 35 ‐‐mismatches 2.

Detection of piRNA reads mapped to telomeric TEs and host genes, piwi, and aub

We sought to quantify piRNA abundance from each telomeric TE family, as well as the host genes that were captured by telomeric TEs, excluding the captured region of the host gene to avoid cross-mapping of TE-derived piRNAs to the host gene. To do this, we masked the telomeric TEs and host genes (i.e., piwi, aub, and CG12520) in the appropriate species genome using RepeatMasker (version 4.1.2) with the parameters -no_is -norna -nolow -pa 4 -engine ncbi. We then created a custom reference genome by appending consensus sequences of each telomeric TE family as well as mRNA sequences of the longest isoform of the captured host gene(s), in which the captured region was masked using BEDTools, version 2.25.0 (Quinlan and Hall 2010), to the masked genome assembly. Next, we used ShortStack (version 3.8.5) with the parameters ‐‐dicermin 20 ‐‐dicermax 35 ‐‐mismatches 2 to align the small RNA reads to the custom reference genome. We retained reads with lengths between 23 bp and 30 bp and used BEDTools (version 2.25.0) to calculate coverage of sense and antisense alignments on the telomeric TEs and the relevant host transcript(s) (i.e., piwi, aub, and/or CG12520).

Phasing and ping-pong signatures

To assess signatures of phasing, we used sense-strand piRNA alignments from three regions of piwi/aub: upstream of the gene capture region, downstream, and the full gene, excluding the capture region. We collapsed exact duplicates, and for each 3′ end, we calculated the downstream distance to the nearest piRNA 5′ end. For each distance between 0 and 20 bp, we calculated the fraction of 3′–5′ piRNAs pairs whose ends were separated by that distance.

To assess ping-pong signatures, we used SAMtools (Danecek et al. 2021) to extract TE-derived antisense piRNAs and piwi/aub-derived sense piRNAs, from alignments to the unmasked captured host gene and unmasked gene capture TE generated using ShortStack (version 3.8.5; ‐‐dicermin 20 ‐‐dicermax 35 ‐‐mismatches 2). We then aligned both TE-derived antisense piRNAs and piwi/aub-derived sense piRNAs to the unmasked piwi/aub host gene and used PingPongPro (Uhrig and Klein 2019) to identify ping-pong signatures.

d_N/d_S (K_a/K_s) analysis

We used MEGA X (Kumar et al. 2018) to align coding sequences and used KaKs_Calculator 2.0 (Wang et al. 2010) to calculate the ratio of the number of nonsynonymous substitutions per nonsynonymous site to the number of synonymous substitutions per synonymous site. The Fisher's exact test was used (within KaKs_Calculator) to test the null hypothesis of neutral evolution (i.e., equal rates of synonymous vs. nonsynonymous substitutions).

HMM searches

We initially used an alignment of all seven unusual ORF2 sequences as a query to perform sensitive HMM–HMM sequence searches against the UniRef30_2023_02 database on the MPI Bioinformatics Toolkit webserver (Zimmermann et al. 2018) using HHblits (Remmert et al. 2012) but found no significant hits outside of our query species. We then aligned the seven unusual ORF2 amino acid sequences using MAFFT (version 7.520) and created an HMM profile from the alignment using hmmbuild from HMMER3 (version 3.4) (Eddy 2011). We compiled a database of all ORF2 peptides in D. virilis, D. americana, D. novamexicana, and D. littoralis telomeric TEs that contained the endonuclease, RT, and X domains and searched this database with our HMM profile as the query using hmmsearch from HMMER3.