Effects of Arboviral Infections on Transposable Element Transcript Levels in Aedes aegypti

Abstract

Transposable elements are mobile repeated sequences found in all genomes. Transposable elements are controlled by RNA interference pathways in most organisms, and this control involves the PIWI-interacting RNA pathway and the small interfering RNA pathway, which is also known to be the first line of antiviral defense in invertebrates. Using Drosophila, we recently showed that viral infections result in the modulation of transposable element transcript levels through modulation of the small RNA repertoire. The Aedes aegypti mosquito is of particular interest because almost half of its genome is made of transposable elements, and it is described as a major vector of viruses (such as the dengue [DENV], Zika [ZIKV], and chikungunya [CHIKV] arboviruses). Moreover, Aedes mosquitoes are unique among insects in that the PIWI-interacting RNA pathway is also involved in the somatic antiviral response, in addition to the transposable element control and PIWI-interacting RNA pathway genes expanded in the mosquito genome. For these reasons, we studied the impacts of viral infections on transposable element transcript levels in A. aegypti samples. We retrieved public datasets corresponding to RNA-seq data obtained from viral infections by DENV, ZIKV, and CHIKV in various tissues. We found that transposable element transcripts are moderately modulated following viral infection and that the direction of the modulation varies greatly across tissues and viruses. These results highlight the need for an in-depth investigation of the tightly intertwined interactions between transposable elements and viruses.

mosquito, DENV, CHIKV, ZIKV, retrotransposon, TE control

Significance

Transposable elements are genomic parasites that are found in all genomes. In insects, similar RNA interference pathways control transposable element activity and fight against viral infections. In this study, we focused on A. aegypti, the yellow fever mosquito, which harbors almost 50% transposable elements in its genome and is known to be a major vector of arboviruses. Using various publicly available RNA-seq datasets, we showed that arboviral infections induce moderate modulations of transposable element transcript levels in mosquito tissues. In addition, we revealed that the direction of the modulation strongly depends on the virus. These results highlight the need for an in-depth investigation of the tightly intertwined interactions between transposable elements and viruses.

Introduction

Transposable elements (TEs) are repeated DNA sequences found in all genomes, in which their amounts can vary due to their capacity to move within chromosomes. For example, TEs constitute approximately 3%, 47%, and 78% of the genomes of the honey bee Apis mellifera (Honeybee Genome Sequencing Consortium 2006; Elsik et al. 2014), the mosquito Aedes aegypti (Nene et al. 2007), and the arctic krill Euphausia superba (Shao et al. 2023), respectively. The ability of TEs to transpose within host chromosomes and their high abundance in many organisms have made TEs drivers of evolution (Biémont 2010). TE insertions are generally neutral or slightly deleterious. They may induce double-strand DNA breaks, ectopic recombination or disruption of coding sequences (Goodier and Kazazian 2008; Zamudio et al. 2015), but they are also involved in adaptation and genetic innovations (Casacuberta and González 2013; Gilbert et al. 2021; Galbraith and Hayward 2023; Lawson et al. 2023). TEs can be classified into two principal classes according to their transposition intermediates. Class I TEs or retrotransposons have an RNA intermediate and transpose via a “copy-and-paste” mechanism. Class II TEs are TEs with DNA intermediates called DNA transposons and are mainly mobilized through a “cut-and-paste” mechanism (Wells and Feschotte 2020). These classes are further subdivided into subclasses, families, and subfamilies according to their replication mechanism, structure, and sequence similarity (Wicker et al. 2007).

Regulatory mechanisms that provide protection from the deleterious effects of TEs have been selected in hosts. Among these, RNA interference (RNAi) is one of the major processes involved in TE control and acts at the transcriptional or posttranscriptional level through the action of small RNAs. In Drosophila, two RNAi pathways participate in TE control and rely on distinct RNA molecules: small interfering RNAs (siRNAs) and PIWI-interacting RNAs (piRNAs). They have their own RNA precursors and involve particular molecular effectors. RNAi is also known to be the main antiviral mechanism in insects (Obbard et al. 2009; Mussabekova et al. 2017).

piRNAs are 23 to 30 nucleotide-long sequences processed from single-stranded RNA precursors that are enriched in TE sequences and transcribed from genomic loci called piRNA clusters (Siomi et al. 2011; Senti et al. 2015; Suzuki et al. 2020). piRNAs are loaded onto PIWI proteins, members of the AGO protein family. The corresponding RNA-induced silencing complex (RISC) subsequently targets TE RNAs through sequence complementarity and induces the transcriptional or posttranscriptional silencing of TE sequences (Ozata et al. 2019). The Dicer-2 endonuclease plays a central role in the siRNA pathway. This protein slices double-stranded RNA (dsRNA) precursors into 21 nucleotides, and one strand is then loaded onto Argonaute 2, constituting the RISC that degrades any TE or viral sequence complementary to the loaded siRNA (Campbell et al. 2008; Palmer et al. 2018). In Drosophila melanogaster, in which they have been extensively studied, piRNAs are involved in TE control in both germline and gonadal somatic cells to prevent TE transposition, whereas siRNAs fight against TEs and viruses both in the soma and gonads (Siomi et al. 2011; Bronkhorst and van Rij 2014).

Given that the siRNA pathway targets both TEs and viruses, it seems legitimate to wonder how organisms cope when they are challenged from both sides. In particular, one can ask whether and how TE regulation is affected by virus infection. Two scenarios can be envisioned (Roy et al. 2020). The first one posits that there is a trade-off between TE control by the host and its ability to fight a viral infection so that an organism cannot achieve both effectively at the same time. If this is the case, then one would expect an increase in TE transcript amounts upon viral infection. This scenario is in line with the idea that the Dicer-2 protein, which is the siRNA producer, can be saturated (Durdevic et al. 2013). The second scenario involves a synergistic interaction between antiviral and TE silencing responses: viral infections can activate small RNA pathways and thus trigger a strong TE silencing response and conversely. In a recent study, Roy et al. (2020) investigated the effects of viral infections on TE transcript levels in Drosophila. They showed that during Sindbis virus (SINV) infection in Drosophila simulans and D. melanogaster carcasses, TE transcript and TE-derived small RNA modulation support the scenario of an amplification between TE control and antiviral immunity: upon infection, TE transcript levels are downmodulated, which is associated with an increase in TE-derived small RNAs. In a subsequent study, Roy et al. (2021) performed similar experiments and showed that when D. simulans was infected with D. melanogaster sigma virus (DMelSV), modulation resulted in a different scenario: TE transcripts levels increased in carcasses, in line with a trade-off scenario. Thus, depending on the Drosophila species and the virus, TE transcripts and associated small RNAs are differentially modulated. Data from other organisms are thus clearly needed to determine whether a general pattern emerges for the impact on both piRNA and siRNA pathways and to determine the extent to which TEs are modulated during viral infection. This comparative genomics approach is all the more needed given that D. melanogaster appears to be an exception rather than the rule among arthropods in terms of its very low abundance of somatic piRNAs (Lewis et al. 2018). Here, we focused on TE transcript modulation upon viral infection by DENV, ZIKV, and CHIKV in A. aegypti.

Aedes aegypti is a dipteran species of the Culicidae family, which diverged from Drosophilidae approximately 260 million years ago (Gaunt and Miles 2002). It is a major vector of arboviruses, such as Zika virus (ZIKV), yellow fever virus (YFV), dengue virus (DENV serotypes 1 to 4), and chikungunya virus (CHIKV). The first three viruses belong to the Flavivirus genus and have single-stranded positive and nonpolyadenylated RNA genomes of approximately 10 kb (Markoff 2003; Kuno and Chang 2007). These genomes encode a unique long open reading frame (ORF) corresponding to three structural proteins and seven nonstructural proteins (Chambers et al. 1990). CHIKV is an Alphavirus of the Togaviridae family that has a single-stranded positive and polyadenylated RNA genome (Strauss and Strauss 1994; Solignat et al. 2009). It encodes two ORFs corresponding to five structural proteins and four nonstructural proteins (Solignat et al. 2009). In addition, some viral products act as viral suppressors of RNAi, such as the nonstructural protein NS2A of DENV-2 (Qiu et al. 2020), NS4B of all DENV serotypes (Kakumani et al. 2013), or NSP2 and NSP3 of CHIKV (Mathur et al. 2016). Flaviviruses are also known to produce subgenomic flavivirus RNAs (sfRNAs) that antagonize the antiviral response of mosquitoes (Schnettler et al. 2012, Moon et al. 2015, Göertz et al. 2019). In addition to being a major vector of arboviruses, A. aegypti is known for its large genome, which carries many TEs (∼47%) (Nene et al. 2007). The most abundant TEs are non-long terminal repeat (LTR) retrotransposons, followed by DNA transposons and LTR retrotransposons (Matthews et al. 2018). Some authors have proposed that this ability to vector viruses and this high TE content are linked with a particular piRNA pathway in Aedes mosquitoes. Indeed, A. aegypti appears to be an exception among arthropods for producing antiviral piRNAs (Hess et al. 2011; Morazzani et al. 2012; Miesen et al. 2015, 2016; Lewis et al. 2018). In addition, in the Aedes genus, there has been an expansion of PIWI genes. Indeed, Campbell et al. (2008) found eight proteins of the PIWI subfamily in A. aegypti with specific localizations and functions, compared to only three in Drosophila (supplementary fig. S1, Supplementary Material online). The authors propose that this expansion could be adaptive for the control of a larger TE burden. Miesen et al. (2015) identified four distinct groups of A. aegypti TE families based on the different PIWI proteins involved in their regulation. In addition, Varjak et al. (2017) showed that Piwi4 interacts with effectors of both the siRNA and piRNA pathways in A. aegypti, establishing a link between these pathways that has remained elusive in Drosophila. Moreover, the A. aegypti genome also contains endogenous viral elements (EVEs), which are remnants of ancient viral infections and behave as piRNA clusters (Whitfield et al. 2017).

For these reasons, we investigated the impacts of viral infections on TE control in A. aegypti. The aim of the present study was to test whether DENV, ZIKV, or CHIKV infection affects TE transcript levels and whether the impacts vary according to the virus and tissue. We exploited RNA-seq data already available in public databases. Our analyses show that these viral infections have only a modest effect on TE transcripts and that the direction of the modulations varies greatly across viruses, tissues, and TE superfamilies. The datasets also allowed us to test the effect of temperature on TE transcript levels upon viral infection and thus to determine the magnitude of the impact of these factors.

Results

TE Transcript Levels Are Modulated upon Arboviral Infection in A. aegypti

We analyzed 11 publicly available RNA-seq datasets, corresponding to 371 samples, which covered a wide variety of tissues (whole body, midgut, ovary, and cell culture), incubation times (from 3 h postinfection [hpi] to 14 d postinfection [dpi]), mosquito strains and serotypes for three viruses (ZIKV, DENV, and CHIKV) (Table 1). Across the different samples, the ratios of TE counts to gene counts were consistently very low, ranging from 0.01 to 0.05, except 0.11 for PRJNA545086 (ZIKV, whole bodies) (supplementary table S1, Supplementary Material online). For each modality, which corresponds to a single combination of viruses, tissues, time after infection, and additional factors, we computed the log2-fold change (lg2fc) of the normalized read counts attributed to TEs in the infected samples compared to those in the corresponding noninfected samples (mock) (supplementary fig. S2 to S16, Supplementary Material online). As an example, we show the results obtained upon DENV-2 infection in samples prepared from whole bodies at 18 hpi (Fig. 1). We detected a global decrease in TE transcript levels upon infection for this particular modality (mean lg2fc = −0.41, Wilcoxon paired test, V = 920,198, P = 1e−54). Separate analyses splitting TEs into three categories (DNA transposons, LTR retrotransposons, and non-LTR retrotransposons) revealed consistent patterns: mean lg2fc = −0.45 for DNA transposons (Wilcoxon paired test, V = 8,068, P = 6e−17), −0.41 for LTR retrotransposons (Wilcoxon paired test, V = 442,841, P = 1e−27), and −0.37 for non-LTR retrotransposons (Wilcoxon paired test, V = 39,807, P = 1e−14) (Fig. 1). This global decrease is, however, associated with a large dispersion of lg2fc values across TE families, especially for LTR and non-LTR retrotransposons.

Fig. 1.

TE transcript modulation upon DENV-2 infection in samples prepared from whole bodies at 18 hpi, PRJNA476553 (Kang et al. 2018). Log2 of the ratio (infected counts/mock counts) for each annotated TE family of DNA transposons (left panel), LTR retrotransposons (central panel) and non-LTR retrotransposons (right panel). P-values are for Wilcoxon paired tests. Dashed lines correspond to the mean of all TE families.

Table 1

Summary of the analyzed datasets

Virus	Tissue	Time post infection	Mosquito strain	Viral strain	Number of analyzed samples	Accession number	Reference
ZIKV	Midgut	1 or 2 dpi	Chiapas, Mexico; field sampling	ZIKV MEX1-44	36	PRJNA615972	Ferreira et al. (2020)
		2 dpi	California	ZIKV BR15	11a	PRJNA818687	Louie et al. (2022)
		7 dpi	Rockefeller	ZIKV FSS13025	3b	PRJNA379149	Angleró-Rodríguez et al. (2017)
	Whole bodies	10 dpi	Colony from Australia	ZIKV MR766	9	PRJNA545086	Slonchak et al. (2020)
		2, 7 and 14 dpi	Galveston	ZIKV Mex 1-7	18	PRJNA399504	Etebari et al. (2017)
DENV	Midgut	1 or 4 dpi	Wild type population from Thailand	DENV-1 KDH0030A	47	PRJNA386455	Raquin et al. (2017); Suzuki et al. (2017)
		7 dpi	Rockefeller	DENV-2 NGC	6	PRJNA379149	Angleró-Rodríguez et al. (2017)
	Ovaries	3 dpi	Rockefeller	DENV-2 NGC	6	PRJNA786000	Feitosa-Suntheimer et al. (2022)
	Whole bodies	3 or 18 hpi	Trinidad field isolate	DENV-2 JAM1409	24	PRJNA476553	Kang et al. (2018)
CHIKV	Aag2 cells	1 dpi	Aag2 cells	CHIKV IND-10-DEL1 E1	4	PRJNA802350	Dubey et al. (2022)
		2 dpi	Aag2 cells	CHIKV LS3	12	PRJNA885496	Rosendo Machado et al. (2022)
	Whole bodies	3 or 7 dpi	Colony from Australia	CHIKV Asian genotype	70	PRJNA630779	Wimalasiri-Yapa et al. (2021)

dpi: days postinfection, hpi: hours postinfection.

^aOnly antibiotic treatments during the pupal stage and blood-fed samples were analyzed.

^bBoth DENV- and ZIKV-infected samples from Angleró-Rodriguez et al. (2017) shared the same mock condition samples, which are only indicated for DENV infection in the table.

We performed the same analyses for all considered datasets (Fig. 2). Wilcoxon paired tests revealed significant, global modulations of TE transcript amounts in most conditions (after Bonferroni correction for multiple tests), but in directions that varied according to the virus and tissue. The mean lg2fc ranged from −0.45 to 0.79 when all TEs were considered simultaneously. The clearest pattern was observed for DENV infections, which induced either no effect or a decrease in TE transcript amounts (mean lg2fc: from −0.41 to 0.06). The greatest increase in TE transcript levels was observed upon CHIKV infection in whole-body samples (mean lg2fc: from 0.24 to 0.79), while CHIKV infection in Aag2 cell cultures was associated with a decrease in TE transcript levels (mean lg2fc: from −0.45 to −0.05). Regarding ZIKV infections, the results showed greater contrast, and the directions of the modulations were highly variable across studies and samples (mean lg2fc: from −0.22 to 0.24) (Fig. 2).

Fig. 2.

Mean TE transcript modulation in all analyzed datasets. The different shades in the background allow to distinguish the different studies. The horizontal line indicates lg2fc = 0, i.e. absence of modulation. Stars at the top (*) indicate significant global TE transcript modulation (Wilcoxon paired tests, Bonferroni correction). Wh. bodies = whole bodies, Ov. = ovaries, Midg. = midgut, Aag2 c. = Aag2 cells, antib. = antibiotic treatment, sfRNA-def. = sfRNA-deficient ZIKV, stress = larval-stage crowding stress, Dhx15 = Dhx15 knockdown. ^1. Ferreira et al. 2020. ^2. Louie et al. 2022. ^3. Angleró-Rodriguez et al. 2017. ^4. Etebari et al. 2017. ^5. Slonchak et al. 2020. ^6. Feitosa-Suntheimer et al. 2022. ^7. Suzuki et al. 2017; Raquin et al. 2017. ^8. Kang et al. 2018. ^9. Dubey et al. 2022. ^10. Rosendo Machado et al. 2022. ^11. Wimalasiri-Yapa et al. 2021.

TE Transcript Levels Are More Sensitive to Temperature Changes than to Arboviral Infections

Some RNA-seq datasets included samples corresponding to different temperatures (for ZIKV and CHIKV). In these different studies, RNA-seq data were obtained at optimal mosquito temperature (28 °C) or at cold (either 18 °C or 20 °C) or hot (either 32 °C or 36 °C) temperatures. Using an ANOVA approach on total TE counts, we found no interaction between temperature and infection status (supplementary table S2, Supplementary Material online). Temperature had a significant effect in all studies (two-way ANOVA P-values for temperature effect: ZIKV 1 dpi: 1.21e−04, 2 dpi: 9.97e−06, CHIKV 3 dpi: 1.21e−08, 7 dpi: 6.53e−05; Fig. 3 and supplementary table S2, Supplementary Material online). Cold temperatures were always found to be associated with the highest total TE counts. In contrast, these ANOVAs of total TE counts revealed a significant impact of the infection status in CHIKV samples only (two-way ANOVA P-values for infection effect: CHIKV 3 dpi: 5.06e−04, 7 dpi: 7.13e−05). The ratios of the effect sizes of temperature to infection status were always greater than 1, indicating that temperature had a much greater effect than infection status (effect size ratio: ZIKV 1 dpi: 5.65; ZIKV 2 dpi: 4.77; CHIKV 3 dpi: 2.64; and CHIKV 7 dpi: 1.37).

Fig. 3.

Effects of temperature and infection status on TE transcript levels. Sum of normalized TE counts. Red points correspond to means across replicates. Cold = 18 °C (CHIKV) or 20 °C (ZIKV), Optimal = 28 °C (CHIKV and ZIKV), Hot = 32 °C (CHIKV) or 36 °C (ZIKV). A) ZIKV infection, midgut, 1 dpi, B) 2 dpi (Ferreira et al. 2020), C) CHIKV infection, whole bodies, 3 dpi, and D) 7 dpi (Wimalasiri-Yapa et al. 2021).

TE Responsiveness to Viral Infection Varies across Superfamilies

Although mosquito strains differed across datasets and TE landscapes are known to vary across natural populations, we investigated whether the patterns of differentially expressed (DE) TE families were consistent across viruses. Because small-RNA pathways target TE RNAs through sequence complementarity (Campbell et al. 2008; Palmer et al. 2018; Ozata et al. 2019), we assume that TE families belonging to the same superfamily are regulated by related mechanisms and thus share the same sensitivity to homoeostasis perturbations. Hence, we tested whether the different TE superfamilies displayed similar distributions of DE versus non-DE families upon infection across the different conditions. We detected significant differences among superfamilies for the three viruses (chi-squared tests, ZIKV infection: chi-squared = 210.65, simulated P < 4.8e−10; DENV infection: chi-squared = 93.49, simulated P = 6.0e−08; and CHIKV infection: chi-squared = 121.73, simulated P < 4.8e−10) (Fig. 4A–C). The largest contributions to the chi-squared statistics came from particularly infection-responsive superfamilies, i.e. more often DE than expected. Moreover, DNAt (superfamily), Loa, RTE, and Sola2 had more DE TE families than expected, whereas families from BEL superfamilies were less frequently DE than expected upon infection with all three viruses. Some superfamilies displayed contrasting patterns across viruses. For instance, Crack TE families were particularly insensitive to flavivirus infections (DENV and ZIKV). CR1 displayed an increased number of DE TE families upon CHIKV infection. TE families from Mariner/Tc1 were more often DE than expected upon flavivirus infection. The same was true for R1 upon DENV or CHIKV infection. In contrast, the numbers of DE TE families from Gypsy, I, and Tx1 followed random expectations in all conditions (Fig. 4A–C). In conclusion, it appears that there is a strong interaction between the TE superfamily and the virus in terms of responsiveness to infection.

Fig. 4.

Observed and expected numbers of differentially expressed TE families within each TE superfamily upon infection with A) ZIKV, B) DENV, and C) CHIKV. Total: Total number of observed TE families (DE and not DE); Obs: Observed number of DE TE families; Exp: Expected number of DE TE families under the homogeneity hypothesis; C: Contribution to chi-squared statistics under the null hypothesis of a homogeneous distribution of DE and non-DE TE families, i.e. (Obs−Exp)²/Exp; R: Relative deviation from homogeneity, i.e. (Obs−Exp)/Exp.

As families could be up- or downmodulated upon infection in different studies, we further investigated the consistency of the response of TE superfamilies across the different modalities by plotting the mean of the mean lg2fc across the different modalities over the variation coefficient of the mean lg2fc across studies for each superfamily. This allowed us to spot superfamilies displaying a strong response in the same direction in all the modalities separately for each virus. Because some superfamilies were found only in a few samples, they could be associated with a small variation coefficient, as was the case for the P element (Fig. 5A). In order to display this bias, we adjusted the size of the dots according to the presence of the TE superfamily across datasets. The global TE superfamily modulation patterns differed across the three viruses, in agreement with our observations (Fig. 2). Considering the most represented TE superfamilies across datasets, Zator, Kolobok, and Harbinger were particularly upmodulated and L2B, Jockey, Outcast, Crack, and ISL2EU were the most consistently downmodulated superfamilies between modalities upon ZIKV infection (Fig. 5A). With regard to DENV infection, Helitron and ISL2EU were globally the most upmodulated TE superfamilies across datasets, whereas Crack, MITE, Outcast, or Zator were generally downmodulated (Fig. 5B). Finally, CryptonI, Nimb, Tx1, Helitron, and Kiri were mostly upmodulated upon CHIKV infection in all the datasets (Fig. 5C).

Fig. 5.

Mean of the TE superfamily mean lg2fc across the different modalities according to the variation coefficient. A) ZIKV, B) DENV, and C) CHIKV. The larger the dot is, the greater the number of families per superfamily. TE superfamilies most consistently responsive to infection across datasets are in bold.

Discussion

TE Transcripts Are Modulated upon Viral Infection

Overall, the analyses that we performed on a collection of diverse publicly available RNA-seq datasets allowed us to uncover significant effects of arboviral infections on TE transcript levels in A. aegypti. We found that TE transcript levels were modulated upon arboviral infection in most studies. Nevertheless, the mean lg2fc values were lower than one, suggesting that the biological impacts of such modulations are expected to be mild, especially since TE read counts account for only a small fraction of the total gene counts (supplementary table S1, Supplementary Material online).

Although direct quantitative comparisons are difficult to perform, we recall that our previous studies of TE transcripts in Drosophila revealed mean values of −0.25 and −0.32 for lg2fc upon SINV infection in D. simulans and D. melanogaster, respectively (Roy et al. 2020), and 0.32 upon DMelSV infection in D. simulans (Roy et al. 2021). Apart from the CHIKV data, the lg2fc values that we measured here were generally lower, between −0.25 and 0.25 (Fig. 2), which suggests stronger TE transcript homeostasis upon viral infection in A. aegypti than in Drosophila. This could be related to the amplification of piRNA pathway genes in the mosquito genome (Campbell et al. 2008) and their involvement in both TE control and antiviral immunity (Morazzani et al. 2012; Lewis et al. 2018). Alternatively, it is possible that we did not detect any strong impact of arboviral infections on TE transcript amounts simply because arboviral infections globally have a modest impact on the mosquito transcriptome, which is consistent with the fact that arboviruses are known to have only a modest effect on mosquito fitness (Lambrechts and Scott 2009). This idea is reinforced by our analysis of the effect of temperature changes, which have greater impacts on the TE transcriptome than the arboviral infection status. Our study revealed contrasting TE transcript modulations and sensitivities depending on viral infection.

DENV Infection

The clearest pattern was observed for DENV infections; in most studies, DENV infection was associated with TE downmodulation (Fig. 2). Because DENV has tissue tropism for the midgut (Salazar et al. 2007; Angleró-Rodríguez et al. 2017), we expected stronger TE modulation in this tissue. However, this is not what we observed. This may be related to the recent observation that the siRNA pathway controls systemic DENV and ZIKV dissemination in mosquitoes but not in the midgut, although the pathway is fully functional against TEs in this tissue (Olmo et al. 2018). This global decrease in TE transcript levels upon DENV infection is similar to the pattern we observed in Drosophila upon SINV infection, for which we proposed that the antiviral response led to an increase in TE control (Roy et al. 2020). Indeed, the dsRNA uptake pathway was demonstrated to participate in the systemic spread of the RNAi antiviral response (Saleh et al. 2009), and we proposed that dsRNA molecules corresponding to TE sequences—resulting from antisense transcription or from inverted repeat pairing—hijack the process, leading to a concomitant increase in virus-derived siRNAs and TE-derived small RNAs. Interestingly, the study by Kang et al. (2018) revealed lower viral replication in individuals subjected to larval-stage crowding stress, which we found to be associated with a less pronounced decrease in TE transcript levels. This highlights the concomitant variation of the host antiviral response and TE control.

CHIKV Infection

With regard to CHIKV infections, TE transcripts were globally downmodulated in somatic Aag2 cells but clearly upmodulated in whole-body samples, with the highest observed mean lg2fc of all datasets. This discrepancy between the two studies is difficult to explain. Cell cultures are known to be subject to frequent genome rearrangements (Lee et al. 2014), which can have impacts on TE landscapes. CHIKV data in whole bodies all come from the same study, and results from additional studies are needed to further confirm our conclusion. However, these findings are in line with what we previously observed in Drosophila upon DMelSV infection: we observed a global increase in TE transcript levels in fly carcasses. The underlying molecular mechanisms are not understood, but a scenario involving a trade-off between antiviral immunity and TE control may be proposed.

ZIKV Infection

ZIKV infection led to more variable patterns across datasets, and we observed that TE transcript levels decreased as well as increased, depending on the study. Thus far, it has been very difficult to identify the factors determining the direction of the modulation. Moreover, the different viral strains used to infect mosquitoes all originated from the Asian lineage (Faria et al. 2016; Lanciotti et al. 2016). In comparative studies, strains from the Asian lineage presented more similar infection characteristics than strains from the African lineage (Roundy et al. 2017; Ou et al. 2021). Thus, even if mosquitoes were infected by different ZIKV strains, it cannot be used as a strong argument to explain the variability of the modulation observed (Fig. 2).

Some TE Superfamilies Are Particularly Infection Responsive

We analyzed in detail the patterns of TE modulation upon viral infection at the superfamily level. Indeed, families belonging to the same superfamily display some degree of structural and sequence similarity, enabling common regulatory mechanisms to some extent. Across the different studies, we found that the TE superfamilies DNAt, Loa, RTE, and Sola2 were the most often impacted by the three viral infections compared to other superfamilies, whereas BEL was very insensitive to arboviral infections (Fig. 4A–C). In addition, our results showed that some superfamilies displayed particularly homogeneous modulations across all the corresponding families in different studies of the same virus, as is the case for Helitron upon DENV or CHIKV infection or L2B upon ZIKV infection (Fig. 5A–C). However, we cannot exclude the possibility that this is due to the lower number of families in these superfamilies. Indeed, a large number of families increases the probability of detecting heterogeneity across study modalities.

Variations Across Studies Are Explained by Various Biotic Factors

Similar to what we described in Drosophila, our results clearly show that the direction and extent of TE transcript modulation depend on the virus. Lambrechts and Scott (2009) showed that Alphaviruses have a greater impact on mosquito fitness than other viral genera, such as Flaviviruses. Consistent with this, in the present analysis, we found that the Alphavirus (CHIKV) had the strongest impact on TE transcript modulation. Finally, the study by Angleró-Rodriguez et al. (2017) allows to compare the impacts of DENV-2 and ZIKV on TEs using the same Rockefeller mosquito strain, and the same infection protocol at 7 dpi. While their data did not allow us to detect any significant effect on TEs upon DENV infection, we observed an upmodulation upon ZIKV infection. Even though ZIKV and DENV are Flaviviruses with shared features, Williams et al. (2020) showed that A. aegypti is less resistant to ZIKV infection than to DENV-2 in terms of viral suppression effectiveness, possibly because ZIKV is less susceptible to the siRNA pathway. We assume that this lower resistance to the virus can be reflected in various host responses, such as stronger modulation of TE transcripts.

As described above, the tissue under consideration may be relevant. Indeed, (i) depending on viral tropism, viral replication is not homogeneous in the mosquito body, and thus, the interaction with TE control, if any, is expected to vary, and (ii) the expansion of the PIWI genes in A. aegypti, as previously mentioned, is accompanied by both spatial and functional specializations of the different genes in the mosquito body, e.g. TE control, antiviral immunity or both. Piwi 1 to 3 are specific to the germline, Piwi 4 to 6 and Ago3 are specific to adult somatic tissues, and Piwi 7 is localized in early embryos (Akbari et al. 2013; Miesen et al. 2015) (supplementary fig. S1, Supplementary Material online). Thus, according to the tissue and function of the Piwi protein, the interaction between TE control and antiviral immunity, if any, is expected to differ. As mentioned above, flaviviruses are known to produce sfRNAs, which facilitate infection, suppress antiviral RNAi, and are required for efficient viral transmission (Schnettler et al. 2012; Moon et al. 2015; Göertz et al. 2019; Slonchak et al. 2020). Thus, it is possible that upon DENV and ZIKV infection, TE-derived piRNAs and siRNAs are also impacted by sfRNAs, which subsequently induces changes in TE transcripts. In line with this idea, the results obtained in whole-body samples at 10 dpi by mutated or not ZIKV for sfRNA production show that sfRNAs could have an impact on TE transcript regulation because when the virus cannot produce sfRNAs, the downmodulation of TE transcripts is milder (Fig. 2) (Slonchak et al. 2020).

Limitations

We have to acknowledge that the present approach suffers from several limitations. In particular, while the analysis of publicly available datasets allows the sensible use of research efforts and the exploration of a wealth of data, the corresponding experimental setups may not be perfectly suited to all research questions. Most notably, cross-study comparisons are not straightforward since many factors can vary or lack control. Of major impact, mosquito genotypes and virus strains are known to determine the outcome of arboviral infections (Lambrechts 2011; Lambrechts et al. 2013; Lequime et al. 2016; Souza-Neto et al. 2019; Ekwudu et al. 2020). In addition, considering that we found that the different TE superfamilies are not evenly sensitive to viral infections and that the TE content is known to vary across natural populations (Vieira et al. 1999, 2002; Goubert et al. 2017; Mérel et al. 2021), it is expected that different mosquito genotypes will display different modulations of their TE transcriptomes depending on the TE composition of their genomes. The tissue under study is also of fundamental importance since viral replication does not occur homogeneously within the whole body; therefore, the molecular interactions with TEs may differ. Moreover, some studies have been performed in cultures of cells, whose genomes are known to suffer frequent rearrangements (Lee et al. 2014), which may have an impact on TE distribution (Potter et al. 1979; Ilyin et al. 1980). In addition, the infection status with regard to bacterial or chronic viral infections is not known in most studies, but it was shown that they could interfere with the antiviral response (Ramirez et al. 2012; Parry and Asgari 2018; Souza-Neto et al. 2019; Olmo et al. 2023). Wolbachia bacteria are known to reduce the replication of many viruses by, e.g. (i) modifying gene regulation (Zhang et al. 2013) or (ii) increasing the level of piRNAs and Ago3 transcripts (Mayoral et al. 2014). Current evidence suggests that A. aegypti does not harbor Wolbachia naturally (Ross et al. 2020); however, the bacteria invaded some A. aegypti populations after release in the wild of artificially infected mosquitoes (Moreira et al. 2009; Hoffmann et al. 2011, 2014; Walker et al. 2011; Ahmad et al. 2021). The infection status of the populations under study is not known, and even though it is highly unlikely, some of them may have been affected because certain releases were carried out before mosquito sampling and in nearby geographical areas, as is the case for Queensland in Australia (PRJNA630779 in Wimalasiri-Yapa et al. (2021) and PRJNA545086 in Slonchak et al. (2020)). Other parameters, such as the route of infection (Mondotte and Saleh 2018) or the viral dose, may also impact the outcome of the infection and therefore the potential effect on the TE transcriptome. However, despite these limitations, the large number of analyzed samples allowed us to identify general trends in the interaction between arboviral infections and TE control.

Conclusion

The results of our analyses show that TE transcript levels are modulated upon arboviral infection in A. aegypti. However, we revealed that the direction of the modulation depends on the virus and potentially on some unidentified factors. Overall, the extent of the modulations that we found appears to be particularly low, raising doubts about their biological impacts. However, further studies are needed to evaluate this phenomenon. Nevertheless, even subtle changes in TE transcript levels may result in differences in TE transposition rates, and, therefore, in somatic genetic diversity (Faulkner and Garcia-Perez 2017; Chang et al. 2019), or even in genome mutation rates if the ovarian transposition rate is affected.

Moreover, it would be particularly interesting to test whether TEs have an impact on viral replication in mosquitoes. This idea is suggested by recent research that revealed the role of EVEs in mosquitoes. These are remnants of ancient viral infections, which are acquired through recombination with LTR retrotransposons, behave as piRNA clusters, and participate in viral immunity against exogenous viruses (Whitfield et al. 2017). EVEs are known to vary across mosquito populations and are proposed to explain differences in vector competence (Suzuki et al. 2017; Whitfield et al. 2017). The obvious similarities between EVEs and TEs suggest that this is a relevant direction for future research.

Materials and Methods

Datasets

We retrieved publicly available RNA-seq datasets corresponding to A. aegypti samples infected with DENV, CHIKV, or ZIKV and their corresponding noninfected samples using the prefetch and fastq-dump tools of the sra-toolkit (http://github.com/ncbi/sra-tools/wiki/01.Downloading-SRA-Toolkit). All datasets were produced following Illumina procedures. We only retained datasets that displayed at least three biological replicates. Sequencing quality was assessed using the fastqc-0.11.9 tool (www.bioinformatics.babraham.ac.uk/projects/fastqc/), and poor-quality samples were removed (SRR17843159 and SRR17843156 in PRJNA802350). When the datasets had technical replicates, they were concatenated into a single fastq file. SRR18435178 (PRJNA818687), SRR11440744 (PRJNA615972), and SRR6818585 and SRR6820666 (PRJNA386455) were also removed because they were potentially corrupted.

We used version 5.0 of the A. aegypti reference genome (AaegL5.0) obtained by long-read sequencing via PacBio technology (Matthews et al. 2018). This reference genome corresponds to the LVP_AGWG strain of A. aegypti. After masking TEs using RepeatMasker-4.1.2-p1 (http://www.repeatmasker.org) to remove them from genes, we retrieved the gff file corresponding to gene models in AaegL5.0 (GCF_002204515.2_AaegL5.0_genomic.gff.gz) and used getfasta-2.30.0 from bedtools (Quinlan and Hall 2010) to create a multifasta file of gene sequences devoid of TEs. For the analysis of TEs, we considered A. aegypti sequences retrieved from the Repbase-25.08 database of eukaryotic repetitive elements (Jurka et al. 2005), corresponding to 3,013 TE families classified into 40 superfamilies. We formatted these sequences to analyze the TE data using the TEcount module of TEtools (Lerat et al. 2017).

RNA-seq Analysis

Sequencing adaptors were removed using Trimmomatic version 0.39 (Bolger et al. 2014). Alignments against the masked multifasta of reference genes and gene counts were performed using kallisto-0.48.0 (Bray et al. 2016). TE transcript counts were obtained using the TEcount module of TEtools (Lerat et al. 2017) and the list of TE sequences described above.

Gene and TE count tables were concatenated and further analyzed using the R package DESeq2-1.34.0 (Love et al. 2014). For each study and tissue, normalization was done globally on the concatenated table using the DESeq function. We considered a minimum threshold of ten normalized counts for TEs in both infected and mock conditions. Global TE transcript modulation was assessed by Wilcoxon paired tests at the TE family or superfamily levels across mean counts in infected and control conditions. Bonferroni correction was applied across the different datasets.

The effect of temperature on the log-transformed sums of the normalized TE counts of all TE families was assessed using two-way ANOVA, and the factors temperature, infection status, and temperature × infection status interaction were considered. The assumptions of normality and homogeneity were met (Shapiro, Bartlett, and Levene tests), except for samples corresponding to CHIKV infection in whole bodies at 3 dpi. The effect size of temperature was determined by averaging the means of both the control and infected conditions for each temperature, followed by computing the difference between the two most extreme means. Similarly, the effect size of infection was assessed across all three temperatures by computing the difference between the averages of all means from the infected condition on one side and all means from the mock condition on the other side.

TE responsiveness to viral infection was achieved after extracting TE families DE or not between infected and mock conditions from DESeq tables. A TE family was considered non-DE if it had an adjusted P > 0.05 across all modalities in the DESeq2 analysis; it was considered DE if it showed an adjusted P < 0.05 in at least one modality. Chi-squared tests for homogeneity with Monte Carlo simulations for the P-value performing 2.1e + 09 iterations were used to assess whether the proportions of DE and non-DE families were the same in all the superfamilies.

The consistency of TE superfamily responses across the different study modalities was further analyzed by computing the mean lg2fc for each TE superfamily expressed in each modality after aggregating the TE family lg2fc. Only families displaying at least ten counts in one condition (infected or mock) were kept for analysis. Variation coefficients were computed from these means for each study modality. The global mean lg2fc for each TE superfamily across all the modalities was also computed. Thus, high magnitude (either positive or negative) for mean lg2fc and a low variation coefficient of lg2fc corresponded to TE superfamilies displaying a strong response in the same direction across all the modalities.

Supplementary Material

Supplementary material is available at Genome Biology and Evolution online.

Acknowledgments

This work was performed using the computing facilities of Laboratoire de Biométrie Biologie Évolutive/Pôle Rhône-Alpes de Bioinformatique (CC LBBE/PRABI). We thank Philippe Veber, Natacha Kremer, Séverine Chambeyron, Guillaume Minard, François Sabot, Olivier Terrier, Alicia Reyes Ramirez, Cristina Vieira, William Vilas Boas Nunes, Camille Mayeux, Daniel Siqueira de Oliveira, Miriam Merenciano, and Anaïs Larue for helpful discussions. We also thank Clément Goubert for providing technical assistance. We thank the editor for her positive feedback that helped us improve the manuscript.

Funding

This work was supported by Agence Nationale de la Recherche (ANR) LongevitY ANR-20-CE02-0015, MosquiTEs ANR-21-CE02-0013, and TEMIT ANR-19-CE02-0022 grants and Universite Claude Bernard Lyon1.

Data Availability

There are no new data associated with this article. The scripts used for the analyses are available in Supplementary material.

Author notes

Contributed equally