Transposable elements contribute to the evolution of host shift–related genes in cactophilic Drosophila species

Genome assemblies and gene annotation

To investigate the potential role of TEs in host shift of cactophilic Drosophila species, we performed Oxford Nanopore Technologies (ONT) long-read sequencing on species that have different preferential cacti as hosts: D. buzzatii (Opuntia sp.), D. koepferae (columnar cacti), D. arizonae (Opuntia sp. or columnar cacti), D. m. mojavensis (F. cylindraceus), D. m. wrigleyi (Opuntia sp.), and D. m. sonorensis (S. thurberi). We obtained high-resolution assemblies for all genomes, with an average of 567 scaffolds, 13.8 Mb of N50 (Supplemental Table S1), and ∼98% of benchmarking universal single-copy orthologs (BUSCO) genes (Supplemental Fig. S1). Furthermore, the gene annotation for these genomes covered 98.63% of all reference genes in the D. mojavensis subspecies and 95.98% in D. arizonae. The de novo gene annotation in D. buzzatii and D. koepferae revealed a total of 18,050 and 17,848 coding genes, for which 63.77% and 64.13% were successfully assigned as one-to-one orthologs with D. mojavensis, respectively. These results represent the higher gene repertoire identified for the D. buzzatii and D. koepferae species, altogether with their ortholog relationship with D. mojavensis and D. arizonae. Because the analysis of HLAU genes is the focus of this study, their efficient annotation is relevant. They are represented by nine gene families associated with the three steps in the use of host resources: (1) localization—OBPs and ORs, (2) acceptance—GRs and IRs, and (3) host usage—ABCs, ESTs, GSTs, UDPs, and CYPs (Supplemental Table S2). In the D. arizonae and D. mojavensis subspecies, we recovered ∼98.21% of the total HLAU genes from the D. mojavensis reference genome. In D. buzzatii and D. koepferae, the number of annotated HLAU genes covered ∼95% in both species compared with D. mojavensis (Table 1).

Minor effects of positive selection in HLAU genes

The use of Opuntia sp. has been proposed as the ancestral state in cactophilic Drosophila species, whereas columnar cacti are the derived state (Oliveira et al. 2012). To investigate whether species that prefer columnar cacti as hosts exhibit signs of positive selection in HLAU genes, we employed a branch-site approach to compare their evolutionary rates with those of species that use Opuntia sp. We also performed selection analysis with branch and site approaches separately. These analyses were assessed without D. arizonae due to its nonpreferential usage between columnar and Opuntia sp. cacti (Oliveira et al. 2012). Moreover, we added Drosophila navojoa (publicly available from the NCBI Genomes database (https://www.ncbi.nlm.nih.gov/home/genomes/) under accession number GCF_001654015.2), which is a basal species of the mojavensis cluster that uses solely Opuntia sp. (Oliveira et al. 2012). Thus, we carried out the branch-site selection test with D. buzzatii, D. m. wrigleyi, and D. navojoa as Opuntia-related branches (background), and D. m. mojavensis, D. sonorensis, and D. koepferae as columnar-related branches (foreground). It is important to highlight that because of the requirement for one-to-one orthologous genes between the six genomes, our analysis reduced the total number of HLAU genes to 42.9% of the initial number presented in Table 1.

The site model revealed two genes with sites under positive selection: Atet (LOC26527843, FDR = 0.0215) with seven sites with ω > 3 and Cyp4ae1 (LOC6585222, FDR = 0.0394) with two sites with ω > 4 (Supplemental Table S3). This result suggests that both genes have sites under diversifying selection across the six cactophilic species. Still, it does not test the hypothesis of positive selection between species with preference columnar cacti versus Opuntia sp. Thus, we applied the branch model with columnar cacti dwellers labeled as foreground. This approach allows us to test the evolutionary rate between species with different host preferences. The analysis revealed that 96.7% of the genes exhibit different selective regimes on the foreground branch compared with the background (FDR < 0.05). Although the test supports a significant shift in ω, no cases of positive selection (ω > 1 in the foreground) were reported (Supplemental Table S3). The average ω in the foreground branch (columnar cacti dwellers) is 0.13, whereas in the background it is 0.11. Some dramatic changes have been observed, as in the Or92a (LOC6585146) for instance. This gene has ω = 0.0001 in the foreground and ω = 0.02694 in the background (FDR = 2.15 × 10⁻⁷) (Supplemental Table S3), suggesting it has evolved under strong purifying selection. Altogether, these results indicate that HLAU genes are under purifying selection in columnar cacti dwellers. Given that the branch model does not differentiate among codons, these results reflect a global change in selective pressure across the entire gene in the foreground. Then, we also investigated whether specific sites could be under positive selection with the branch-site model, but we did not observe any significant difference (FDR < 0.05) between the foreground and background branches (Supplemental Table S3). This result is consistent with the branch model, in which we did not obtain any signature of positive selection.

Evolutionary dynamics of TEs in cactophilic Drosophila species

We used a combination of automatic and manual methods to build a curated TE library for each species based on the classification of TEs (Wicker et al. 2007). The method was designed to perform polishing steps in both TE consensus and copies to remove artifacts (see Methods). The total number of TE consensuses (LTRs separated from internal sequences) ranged from 255 in D. buzzatii to 372 in D. m. mojavensis (available in Zenodo repository: https://zenodo.org/records/16501018). The TE content in the genomes ranged from 15.04% in D. buzzatii to 22.01% in D. m. mojavensis (Fig. 1A). In addition, we observed a positive correlation between TE content and genome size estimated from the assemblies (Pearson's correlation: P = 0.0082; R = 0.93) (Supplemental Fig. S2). Both species from the buzzatii cluster had the lower genome sizes and TE contents compared with species of the mojavensis cluster (Fig. 1A). Regarding the proportional load of TE orders in the genomes, LTRs and Helitrons contribute on average with ∼70% of total TE content in all species (Fig. 1B).

Figure 1.

TE annotation in the six genomes of the cactophilic Drosophila species. (A) Genome size of the six cactophilic Drosophila species sequenced in this work, and their respective TE content (blue). (B) Proportional distribution of the total TE content in the genomes. (C) Phylogenetic relationship among species (adapted from Matzkin 2008; Oliveira et al. 2012), alongside TE order landscapes representing total copy number (y-axis) and the divergence of TE insertions compared with the consensuses (x-axis), based on the Kimura two-parameter (K2P) model corrected with CpG. The copies located on the left of the graph have low divergence with their respective consensus, hence inferring that they might be conserved/recent/active, whereas copies on the right side represent copies with high divergence owing to the accumulation of mutations, corresponding to old/inactive copies.

Open in new tabDownload PowerPointLink to figure

In the mojavensis cluster (D. arizonae and the three D. mojavensis subspecies), the abundance of TEs in the genome varied substantially (Fig. 1B), despite their recent divergence time of 1.5 Mya (Sanchez-Flores et al. 2016). D. m. mojavensis has the TE-richest genome (22.01%), with the abundance following from the highest to the lowest: LTR > Helitron > LINE > DNA (Fig. 1B). The distribution of TEs in the other two D. mojavensis subspecies and D. buzzatii does not follow such order. LTRs and Helitrons have the highest abundance, followed by DNA and LINEs (Fig. 1B). In D. arizonae and D. koepferae, Helitrons are the most abundant, followed by LTR, DNA, and LINEs (Fig. 1B). The observed differences are often associated with the divergence in the evolutionary patterns of the TE orders between the genomes. We used copy number as a proxy to measure TE expansion in the genomes, and divergence Kimura two-parameter as age measurement. In the mojavensis cluster, we observed that LTRs and Helitrons had a higher expansion in D. m. mojavensis compared with the other genomes (Fig. 1C). The higher copy number in D. m. mojavensis is mostly from young TEs (K2P = 0–5), which might explain the higher TE content in this species. DNA elements have a constant evolutionary dynamics in the genomes, except in D. m. wrigleyi. In this species, we observed a peak in copy number with K2P between 15 and 20, considered old insertions. We analyzed whether these insertions were distributed across the genome or had local accumulation. The result revealed an even distribution across the longest scaffolds of the genome (Supplemental Fig. S3), following similar patterns of both younger and older insertions. LINEs had a similar evolutionary landscape in D. mojavensis subspecies, with a higher copy number from young copies compared with D. arizonae. In the buzzatii cluster, Helitrons show higher copy number from young insertions in D. koepferae than in D. buzzatii (Fig. 1C), whereas LINEs show a past transposition burst in D. buzzatii than in D. koepferae. The antagonistic patterns of TE dynamics observed in the six genomes suggest that TEs had different successful lineages across these species.

HLAU genes are enriched in TEs within regulatory regions

TE insertions may be selected in the promoter of genes if they provide a benefit to the individual fitness, through the spreading of histone marks, or cis-regulatory elements (Daborn et al. 2002; Guio et al. 2018). To investigate whether TEs are enriched in the promoters (2 kb upstream of the transcription start site [TSS]) from the nine HLAU families, we performed permutation tests in the six genomes (see Methods). To compare our findings with other non-host shift–related genes, we also included zinc-finger (ZF) and serine/threonine kinases (SER) as two non-HLAU gene families. We observed that ORs have enrichment of TEs at the promoters in all genomes from species of the mojavensis cluster (P-value < 0.05) (Fig. 2A). For acceptance genes, D. koepferae had enrichment of TEs in the promoters of GRs (Fig. 2A). Finally, the enrichment of TEs within promoters of detoxification genes was observed in D. arizonae (CYPs), D. buzzatii (GSTs), and D. koepferae (ABCs). The results depicting the histograms with frequency of genes with TEs for each significant permutation test, and P-values are shown in the Supplemental Figure S4 and Supplemental Table S4. Overall, detoxification genes had z-scores toward depletion of TEs on promoters, with significant depletion for UGTs in D. buzzatii (P-value = 0.033) and D. arizonae (P-value = 0.0135) and for GSTs in D. koepferae (P-value = 0.0205) (Fig. 2A). In the non-HLAU gene families, we observed the expected result by chance: either the z-score indicated the expected number of genes with TEs on promoters in a genome-wide comparison (z-score near zero), or it indicated depletion of TEs, as observed for SER genes in D. buzzatii. This result reinforces our findings, indicating the significant and evolutionary unexpected enrichment of TEs in HLAU genes, especially on the promoters of ORs.

Figure 2.

Enrichment of TEs in HLAU genes. (A) Enrichment of TEs in the promoter region (2 kb upstream) of the HLAU gene families. Z-score > 0 indicates enrichment, and Z-score < 0 indicates depletion. Red dots indicate P-value < 0.05; gray dots, P-value > 0.05. (B) The same enrichment analysis as depicted in A, but for non-HLAU gene families, represented by zinc-finger proteins (ZF) and serine/threonine kinases (SER). (C) Proportion of TE orders located upstream and downstream ORs. (D) Transcription factor binding site (TFBS) frequency on TEs located upstream of OR genes in all species under study. (E) Proportional frequency of TE families with TFBSs embedded in their sequences from TE copies located upstream of OR genes. (F) Significant Fisher's exact test between OR genes with OR-related TFBSs from TE insertions compared with the frequency observed in the non-HLAU gene families SER and ZF. (G) Proportion of TE orders carrying OR-related TFBSs in non-HLAU gene families.

Open in new tabDownload PowerPointLink to figure

To obtain further insights into the potential regulatory role of TEs, we also analyzed the enrichment at the 2 kb downstream region. Among all species, the only significant enrichments were observed in ORs from species of the mojavensis cluster (Supplemental Table S4), as we observed in the promoters. The TE enrichment upstream of and downstream from ORs indicates a possible coevolution of this gene family with TEs. Alternatively, it may also suggest that ORs are located in TE-rich regions of the genomes. To investigate this hypothesis, we identified TE-rich regions (Supplemental Fig. S5) and analyzed the frequency of ORs. We observed a single OR in TE-rich regions in D. arizonae, D. m. mojavensis, and D. m. wrigleyi. (Supplemental Table S5). Therefore, we rejected the hypothesis that ORs are enriched with TEs in their vicinity owing to their location in TE-rich regions. Other factors might be involved in the evolutionary maintenance of TEs near ORs, such as the co-option of regulatory motifs.

A few TE families are more prone to contribute to cis-elements than others, such as retroelements in D. melanogaster (Chang et al. 2022). Therefore, we investigated whether the observed enrichment of TEs in the HLAU gene families is TE specific. Considering all species, 45.36% and 50.36% of the insertions located upstream of and downstream from, respectively, ORs correspond to Helitrons (Fig. 2C). The frequency of Helitrons upstream ORs is significantly higher than expected given their genomic distribution in species of the mojavensis cluster (Supplemental Table S6). We decided to analyze whether these Helitron insertions upstream of OR genes could contain gene fragments. Helitrons carrying host gene exons have been reported previously in other species than Drosophila, resulting in a high copy number of insertions overlapping with genes (Pritham and Feschotte 2007; Thomas et al. 2014). By intersecting exons from genes enriched with their respective upstream Helitrons, we observed that Or85c (overlap, 772 nt) and Or22c (52 nt) had overlap with Helitrons in D. m. mojavensis, whereas D. m. sonorensis had only Or22c (58 nt). Although the insertion in Or85c suggests a subspecies-specific insertion, the insertion in Or22c is likely to be inherited from the species’ common ancestor (Supplemental Table S7). At the ORs’ downstream region, we observed four genes with Helitrons overlapping exons: Or71a, Or10a, Or13a, and Or82a. The same Helitron in the Or82a was observed in the three D. mojavensis subspecies, highlighting the inheritance from the common ancestor and maintenance of the insertion over time. Taken together, our results rejected the hypothesis of spread Helitron copies with OR fragments, which could explain the high frequency of these elements upstream of and downstream from ORs. However, these elements could still act on gene regulation if they carry cis-elements as transcription factor binding site (TFBS).

TEs located near genes can modulate gene expression by donating TFBSs or poly(A) signals to the ancestral regulatory motifs. To investigate the potential functional role of the observed enrichment of TEs in ORs, we performed TFBS identification on all TE insertions located in OR promoters. We used the transcription factors (TFs) previously described to participate in OR's transcription in D. melanogaster: acj6, onecut, xbp1, fer1, and zf30C (Jafari et al. 2012). Taking all species together, the most frequent TFBS found was for the TF onecut (46.15%), followed by zf30C (33.84%) and Xbp1 (7.69%) (Fig. 2D). Helitrons were the TEs with the highest frequency, representing 50% of all insertions carrying ORs’ TFBSs (Fig. 2E). Our analysis revealed that 39.28% of the ORs with TEs in promoters have the TFBS located within insertions (Supplemental Table S8). To test whether these results are specific to ORs, we also performed the same detection of TFBSs and calculation of TE frequency in the non-HLAU families ZF and SER. Our results show that ZF and SER have, on average, 15.55% and 17.14% of the genes with OR-related TFBSs, respectively. The detection of these TFBSs in ZF and SER is likely owing to chance, given the short length of the motifs. Comparing these proportions with the ones found in ORs, we observed significantly more TFBSs compared with ZF (Fisher's exact test: P-value = 0.0004341) and SER (Fisher's exact test: P-value = 0.0008519) (Fig. 2F; Supplemental Table S9). In addition, Helitrons also had higher prevalence than other TEs upstream of SER and ZF genes (Fig. 2G), demonstrating a potential insertion preference for promoter regions. Although similar Helitron frequencies between OR, SER, and ZF genes, our data suggest that Helitrons might have been retained in OR promoters owing to the presence of functional TFBSs (Fig. 2F).

HLAU genes are differentially expressed between Opuntia sp. and columnar cacti dwellers

Differential expression of genes has been associated with adaptation to new hosts in insects (Matzkin and Markow 2013). Here, aiming to test whether gene expression divergence is associated with adaptation to different cactus hosts, we performed differential expression analysis between species with Opuntia sp. preference versus columnar cacti preference, in adult heads and larvae. Because our fly stocks were reared in the laboratory with corn media, the observed transcriptome response is constitutive rather than triggered by the native cacti. Overall, the principal component analysis from HLAU gene expression (normalized read counts) revealed remarkable differences between species in both head and larval tissues (Fig. 3A). The variance of gene expression resembles the phylogenetic divergence of the species, with the highest variance between the buzzatii and mojavensis clusters (PC1) followed by divergence between D. buzzatii and D. koepferae (PC2). The species of the mojavensis cluster were clustered together with a slight separation between D. arizonae and the subspecies of D. mojavensis.

Figure 3.

Differential expression analysis of HLAU genes. (A) Principal component analysis of HLAU genes’ expression in head and larval tissues. Both tissues had HLAU genes’ variance representing the phylogenetic distance between species: D. buzzatii and D. koepferae, compared with the mojavensis cluster corresponds to the PC1 (45%), and 26% between them in PC2. Species from the mojavensis cluster had strong similarity, with D. arizonae slightly separated from D. mojavensis. (B) HLAU genes’ expression in the Opuntia sp. dweller D. m. wrigleyi versus its sister species D. arizonae, which does not have preference between Opuntia sp. and columnar cacti. The left plot represents expression in the head; the right plot, in larval tissue. (C) HLAU expression in head: differential expression analysis of HLAU genes from heads in the Opuntia sp. dweller D. m. wrigleyi versus the three columnar dwellers D. m. mojavensis, and D. buzzatii versus D. koepferae. (D) HLAU expression in larvae: the same pairwise species as in C with Opuntia sp. dweller versus columnar dwellers, but from larvae tissue. (B–D) Filled dots represent HLAU genes with log₂ fold-change > |1| and FDR < 0.05; transparent dots represent nonsignificant differential expression.

Open in new tabDownload PowerPointLink to figure

The PCA variance in the heads was highly correlated with ABC, UDP, and GST expression for PC1, whereas EST had a high correlation with PC2 (Supplemental Fig. S6A). In larvae, we observed that ABC, UDP, EST, and CYP are correlated to PC1, whereas GR, IR, and OBP are positively correlated to PC2 (Supplemental Fig. S6B). These results suggest that genes associated with host usage contribute to the variance of expression patterns between species in heads, whereas variance in larvae is represented by families from location, acceptance, and host usage.

We performed pairwise analysis of HLAU gene expression between species with Opuntia sp. preference (D. m. wrigleyi and D. buzzatii) versus species with preference for columnar cacti and with D. arizonae as nonpreferential cacti. All differentially expressed HLAU genes, as well as the presence of TEs in their promoters, can be accessed in Supplemental Table S10. In heads, D. m. wrigleyi had 27 and 19 HLAU genes with higher expression compared with D. m. mojavensis and D. m. sonorensis, respectively (Fig. 3C). Four genes were found with higher expression in the heads of D. m. wrigleyi compared with the other two subspecies: UDP-Ugt4, Gst-theta3, Abc-a3, and Cyp6a21. In the subspecies with columnar cacti preference, we observed, on average, 26 HLAU genes with higher expression compared with D. m. wrigleyi. Among them, seven were recurrent in the two subspecies: Ir21a, EstB1, UDP2, UDP-Ugt5, Cyp9b2, UDP-Ugt5-1, and Cyp12a2. Because they were observed in the two subspecies with columnar cacti preference, we propose that they are candidates to be related to the adaptation to columnar cacti. Finally, the same analysis on larvae revealed that D. m. wrigleyi has three genes with recurrent higher expression compared with the other subspecies—Or82a, Obp59a, and Ir21a; whereas in the other two D. mojavensis subspecies, we observed 17 common upregulated genes compared with D. m. wrigleyi, which comprise seven CYP genes, three GSTs, four UDPs, one EST, and one OBP (Fig. 3D). Regarding D. buzzatii and D. koepferae, we observed 27 HLAU genes in heads and 54 in larvae with higher expression in D. buzzatii, and 63 in heads and 29 with higher expression in D. koepferae. We also compared the head expression of HLAU genes between D. m. wrigleyi and its sister species, D. arizonae, without cacti host preference: D. arizonae. In D. m. wrigleyi, we found 31 HLAU genes with higher expression compared with its counterpart D. arizonae. In larvae, D. m. wrigleyi had 43 HLAU genes with higher expression compared with D. arizonae (Fig. 3B). This result demonstrates that HLAU genes have divergent differential expression in both heads and larvae.

To further investigate whether differential expression of HLAU genes has a higher prevalence compared with non-HLAU genes, we compared the proportion of differentially expressed genes associated with localization, acceptance, and usage with the proportion of differentially expressed genes from the non-HLAU gene families SER and ZF. In heads, our analysis revealed that localization and usage have higher differential expression than non-HLAU genes between D. buzzatii and D. koepferae (Supplemental Fig. S7). In D. mojavensis subspecies, we observed higher differential expression for localization genes between D. m. wrigleyi and D. m. mojavensis. In larvae, we did not observe any significant difference between non-HLAU and HLAU genes in D. buzzatii versus D. koepferae, but the pairwise differential expression between D. m. wrigleyi and two subspecies and D. arizonae revealed a higher prevalence of gene families associated with usage (detoxification) compared with non-HLAU (Supplemental Fig. S7). In D. m. wrigleyi versus D. arizonae, we also observed higher differential expression of localization genes than in non-HLAU. Taken together, the quantitative results suggest that adaptations that have allowed the shift from Opuntia sp. to columnar cacti may arise from the evolution of gene expression.

The enrichment of TEs on HLAU genes is not associated with changes in gene expression

Because we observed enrichment of TEs on the regulatory regions of a few host shift–related families, mainly ORs, and we also found differentially expressed HLAU genes between species, we tested whether the presence of TE insertions is associated with differential expression. To do so, we selected expressed HLAU genes (see Methods) and performed a X² independence test, splitting them into two variables with two categories each: (1) expression level—differentially expressed, and nondifferentially—and (2) TE presence—with TEs at 2 kb upstream and without TE at 2 kb upstream. Considering all pairwise analysis (Fig. 3B–D), the X² results demonstrated that the differential expression of HLAU genes is not associated with the presence of TEs on the promoter region (Supplemental Table S11). Because we did not observe overall effects of TEs on the expression of HLAU genes together, we investigated the association of TEs with differential expression only for HLAU gene families that we found TE enrichment (Fig. 2A). None of the enriched gene families had a significant association between the presence of TEs and their expression level (Supplemental Table S11). Thus, a global enrichment of TEs in the regulatory region is not likely to be associated with the differential expression of HLAU genes.

Make or break: gene–TE chimeras contribute either predominantly or marginally to gene expression

Although highly deleterious, a small fraction of TEs generating chimeric transcripts can be positively selected, giving rise to domestication and exaptation events (Capy 2021). Aiming to identify the extent of the transcriptome variability generated by gene–TE chimeric transcripts, we analyzed the transcriptome of head and larval tissues in the six cactophilic Drosophila species. The method employed in this work considered paired-end reads spanning from TEs to exons. Although systematically found in all replicates, many chimeras are likely to be products of pervasive transcription, generating ectopic transcripts whose expression levels are insufficient to play a functional role. Therefore, we developed an automatic method to assemble chimeric transcript isoforms and compute their expression contribution relative to the total gene expression. This method is implemented in the pipeline ChimeraTE v2.0 (Oliveira et al. 2023) as an optional downstream analysis to the total chimeras detection.

Our results revealed that, on average, 715 and 574 genes produce chimeric transcripts across all species, corresponding to 7.51% and 5.04% of the expressed genes in heads and larvae, respectively. In both, TE-exonized transcripts represented the highest frequency, with 87.59% in heads and 89.09% in larvae (Fig. 4A), followed by TE-terminated transcripts with 10.63% and 9.45% and in TE-initiated transcripts with 1.79% and 1.51%. All species have a similar number of expressed genes (heads: x¯ = 9508, σ = 300; larvae: x¯ = 11,375, σ = 117), but we observed on average 11-fold fewer chimeric transcripts in D. buzzatii and D. koepferae compared with the D. arizonae and D. mojavensis subspecies. This result suggests that either species from the mojavensis cluster are more prone to generate chimeras or that D. buzzatii and D. koepferae have a strong purifying selection against gene–TE chimeric transcripts.

Figure 4.

Chimeric transcripts in cactophilic Drosophila species. (A) Proportion of each chimeric transcript category based on the position of TE insertions relative to the gene structure. (B) Histograms with contribution of chimeric transcript isoform relative to total gene expression (x-axis). Chimeric transcripts without assembled sequence or contribution <25% were removed (gray bars). (C) Number of common and unique chimeric transcripts from HLAU genes detected in heads and larvae. (D) The Eato/LOC6578120 (ABC transporter) gene has a conserved TE-exonized transcript with a LINE/R1 element located in its seventh intron. The depth in exons, including the seventh intron, was assessed with RNA-seq from head tissue.

Open in new tabDownload PowerPointLink to figure

In heads and larvae, ∼49% of the detected chimeric transcript isoforms were not assembled or had a contribution to gene expression <25% (Fig. 4B). Among these discarded chimeras, the vast majority were not detected by the transcriptome assembly approach (0% contribution) (Fig. 4B), which is likely to be associated with negligible expression. Considering chimeric transcripts with contribution ≥25%, >95% have contribution equal to 100% (Fig. 4B). This result indicates that almost all genes producing chimeras are dependent on the TE sequence to constitute their complete mRNA sequence. Because ∼88% of all chimeras detected with >25% contribution are from TE-exonized transcripts, for which 72.35% are from TEs embedded in exons, a high frequency of chimeras with total gene expression could be expected. Notably, this pattern is also observed in TE-initiated and TE-terminated transcripts (Supplemental Fig. S8), despite the TEs being located upstream of and downstream from the gene body. Our results reveal a clear bimodal trend: Either chimeric transcripts have low expression and are unlikely to be biologically relevant, or they account for the total gene expression, suggesting a functional role.

HLAU genes generate chimeric transcripts

To investigate cases in which TEs produce chimeric transcripts from HLAU genes, as well as their possible role in the evolution of these genes, we identified all chimeras derived from HLAU genes in head and larval transcriptomes. We found a total of 28 (head) and 24 (larvae) HLAU genes producing chimeric transcripts. Although D. buzzatii had 11-fold fewer chimeras than did species from the mojavensis cluster, we observed 29 gene–TE (a few genes with more than one chimera) chimeras derived from HLAU genes in heads and larvae, the highest number among all species. Its sister species D. koepferae had only four chimeras, which are the same two ABC genes in heads and larvae (Supplemental Table S12). In D. buzzatii, ∼93% of all chimeras are from GSTs, expressed in both larvae and heads. All of them represent cases of Helitron exonizations (Supplemental Table S12). We did not find any GST with chimeric transcripts in the other transcriptomes, revealing an unprecedented recent evolutionary novelty in GSTs’ evolution in D. buzzatii.

Chimeric transcripts from HLAU genes are represented by 14.81% TE-terminated transcripts and 85.19% TE-exonized transcripts. We did not observe any TE-initiated transcript from HLAU genes. Despite the compelling evidence of TEs enrichment in the promoter of OR genes from species of the mojavensis cluster, we did not observe chimeric transcripts from ORs at high frequency. Only Or45b (LOC6578440), Or13a (LOC26527676), and Or85a (LOC26528431) had chimeric transcripts in D. m. sonorensis, D. m. wrigleyi, and D. arizonae, respectively (Supplemental Table S12). This result reveals that TEs located in the promoter of ORs are not likely to be co-opted as alternative TSSs. In the species of the mojavensis cluster, we observed ABCs as the gene family with the highest frequency of chimeras (Supplemental Table S12), despite not having enrichment (Fig. 2A). We observed, on average, five ABC genes with exonization of multiple TE copies in this group of species. This result indicates that the genomic distribution of TEs is not necessarily associated with their transcriptional interaction with genes.

The quantitative analysis on chimeric transcripts revealed transcriptional polymorphism in HLAU genes owing to TEs. We investigated common chimeras present in more than one species, highlighting their retention over evolutionary time, as well as species-specific chimeras that may be associated with recent insertions. Our results revealed that 94.44% and 88.57% of the chimeric transcripts derived from HLAU genes are species specific in heads and larvae, respectively (Fig. 4C). In heads, the Opuntia sp. dwellers D. buzzatii and D. m. wrigleyi were the ones with the highest number of HLAU genes with species-specific chimeric transcripts (Fig. 4C). The 15 chimeras uniquely found in D. buzzatii were derived from 12 genes (two ABCs and 10 GSTs), whereas the 14 chimeras from D. m. wrigleyi were observed from 10 genes (five ABCs, three IRs, one CYP, and one OBP) (Supplemental Table S13). In larvae, D. buzzatii follows a similar result as observed in heads, 15 chimeras derived from 14 genes (two ABCs and 12 GSTs), whereas D. m. wrigleyi had only four unique chimeras (two ORs, one IR, and one CYP).

Cases of common chimeric transcripts between populations/species may indicate exaptation/domestication events. In the six transcriptomes analyzed here, we found eight genes in the heads and four in larvae present in more than one genome (Supplemental Table S13). These common chimeric transcripts occurred only between D. mojavensis subspecies, revealing that TEs interacting with HLAU genes are likely to be lost over higher evolutionary scales. Their inclusion into the mRNA sequence and high contribution can be observed. For instance, the ABC transporter (Eato gene) has an exonized LINE/R1, which is located in the seventh intron (Fig. 4D), with an average contribution of 84.25% to the gene expression. Taken together, these results demonstrate unprecedented transcriptional interaction between TEs and HLAU genes in adult heads and larvae.

No positive selection in columnar cacti dwellers versus Opuntia sp. dwellers

The necrotic tissue of each cactus-host used by cactophilic Drosophila species differs in its composition in terms of yeasts and alkaloids (Fogleman et al. 1982; Fogleman and Danielson 2001). Therefore, species have adapted to feed and breed on different hosts, as demonstrated by behavioral and physiological response when species feed and breed in nonpreferential hosts. Positive selection in HLAU genes has been proposed to allow the recognition of a broader spectrum of odors, to metabolize nutrients, and to activate detoxification pathways (Markow 2019). Indeed, specialist species are likely to lose some of their genes related to sensory pathways, whereas the genes that are retained are subject to positive selection (McBride 2007; McBride and Arguello 2007). A preliminary genome-wide analysis of signatures of positive selection in the four D. mojavensis subspecies revealed high ω rates on genes associated with chemosensation, perception, immunity, behavior, detoxification, and reproduction (Allan and Matzkin 2019). Although a previous genome-wide analysis has identified 1294 genes under positive selection (Guillén et al. 2015), here we focused the analysis to specifically investigate HLAU gene families. Although the site model identified diversifying selection at specific sites of Atet and Cyp4ae1 genes, the branch and branch-site models revealed no evidence of diversifying selection. The branch and branch-site models analyze the common differences of species adapted to columnar cacti, rather than independent evolutionary events, making it more stringent. The lack of genes with statistical support (FDR < 0.05) for positive selection in the branch analyses may also be associated with technical limitations, because we based our analysis on one:one orthologs, and a certain proportion of the genes evolve too fast to be correctly assigned during ortholog identification (Kristensen et al. 2011).

OR promoters are enriched with TE insertions

The literature has recurrently demonstrated that TEs can create new regulatory elements, including gene promoters (Faulkner et al. 2009; Emera et al. 2012; Kapusta et al. 2013), enhancers (Xie et al. 2013, 2; Ye et al. 2020), and insulators (Schmidt et al. 2012; Wang et al. 2012). In mammals, TE-derived sequences contribute up to 40% of genome-wide binding sites for TFs located in TE-rich regions (Sundaram et al. 2014). It has been proposed that such co-option of TE-regulatory elements for gene regulatory networks can provide substantial modification of gene expression over short timescales, contributing to genetic variability at the transcriptome level (Mateo et al. 2014; van't Hof et al. 2016; Gill et al. 2021). The colonization of new habitats and hosts is accompanied by behavioral changes, mainly through alterations in environmental perception (Dekker et al. 2006; Matsuo et al. 2007). Here, we observed that ORs are enriched by TEs in the three subspecies from the mojavensis cluster. In ants, ORs have also been found enriched in TEs. In this case, TEs are associated with tandem duplication by crossing over, allowing massive duplication events (McKenzie and Kronauer 2018). Because there is no expansion of ORs in cactophilic species, the enrichment in TEs is likely to have another function, probably related to gene regulation owing to its specific location in promoters. Although we did not observe a significant enrichment of TEs in differentially expressed genes, we argue that TEs might contribute to the regulation of specific HLAU genes. In ORs, ∼39% of the genes have the OR's TFBSs located within TEs. We propose that this subset of genes is not sufficient to have a significant signal when tested with ORs without TEs, but they may have functional TE-derived cis-regulatory elements.

Helitrons have been extensively reported on the promoter of other eukaryotic species (Gupta et al. 2005; Cultrone et al. 2007; Thomas et al. 2014; Huang et al. 2016), in which a few cases of embedded TFBSs were reported (Carareto et al. 2014; Batista et al. 2019). Here, we demonstrated that ∼50% of the enrichment in TEs observed at the promoter region of ORs might be owing to the exaptation of TFBSs derived from Helitrons. In addition, the recurrent enrichment in species of the mojavensis cluster might demonstrate potential cases of TE co-option inherited by the common ancestor. Our genomic analysis provides a robust list of potential candidate genes in which Helitrons might play a regulatory role. This is an interesting hypothesis that can be determined through cutting-edge technologies such as CRISPR methods to delete insertions with site-specific resolution (Merenciano et al. 2023). Our results suggest a putative role of TEs in the regulation of HLAU genes across multiple species/subspecies.

Head transcriptome reveals subspecies-specific patterns

The differential expression of HLAU genes in heads revealed that localization and host usage (detoxification) gene families have the highest divergence across the subspecies, because we observed more differentially expressed genes for these families. In D. m. wrigleyi, which is the only D. mojavensis subspecies with a preference for Opuntia sp., we observed that CYP genes had on average 12 upregulated genes compared with the other subspecies. This result highlights that the divergence of cacti preference between D. mojavensis subspecies might be explained by modifications in gene expression of HLAU genes in the head. Although Opuntia sp. has lower toxicity than columnar cacti (Wright and William 2014), we found a similar number of upregulated detoxification genes in D. m. wrigleyi and compared to the other subspecies. This suggests that D. m. wrigleyi has a constitutive detoxification response that is independent of host toxicity. This is a plausible explanation because Opuntia sp. produces fewer benzaldehyde volatiles than the columnar cacti S. thurberi, but both cacti are lethal to D. melanogaster after 48 h (Wright and William 2014). Therefore, the adaptation to use Opuntia sp. must have a constitutive gene expression of the detoxification pathways in D. m. wrigleyi.

Comparative transcriptomic studies have revealed divergence in gene expression between D. buzzatii and D. koepferae, particularly in head tissues associated with sensory perception and detoxification. For instance, Guillén et al. (2015) reported significant differences in the expression of genes involved in olfaction and xenobiotic metabolism, suggesting that these species have adapted to their respective cactus hosts (Opuntia for D. buzzatii and columnar cacti for D. koepferae). Similarly, divergence in expression patterns has been reported in these species in ORs, OBPs, and CYP genes (Hasson et al. 2019). These differences are consistent with our results in heads and larvae, reinforcing the hypothesis that ecological specialization and divergence in head transcriptomes play a role in host shift between D. buzzatii and D. koepferae. Further studies must be conducted to explore the differential expression of HLAU genes in the gut, in which a significant portion of the detoxification occurs in Drosophila (Wu et al. 2015).

Previous studies with D. melanogaster have shown that sequence divergence in ORs and OBPs significantly modulates odor sensitivity (Matsuo et al. 2007; Hickner et al. 2016), potentially leading to divergence of host location. In addition, many ORs are involved in a broad range of ecological interactions, influencing behaviors associated with oviposition and feeding (Mansourian and Stensmyr 2015). For instance, OR evolution has mediated the herbivorous nature of leaf-mining specialist Scaptomyza flava (Goldman-Huertas et al. 2015). Furthermore, transcriptomic analyses have identified 18 ORs that are differentially expressed between heads of D. m. wrigleyi and D. m. mojavensis, demonstrating that modifications in the sensory transcriptome have evolved between these species (Crowley-Gall et al. 2016). In this study, we expanded the understanding of this divergence by comparing species with a preference for Opuntia sp. and columnar cacti, and showed that each species has a distinct set of differentially expressed HLAU genes.

Gene expression in larvae is one of the factors shaping the adaptation to new hosts

Differences in larval behavior and survival have been observed when Opuntia-related species are reared in columnar cacti (Ruiz and Heed 1988; Mateus et al. 2019). A previous study showed that D. m. mojavensis, which primarily uses S. thurberi as its host, has a 60% reduction in larval viability when reared on S. gummosus–based media (Matzkin 2012). The authors also proposed that 21% of the genes are differentially expressed between larvae reared in the primary and alternative cacti hosts. Many genes were reported with detoxification and energy production, among others. Another similar study with larvae transcriptome, but with D. buzzatii and D. koepferae, proposed that D. buzzatii flies reared in columnar cacti media have a stronger detoxification response in comparison with D. koepferae flies, suggesting that D. koepferae has a canalized transcriptome response to the toxic alkaloids from columnar cacti (De Panis et al. 2022).

Here, we observed that the majority of differentially expressed HLAU genes are from gene families associated with detoxification. Because D. m. wrigleyi and D. buzzatii use primarily Opuntia sp., we expected to find a higher number of detoxification genes in the other species/subspecies with columnar cacti preference. However, we observed a similar number of differentially expressed detoxification genes regardless of the host preference. CYPs were overrepresented in our results, but they represent a diverse gene family that can have functions other than detoxification, such as those related to developmental pathways in larvae (Tijet et al. 2001; Feyereisen 2012). Thus, further analysis of the biological role of upregulated CYP genes in D. m. wrigleyi and D. buzzatii might be assessed. Alternatively, as we observed in the head transcriptome, although Opuntia sp. has lower toxic compounds compared with columnar cacti, D. m. wrigleyi and D. buzzatii may also have a resistance to toxicity (Wright and William 2014). Thus, any cactophilic species is expected to have evolved a certain level of detoxification ability, which explains our observations. In this study, we propose a set of common HLAU genes in columnar cacti dweller species associated with detoxification. We observed that the three D. mojavensis subspecies have the same seven genes with higher expression compared with D. m. wrigleyi. Most of them are associated with detoxification (EstB1, UDP2, UDP-Ugt5, Cyp9b2, UDP-Ugt5-1, Cyp12a2), except Ir21a. It is noteworthy that other evolutionary forces than gene expression might be associated with host shift. For instance, gene duplication events have been reported in CYPs and GSTs in cactophilic Drosophila species (Rane et al. 2019). These duplications are thought to enable flies’ resistance to different toxic environments. Similarly, D. buzzatii and D. koepferae show both copy number variation and expression divergence in detoxification genes, indicating that gene duplication followed by regulatory or coding sequence divergence contributes to host-specific physiological adaptations (Guillén et al. 2015; Rane et al. 2019). These findings support the view that gene duplication serves as a key evolutionary mechanism facilitating ecological specialization in cactophilic Drosophila. Our results show that divergence in the regulation of gene expression in larvae is one of the factors likely to shape the adaptation to new hosts.

TEs are a source of transcriptome variability between cactophilic Drosophila species

TEs contribute to genome and transcriptome evolution by rewiring regulatory networks and by incorporating protein domains into gene transcripts. In both cases, TE copies that are inserted within genes can integrate their sequences into the mRNA, generating chimeric transcripts. A recent study with different ecotypes of D. melanogaster has shown that these chimeras have the potential to generate transcriptome novelties (Oliveira et al. 2023). These chimeras have also been demonstrated to be active in a tissue-specific manner (Coronado-Zamora and González 2023). Although the production of chimeric transcripts derived from polymorphic TE insertions is well known, there is no study addressing the question of how TEs contribute to the evolution of genes associated with host shift. Here, we took advantage of genomic and transcriptomic data to analyze the transcriptional interaction of TEs with HLAU genes through gene–TE chimeric transcripts. Our finding demonstrated that around one-half of the chimeras have low expression levels (0%–5%) and are likely to represent pervasive transcription of TEs. Notably, >95% of the chimeras overpassing the 5% contribution to gene expression constitute nearly the total gene expression. In this analysis, we included only TE copies located within 3 kb upstream and 3 kb downstream regions, although longer-range gene–TE interactions have been reported (Chuong et al. 2016). Further studies integrating RNA-seq with 3D chromatin folding approaches such as Hi-C could help identify chimeric transcripts from TEs in trans.

Considering previous findings in different tissues of D. melanogaster, in which 38% of all chimeric transcripts express only the chimeric isoform (Coronado-Zamora and González 2023), our result shows a higher prevalence of chimeras in the repertoire of gene splicing. In tetrapods, the exonization of TE-derived transposases has been shown to constitute the unique gene isoform (Cosby et al. 2021). In HLAU genes, all chimeric transcripts represent TE-exonized transcripts. We suggest that TEs in the promoters of ORs might still recruit TFs at their terminal sequences, which would prevent their inclusion in the gene mRNA (chimera). In addition, the maintenance of TEs in the promoters of ORs might also be owing to epigenetic marks. TEs can have chromatin marks (Lippman et al. 2004; Pal et al. 2023), making them versatile regulatory motifs acting in the regulation of nearby genes (Guio et al. 2018; Coronado-Zamora and González 2023). Further experimental analysis, as ChIP-seq, may be performed to test the TEs’ function in OR's transcription. Furthermore, TE annotation remains challenging in nonmodel organisms. Although we performed several steps to remove potentially false TE copies from the genome, we can not rule out the possibility of remaining annotation errors in the exons of HLAU genes.

A comparative analysis between generalist and specialist Drosophila species revealed that niche amplitude is not likely to play a role in TE dynamics in the genomes (Fonseca et al. 2019). Nevertheless, several studies have proposed that the successful adaptation to new environments in invasive species might be associated with TE activity, despite reduced genetic diversity caused by founder effects (Marin et al. 2021; Mérel et al. 2021). Here, we demonstrate that TEs interact with genes associated with host shift, providing a transcriptomic variability in six cactophilic species with different host preferences. Our results show that ABC transporters produce chimeric transcripts in all cactophilic species, highlighting a compelling mechanism of structural innovation. In insects, MITEs (nonautonomous DNA transposons) have been repeatedly identified within defensome genes, including ABC transporters, such as in Helicoverpa spp., suggesting a potential causal link to insecticide resistance (Zidi et al. 2022). Functionally, exonized TEs are stably expressed and translated, often producing altered protein domains or localization signals, thereby expanding proteomic repertoires while providing raw material for adaptive selection (Arribas et al. 2024). Taken together, our findings suggest that TE exonization in ABC transporters may offer cactophilic Drosophila a versatile genomic strategy to fine-tune transporter function, potentially under the selective pressures of xenobiotics. Although we found a high contribution of the chimeric isoform to total gene expression, future work should be conducted to elucidate which exonized TEs are translated and whether they modulate substrate specificity, protein stability, or cellular localization, thereby affecting organismal fitness.

Our discovery of 13 GST genes in D. buzzatii harboring exonized Helitrons underscores an unprecedented finding of TE-driven functional diversification in detoxification enzymes. Helitrons have been shown to supply novel exons, splice sites, and promoters across eukaryotes (Thomas et al. 2014). This finding, unique to D. buzzatii, suggests a lineage-specific co-option of Helitrons to expand and diversify GST isoforms. GSTs are central to xenobiotic metabolism, and exonization of Helitron-derived sequences could introduce novel domains or alter enzymatic properties, potentially conferring adaptive advantages under environmental pressures such as exposure to cactus-derived toxins. Future work should be performed to determine the function of these Helitron-derived motifs in the host preference of D. buzzatii.

It is important to state that the identification of chimeric transcripts through transcriptome analysis is an important, but limited step in terms of function. These mRNA molecules might be degraded by surveillance pathways that control for aberrant mRNA production, such as no-go decay (Harigaya and Parker 2010), nonstop decay (Vasudevan et al. 2002), and nonsense-mediated RNA decay (Hug et al. 2016). Our study provides a set of HLAU genes producing chimeric transcripts, in which part of them might have a phenotypic impact and ultimately be associated with host preference. However, it is fundamental to perform comparative phenotypic studies of host adaptation cues with mutant strains for the TEs interacting with HLAU genes identified in this study. The application of techniques of RNAi aiming at chimeric transcripts, as well as the deletion of TE insertions related to chimeras with CRISPR-Cas9 (Merenciano et al. 2023), would be crucial to confirm the predictions of our results.

Potential role of TEs in host shift

Host specificity is a key mechanism of reproductive isolation in plant pathogens and is regulated by the repertoire of effector genes within each pathogen. In Phytophthora sp., the host specificity is partly regulated by RXLR class effectors that facilitate host exploitation. Notably, synthetic chimeras of a short interspersed element (SINE) linked to an effector gene in Phytophthora infestans induced their silencing (Whisson et al. 2012). This silencing likely occurs naturally, as transcriptional inactivation of effectors has been reported, and more than one-half of the RXLR effectors in the P. infestans genome are located in TE-rich regions (Raffaele et al. 2010). Consequently, TEs inserted near these genes may have influenced host shift in P. infestans. Our findings demonstrate multiple interactions between genes and TEs, highlighting the contribution of TEs to the genome and transcriptome evolution of cactophilic species. Although our data did not allow the quantification of the TE insertions frequency in natural populations, our results demonstrate the first transcriptome-wide evidence of TE co-option in cactophilic Drosophila species. Genes producing chimeras with TEs require further investigation, as gene–phenotype associations involving HLAU genes have not been fully characterized, preventing us from establishing a specific link between chimeras and potential host preference. In addition, our study focused exclusively on nine gene families that have been shown to evolve during host shift events (Markow 2019). However, host shift is undoubtedly a complex trait to study, and many other genes can be associated with the adaptation to new hosts. For instance, differential expression of genes associated with development and neurological processes has been identified in cactophilic Drosophila species reared in different cacti media (De Panis et al. 2016). Behavioral changes are also based on host shift evolution and have been reported in many systems (Knolhoff and Heckel 2014). Despite focusing on nine relevant gene families, we suggest further studies including other aspects of host shift than localization, acceptance, and host usage.

The existence of reproductive barriers between populations that undergone a host shift is the key process to ecological speciation. In our model, with species from the buzzatii and mojavensis clusters, incipient reproductive isolation has been reported between D. m. sonorensis and D. m. baja (Zouros and d'Entremont 1980; Markow 1991; Pfeiler et al. 2009). The premating isolation between them occurs owing to the difference in cuticular hydrocarbons, which participate in the signaling pathway of mate recognition (Howard and Blomquist 2005). Importantly, artificially shifting the cactus host was the main cause of changes in the hydrocarbon profiles. Thus, the source of reproductive barriers is dependent on the cacti species on which the flies feed and breed (Stennett and Etges 1997; Havens and Etges 2013). The close relationship between HLAU genes to the host shift process, and the observed contribution of TEs to their evolution, supports a possible role of TEs in the host preference of cactophilic species. However, it is important to highlight that HLAU gene families are not likely to represent the only genes involved in host shift. Other gene families can also contribute to this process, such as genes associated with development (Savković et al. 2023) or with decision making (Williams-Simon et al. 2019). In addition, we treated all genes from each gene family as having functions directly associated with location, acceptance, or detoxification. Nevertheless, some genes may have different functions than those, such as CYP genes in hormone metabolism (Sutherland et al. 1998). Because we do not have experimental data showing these functions in D. mojavensis, we maintained all genes. Further studies must be conducted to test the function of the genes reported in our work as HLAU-related.

Fly stocks and genome sequencing

Most of the strains were obtained from the UC San Diego Drosophila Stock Center (D. m. mojavensis: Anza [01], Anza Borrego Desert, California, Stock Center 15,081–1352.01; D. m. wrigleyi: CI [22], Catalina Island, California, Stock Center 15,081–1352.22; and D. m. sonorensis: AG [26], Agiabampo Bay, Sonora, Mexico, Stock Center 15,081–1352.26) and D. arizonae (HI [17], collected in Metztitlan, Hidalgo, Mexico, Stock Center 15,081–1271.17). The D. buzzatii (Bu28) and D. koepferae (Ko2) stocks correspond to two laboratory lines used in previous works (Romero-Soriano et al. 2016; Bodelón et al. 2022). Both stocks were maintained by brother–sister mating for more than a decade and then kept by mass culturing.

We used the same protocols of DNA extraction and genome sequencing for all samples. DNA was extracted from 10 males and 10 females from each species using the Qiagen DNeasy blood & tissue kit. Then, we evaluated the genomic DNA amount and quality with NanoDrop One UV-vis spectrophotometer (Thermo Fisher Scientific) and Qubit 1.0 fluorometer (Invitrogen). Three micrograms of DNA was repaired using the NEBNext FFPE DNA repair mix (NEB M6630). We performed end repair and dA-tailing using the NEBNext end repair/dA-tailing module (NEB E7546). Ligation was then assessed with the ligation sequencing kit 1D. MinION sequencing was performed according to the manufacturer's guidelines using R9.4.1 flow cells (ONT FLO-MIN106) and a Nanopore MinION Mk1b sequencer (ONT) controlled by ONT MinKNOW v.18.3.1. Base calling was performed after sequencing using Guppy v.4.0.5 in high accuracy mode for isogenic wild-type strains v.3.1.5.

Genome assemblies and gene annotation

The quality control for Nanopore reads was performed with NanoPlot v.1.10.2 (https://github.com/wdecoster/NanoPlot). Reads with quality lower than seven were removed from downstream analysis. To assemble contigs, Flye v.2.8 (Kolmogorov et al. 2019) was used with default parameters, except –plasmids. All contigs were aligned with minimap2 v2.16 (Li 2018), with the –x map-ont option. The alignment was used to perform the polishing of contigs with four rounds of Racon v1.3.2 (Vaser et al. 2017), with default parameters. The assembly quality metrics were assessed with Assembly Stats v1.0.1 (https://github.com/sanger-pathogens/assembly-stats). Assembly incongruences were manually visualized with D-GENIES v1.2.0 (Cabanettes and Klopp 2018) and corrected with SAMtools v1.9.0. (Li et al. 2009), function faidx, as well as Gepard v1.4.0 (Krumsiek et al. 2007) for the determination of breaking points. The superscaffolding of all the corrected assemblies was performed with RaGOO v.1.1 (Alonge et al. 2019), with –s and –t 4 parameters, using the respective reference genome assembly for each species. Subsequently, the chromosome-scale assemblies were submitted to a BUSCO v.5.4.3 analysis (Simão et al. 2015) with lineage diptera_odb10 (-l parameter). The gene annotations for D. mojavensis subspecies and D. arizonae were performed with the software Liftoff v.1.6.3 (Shumate and Salzberg 2021), with the reference genome of D. mojavensis (GCA_018153725.1) (Kim et al. 2021). We have used the parameters –exclude_partial; -a 0.5 and –s 0.75 to keep annotated genes with at least 75% of the reference gene length. We removed from our gene annotation HLAU genes flagged as “low quality protein” and “unestablished gene function” in the reference D. mojavensis genome.

D. buzzatii and D. koepferae genomes were annotated with a de novo strategy with BRAKER3 v.3.0.8 (Gabriel et al. 2024). To optimize the prediction of gene structures and repertoire of isoforms, we used a mix of RNA-seq data from different tissues publicly available and produced in this work (Supplemental Table S14), in addition to protein sequences from D. mojavensis as reference. Ortholog genes between D. buzzatii, D. koepferae, D. mojavensis, and D. arizonae were obtained with OrthoFinder v.2.5.4 (Emms and Kelly 2019) using DIAMOND v.0.9.14 (Buchfink et al. 2015) as search engine. One-to-one orthologs were obtained directly from OrthoFinder, and remaining proteins without an ortholog pair were analyzed with BLASTP v.2.5.0 (Altschul 1997), for which we considered as one-to-one proteins only those with a single homologous with the parameter -perc_identity 90 and with at least 95% of the reference protein length. Because OrthoFinder performs its search using protein sequences, the identification of orthologs from noncoding RNA genes was assessed with Liftoff (Shumate and Salzberg 2021) using all the noncoding RNA genes from D. mojavensis as a reference to maximize our homology-based gene annotation.

Molecular evolution analysis

To identify HLAU genes with signatures of positive selection in species that prefer columnar cacti compared with those that prefer Opuntia sp., we selected the three D. mojavensis subspecies, D. navojoa, D. koepferae, and D. buzzatii. We removed D. arizonae from this analysis owing to its generalist preference for both columnar and Opuntia cacti. The D. navojoa genome was retrieved from NCBI (GCF_001654015.2) (Vanderlinde et al. 2019). For the other Nanopore genomes sequenced in this work, we performed a polishing step with Illumina paired-end reads. We used the tool NextPolish (Hu et al. 2020) to fix base errors in the assembled genomes. The three D. mojavensis subspecies were polished with genomic Illumina reads sequenced previously (Banho et al. 2021), whereas for D. koepferae and D. buzzatii, we used publicly available genomic Illumina reads (Moreyra et al. 2023) and RNA-seq reads (Guillén et al. 2015), respectively. For RNA-seq data, to increase the read coverage and optimize the polishing to the highest number of genes and exons as possible, we merged RNA-seq reads from adult males and females (Guillén et al. 2015). All genomes were corrected using two rounds of polishing, as recommended (Hu et al. 2020). Then, we extracted the longest CDS sequences for all genes, combining AGAT (https://github.com/NBISweden/AGAT) agat_extract_sequences.pl parameters -t mrna, SAMtools faidx (Li et al. 2009), and seqtk (https://github.com/lh3/seqtk). Posteriorly, we intersected the occurrence of orthologous genes between the six genomes, recovering a list of genes that were present in all of them. The sequences were aligned using MAFFT v.7.487 (Katoh and Standley 2013), and the non-aligned codons were removed using PAL2NAL (Suyama et al. 2006). To obtain a strongly supported tree for this analysis, we inferred the species phylogeny tree with 2043 shared BUSCO genes. The phylogeny was estimated using the concatenated alignment (3,730,914 sites) with IQ-TREE v. 2.2.2.6 (Minh et al. 2020), using the best model for each gene, 10,000 replicates of ultrafast bootstrap, and 1000 replicates for an approximate likelihood-ratio test (ALRT) (parameters -m MFP -bb 10000 -alrt 1000). Finally, we carried out an analysis to detect signatures of positive selection for each gene individually with paPAML v.2.13 (Steffen et al. 2022), which uses CODEML from PAML v.4.10.7 (Yang 2007). We used the models branch-site, site, and branch. To find signatures of positive selection in species that use columnar cacti with branch and branch-site, we marked the respective branches of D. koepferae, D. m. mojavensis, and D. m. sonorensis with the same label (#1) to compare them with the background species that use Opuntia sp. (D. buzzatii, D. m. wrigleyi, and D. navojoa). Thus, we tested two hypotheses: (1) the null hypothesis (H1), in which all branches evolved without selection pressure (model = 2, Nssites = 2, fix_omega = 1, omega = 1) and (2) the alternative hypothesis (H2), in which selection pressure can vary between labeled branches (foreground) and background (model = 2, NSsites = 2, fix_omega = 0, omega = 1.5). The hypotheses were tested through a χ² test, calculated by the log-likelihood (lnL) ratio between H1 and H0 (2ΔlnL) and degrees of freedom obtained from the difference between the number of estimated parameters for H1 and H0 (Δnp). Sites under positive selection were identified by a Bayes empirical Bayes and posterior probability, implemented in CODEML. The P-value was adjusted by Benjamini–Hochberg procedure to control the false-discovery rate (FDR) (Benjamini and Hochberg 1995). Only genes with FDR < 0.05 were considered significant.

RNA extraction and sequencing

We performed RNA-seq for larvae and adult heads from D. buzzatii, D. koepferae, D. arizonae, and the three D. mojavensis subspecies. The RNA extraction was assessed with the RNeasy plus kit for 10 female heads (10 days old) and 15 larvae, with three technical replicates. Female flies were selected because, unlike males, they seek out the cacti to oviposit. This additional behavior might require fine-tuning adaptations in the environment perception that could potentially be absent in males. The library preparation was done with stranded mRNA-seq/standard quantity (ligation), and 100 bp paired-end reads were produced in an Illumina HiSeq 4000. Illumina adaptors and low-quality reads were removed for downstream analysis with Trimommatic v.0.39 (Bolger et al. 2014), witt the following parameters: LEADING:3; TRAILING:3; SLIDINGWINDOW:4:15; HEADCROP:10; CROP:90.

TE annotation

We developed a novel pipeline combining automatic and manual approaches to build TE consensuses for each sequenced genome. First, we used Earl Grey v.2.2 (Baril et al. 2024) to build de novo TE libraries with RepeatModeler2 v.2.0.4 (Flynn et al. 2020), which uses RECON v.1.08 (Bao and Eddy 2002), and RepeatScout v.1.0.6 (Price et al. 2005). We used the parameters -r 32281 (NCBI taxon ID for Drosophila), -c yes (cluster redundant consensus), -m yes (remove consensus with <100 nt), and -e yes (enable HELIANO v1.2.1 to identify Helitrons) (Li et al. 2024). In short, Earl Grey automatically retrieves TE consensus and performs an extension of each sequence using a method named “BLAST, extract, extend” (BEE) (Platt et al. 2016). This method extracts the sequences of the 20 full-length (or longest) TE insertions for each consensus, and then it retrieves the 1 kb flanking regions from them. Posteriorly, insertions altogether with flanking regions are aligned with MAFFT v.7.487 (Katoh and Standley 2013), and the alignment is trimmed with trimAl v.1.4 (Capella-Gutiérrez et al. 2009), with the parameters –gt 0.6 and –cons 60. Subsequently, the alignment is analyzed to check whether the flanking regions from the full-length insertions have high quality, meaning that they are part of the TE rather than genomic regions. Finally, the consensus sequence for each TE family is updated with EMBOSS v.6.6.0.0 (Rice et al. 2000). The BEE method is repeated up to 10 times for each consensus.

Once the extended TE libraries were obtained for the six genomes with Earl Grey, we implemented our pipeline into eight steps. First, we removed all sequences other than TEs from Earl Grey's library: satellites, low complexity repeats, tRNA, and rRNA. At step 2, we checked the presence of redundant consensus respecting the 80–80 (Wicker et al. 2007) rule with all-versus-all BLASTN (-qcov_hsp_perc 80 and -perc_identity 80). We found that ∼45% of the consensuses were duplicated. Some of them were identical, or nearly identical, but in opposite strands. We selected the longest consensus from the ones considered redundant. At step 3, we used Tandem Repeats Finder (Benson 1999) to identify consensuses that had >50% annotated as tandem repeats. Consensus falling in this cutoff were removed as they were indistinguishable from tandem repeats. At step 4, we removed all consensus with high similarity with D. m. wrigleyi CDSs using BLASTN, with the parameters -qcov_hsp_perc 80 and -perc_identity 80. TE consensuses that had passed our filtering thresholds were then submitted to the classification analysis. At step 5 of the pipeline, TE consensus was classified at the family level based on similarity analysis of nonreference consensus with a database provided by the user. In our case, we used all consensuses from Drosophila species available on Dfam v3.7 (Storer et al. 2021). Because the “curated” Dfam library contains only D. melanogaster consensus, we used the “uncurated” library from Dfam. We extracted Drosophila TEs with famdb.py v.0.4.3, with the parameters -i Dfam.h5 families ‐‐include-class-in-name -f fasta_name -ad “Diptera”. Then, we filtered out non-Drosophila sequences and “unknown” TEs with the Unix command line grep -i “drosophila” file.fa | grep -v -i “unknown” and obtained 49,916 out of 344,562 consensus. The pipeline firstly classifies at the family level by doing a BLASTN (-qcov_hsp_perc 80 and -perc_identity 80) between the Drosophila Dfam consensus and the provided de novo consensus. Approximately 50% of the identified consensus were classified on the family level based on this BLASTN. The consensus without matches from this stringent BLASTN and <200 nt were filtered out. Because of the stringency of the matches with BLASTN, we also masked the remaining unclassified consensus with the Drosophila Dfam consensus with RepeatMasker v4.1.8 (Smit 2004). Based on RepeatMasker output, we computed the length of all matches from a given family in the unclassified consensus and how much it has been covered based on the unclassified length. If the coverage is >80% of the total length, we assume that the unclassified consensus has diverged from that of the database, but it is still classified as the same matched family. This step classified ∼25% of the consensus. Finally, the remaining unclassified sequences that could not be classified at the family level are classified at the class/order level. We used a cutoff of 50% of covered matches in the unknown consensus to classify their classes. As the family name, we substituted the pattern from RepeatModeler2 as (i.e., rnd-1_family-1) to a pattern representing “species prefix,” “author initials,” and “homology percentage” (i.e., Dari_DSO_86.0275p#DNA/TcMar-Tc1). In this example, “Dari” stands for D. arizonae, “DSO” is the author's initials, and “86.0275p” represents the percentage of the total length that this consensus has been matched with a DNA/TcMar-Tc1. The remaining putative TE consensuses with <50% of coverage with known Drosophila TEs were removed. At the end of this step, we obtained the polished libraries of TE consensus.

In step 6, we checked the strandness of TE consensus sequences based on stranded RNA-seq data. During the clustering of similar sequences and removal of redundancy, we noticed many consensuses with nearly identical sequences but in different senses. To correct TE strandness, we used our stranded RNA-seq data. We annotated TE insertions in the genome with RepeatMasker v.4.1.4 (122), with the parameters -cutoff 250, -s, and -norna. Then, we mapped the reads from a mixed library from head and larvae with STAR v2.7.11b (Dobin et al. 2013). To be certain about the strandness origin of the reads, we removed all insertions that have any overlap with other TEs and/or exons. Then, for each consensus, the pipeline identifies the three copies in the genome with higher expression. Based on these copies, the number of uniquely mapped reads to the same strand of the copy is quantified. If >50% of the reads align to the same strand of the copies, we assume the consensus strandness was right. Otherwise, the pipeline converts the consensus sequence to its reverse complement. This step will add a strandness tag to the consensus ID. For instance, if 75% of the reads support a forward strand, a tag “75st” will be added to the family name. This is useful to provide a metric of certainty about the strandness of each consensus. Such information can be valuable when working with piRNAs, in which alignment strandness is crucial to interpret the results. Then, we used RepeatCraft (Wong and Simakov 2019) with LTR_Finder v.1.2 (Xu and Wang 2007) results to merge LTRs and their respective internal regions. Subsequently, we removed all TE insertions <80 nucleotides, following the classification rule 80–80–80 (Wicker et al. 2007). Finally, we performed a last filtering step in the TE annotation to remove insertions composed of simple sequence repeats. We used Tandem Repeats Finder v.4.09 (Benson 1999) in all TE insertions from each genome, with the parameters 2, 5, 6, 75, 20, 50, 500, -m, -h, to remove TEs with masked repeats over >50% of their sequences. This method has been previously used to remove TE copies that are indistinguishable from simple repeats (Flutre et al. 2012; Oliveira et al. 2023).

TE enrichment at regulatory regions

To identify TE insertions within and nearby genes for each genome, we created BED files with the genomic positions of all insertions. Then, we created BED files for the gene regions corresponding to the 2 kb upstream of the TSS. To test whether the prevalence of TEs in each gene region was significantly higher than expected (enriched) compared with random genomic distributions of TEs in the same region of other genes, we performed permutation tests with the R package regioneR (Gel et al. 2016), with the parameter randomization = resampleRegions, with 2000 random samples. For each gene family with n genes, the number of genes with TEs was compared with other random n genes 2000 times. The n depends on the size of the HLAU gene family used in the input. For instance, the permutation test identifies how many ORs have TEs upstream out of 59 genes in total (family size). Then, another 59 genes are picked randomly, and the TE frequency is identified. This process creates a normal distribution with 2000 observations, allowing us to estimate if the number of OR genes with TEs is higher or lower than expected considering the genomic background. We also selected ZF proteins and SER as two non-HLAU gene families to compare our findings. Taking into account that TEs might harbor regulatory and/or coding motifs in both strands (Mateo et al. 2014), our analysis considered TEs within the investigated gene regions regardless of the strand.

Identification of TE-rich genomic regions

We split the genomes into bins of 200 kb, generating a BED file with such intervals across the genomes. Then, we used the function intersect from BEDtools v.2.30.0 (Quinlan and Hall 2010) to identify the number of TE insertions located in each genomic bin. Posteriorly, we merged all overlapped TE insertions into single ranges. These ranges were summed to obtain the total TE content in base pairs for each 200 kb bin. By plotting the TE content of each genome, we noticed that bins with >50 kb of TE content were outliers of the normal distribution. Therefore, we defined TE-rich regions as all bins with >50 kb of TE content, which represents 25% of the bin size. We used the same function from BEDtools to identify whether OR genes were present in the TE-rich regions or the TE-poor (bins with <50 kb of TE content).

TF binding sites on TE insertions

We assessed the identification of TFBSs with FIMO (Grant et al. 2011) from the MEME v. 5.0.5 suite tools (Bailey et al. 2015). We extracted the sequences of all TE insertions located 2 kb upstream of ORs with BEDtools v. 2.30.0, function getfasta, using the parameter -s (Quinlan and Hall 2010). We obtained the TFBS motifs described for D. melanogaster from JASPAR database (Fornes et al. 2019) for the TFs acj6, zf30c, fer1, onecut, and xbp1, which have been demonstrated to activate ORs (Jafari et al. 2012). To compare our findings with OR genes, we performed the same analysis with SER and ZF as non-HLAU gene families. Only TFBSs with P-value < 0.01 were considered. The enrichment of OR TFBSs in OR genes compared with non-HLAU gene families was assessed with the four frequencies: (1) ORs with TFBSs, (2) ORs without TFBSs, (3) non-HLAU with TFBSs, and (4) non-HLAU without TFBSs.

Differential expression of genes

We aligned the RNA-seq reads with STAR v.2.7.3 (Dobin et al. 2013), using the assembled genomes and annotations of their respective species. The alignments were used to create the count matrix reporting the expression of each gene, with HTSeq v.0.13.5 (Anders et al. 2015), with the parameters –order pos; -s reverse; –i gene; -t exon; the latter indicates that the count for each gene corresponds to the sum of all reads aligned into exons (CDSs and UTRs). We merged the counts from species of the buzzatii and mojavensis clusters separately, maintaining only ortholog genes that have at least 90% of the length of the annotated gene in the reference genome. This threshold allowed us to avoid misinterpretation of the expression level by comparing genes with different lengths. Then, we performed independent statistical analysis for species of the buzzatii and mojavensis clusters with DESeq2 v.1.28.1 (Love et al. 2014), using a single-factor experimental design, which considers only the factor that different species may have differentially expressed genes. Because of the phylogenetic distance between species of the buzzatii and mojavensis clusters, we performed the differential expression analysis separately for each group. We only considered genes with an adjusted P-value < 0.05 and |log2 fold-change| > 1 as differentially expressed. To test for association between TE enrichment and differential expression, we selected HLAU genes with the mean of normalized counts between replicates ≥100. Then, we built the contingency table to apply the X² independence test, splitting them into four categories: differentially expressed, nondifferentially expressed, with TEs at the promoter region, and without TEs at the promoter region. We carried out the test for all pairwise analyses between cactophilic species with preference for Opuntia sp and columnar cacti. A X² with P-value < 0.05 was considered significant. The proportional differences between localization, acceptance, and host use genes compared with non-HLAU genes (SER and ZF gene families) were assessed with a Fisher's exact test, with a contingency table containing the average number of genes per group: (1) non-HLAU differentially expressed genes, (2) non-HLAU nondifferentially expressed genes, (3) HLAU differentially expressed genes, and (4) HLAU nondifferentially expressed genes. A Fisher's exact test < 0.05 was considered significant.

Identification of gene–TE chimeric transcripts

We detected gene–TE chimeric transcripts with ChimeraTE v.2.0.0 mode 1 (Oliveira et al. 2023). Each species was analyzed using its respective genome, gene annotation, and TE annotation. Shortly, ChimeraTE identifies chimeras based on paired-end reads spanning from exons to TE sequences. Based on the TE position compared with the gene (upstream, exon, intron, downstream), chimeric transcripts are classified into categories. We have used default parameters, except –replicate 3, which considers true chimeras only those identified in the three RNA-seq replicates, and –coverage 10, used to obtain chimeras with at least 10 chimeric reads on average between replicates. To calculate the contribution of the chimeric transcript expression relative to total gene expression, we developed a new module for ChimeraTE. This module can be used downstream to ChimeraTE mode 1 analysis. In the first step, all RNA-seq replicates will be merged into a single FASTQ library. Then, the merged reads are aligned to the genome with STAR (Dobin et al. 2013). Next, StringTie2 (Kovaka et al. 2019) is used to assess transcript expression, including both reference and nonreference isoforms. All TEs generating chimeras are intersected with exons predicted by StringTie2 with PyRanges (Stovner and Sætrom 2020). Only overlaps with 80% of the TE within the exon are considered. Finally, the expression of the TE-containing isoforms is compared with the total gene expression. In our analysis, we removed all chimeric transcripts for which StringTie was not able to build the chimeric isoform. In addition, aiming to keep chimeras with putative function, we selected only the ones with a contribution ≥25%. This analysis can be reproduced with any output from ChimeraTE with the code “mode1_contrib_chimeras.py” (https://github.com/OliveiraDS-hub/ChimeraTE/blob/main/scripts/mode1_contrib_chimeras.py).