Infections by Transovarially Transmitted DMelSV in Drosophila Have No Impact on Ovarian Transposable Element Transcripts but Increase Their Amounts in the Soma.

Transposable elements (TEs) are genomic parasites, which activity is tightly controlled in germline cells. Using Sindbis virus, it was recently demonstrated that viral infections affect TE transcript amounts in somatic tissues. However, the strongest evolutionary impacts are expected in gonads, because that is where the genomes of the next generations lie. Here, we investigated this aspect using the Drosophila melanogaster Sigma virus. It is particularly relevant in the genome/TE interaction given its tropism to ovaries, which is the organ displaying the more sophisticated TE control pathways. Our results in Drosophila simulans flies allowed us to confirm the existence of a strong homeostasis of the TE transcriptome in ovaries upon infection, which, however, rely on TE-derived small RNA modulations. In addition, we performed a meta-analysis of RNA-seq data and propose that the immune pathway that is triggered upon viral infection determines the direction of TE transcript modulation in somatic tissues.

Significance

Transposable elements (TEs) are genomic parasites that are found in all genomes. However, the determinism of their activity is still incompletely known. Here we show that infections by the Drosophila melanogaster Sigma virus (DMelSV) increase somatic TE transcript amounts. A survey of other RNA-seq data sets suggests that somatic TE transcript modulation upon viral infection is a general rule and that the direction of the modulation depends on the host immune pathways that are triggered. In addition, although DMelSV is known to display tropism to ovaries, this study uncovers a strong homeostasis of ovarian TE transcript amounts upon infection, whereas TE-derived small RNAs increase. This potential modulation of the maternally transmitted TE-derived small RNA repertoires could have evolutionary impacts on genomes.

Introduction

Both quantitatively and qualitatively, transposable elements (TEs) are major components of genomes. Quantitatively, they are repeated sequences that can make a few percents of chromosome sequences up to 90% in some plants (Biémont 2010; Wells and Feschotte 2020). Qualitatively, they are mobile sequences that may modify DNA sequences and epigenetic landscapes (Rebollo et al. 2012; Lee and Karpen 2017; Quadrana et al. 2019). Research work of the last decade showed that one of the ways of controlling TE activity relies on epigenetic mechanisms (Slotkin and Martienssen 2007). In insects, these are essentially RNA interference (RNAi) pathways. In Drosophila, the activity of piRNAs (PIWI-interacting small interfering RNAs) is predominant in gonads. Drosophila ovaries are composed of germline tissue and somatic support cells; piRNAs and functional piRNA pathways are present in both ovarian tissues, the germline and somatic follicle cells. piRNAs are 23–30 nt-long single-stranded RNA molecules that recognize TE transcripts by sequence similarity, and trigger their slicing thanks to the RNAse activity of the PIWI proteins they are associated with. Some of these loaded piRNAs also translocate into the nucleus and induce the heterochromatic silencing of the corresponding TE sequences (Siomi et al. 2011; Czech and Hannon 2016). In somatic tissues, the small interfering RNA (siRNA) pathway is also at play in the transcriptional control of TEs. In addition, this RNAi pathway is the first defense line against viral infections (Galiana-Arnoux et al. 2006). siRNAs are 21-nt-long single-stranded RNA molecules, which are processed by Dicer-2 from double-stranded viral RNA (dsRNA) intermediates. These siRNAs are then loaded onto the Ago2 protein, which cleaves viral fragments sharing sequence complementarity to the siRNA (Galiana-Arnoux et al. 2006; van Rij et al. 2006; Mussabekova et al. 2017).

TE transcriptional and transpositional activities vary a lot across the tree of life (Wells and Feschotte 2020), and the determinisms of these activities are still incompletely known. Using a Drosophila system, we have recently shown that viral infection modulates somatic TE transcript amounts, through changes of TE-derived small RNA repertoires (Roy et al. 2020). Infections of wild-type flies by the Sindbis arbovirus (SINV) resulted in a decrease in somatic TE transcript amounts. However, infections of loss-of-function mutants of the siRNA pathway induced a strong increase in somatic TE transcript amounts. In both cases, ovarian TE transcript amounts were virtually unmodulated. The impacts of viral infections on germline TE transcripts is still a matter of debate because SINV does not enter ovaries (Roy et al. 2020). Indeed, contrary to somatic TE modulation, only germline TE modulation is expected to have evolutionary genomic consequences. These are the reasons why we used the same Makindu Drosophilasimulans strain and infected it using another virus: the Drosophila melanogaster Sigma virus (DMelSV). DMelSV is a natural pathogen of D.melanogaster, and it is one of the few characterized Drosophila viruses known to infect ovarian tissues (L’Héritier 1958; Jousset and Plus 1975; Louis et al. 1988), which is the very place where TEs are strongly silenced, due to the piRNA pathway. DMelSV belongs to the Rhabdoviridae. It is an enveloped virus, displaying a single-stranded RNA genome of negative polarity. It was discovered in D. melanogaster (L’Heritier 1948), and found to naturally infect different species of Drosophila; however, D. simulans is not one of these (Webster et al. 2015). DMelSV is particular among Drosophila viruses in that it is vertically transmitted through eggs and sperm, whereas horizontal transmission does not seem to occur (L’Héritier 1958; Fleuriet 1988). Using a previously characterized D. simulans strain, we show that TE transcript amounts increase in somatic tissues upon DMelSV infection. We observe only very weak modulation of ovarian TE transcript amounts, reinforcing the idea that ovaries are strongly protected against modulations of TE transcript amounts. However, ovarian TE-derived small RNAs are significantly more abundant upon infection, suggesting that TE control may vary in the zygotes. In addition, we performed a meta-analysis of previously released RNA-seq data upon viral infection, and we uncover that the direction of the somatic TE modulation is virus-specific. We propose that it relies on the particular immune pathways that are triggered upon infection. Indeed, it is known that they vary a lot regarding the virus and the route of infection (Mondotte and Saleh 2018; Palmer, Varghese, et al. 2018).

Results

DMelSV Replicates in D. simulans

In order to go deeper into the understanding of the impact of viral infections on the germline activity of TEs, we considered DMelSV, which is one of the few characterized Drosophila viruses known to replicate in ovaries. In addition, in order to take advantage of our previous thorough analyses (Roy et al. 2020), we decided to use again the Makindu wild-type strain of D. simulans. As DMelSV is not naturally found in D. simulans (Webster et al. 2015), we first checked whether it replicated in Makindu. Similarly to the test used in D. melanogaster, we anesthetized flies using carbon dioxide (CO2), and we computed the frequency of infected flies as the frequency of flies which did not recover from anesthesia (within 24 h). All flies were dead 11 days after infection (fig. 1), indicating that DMelSV indeed replicated in this experimental system, and that 100% of the flies were successfully infected. The presence of DMelSV genome in the ovaries of infected flies was confirmed using RT-PCR, up until 14 days post infection (dpi) (fig. 1). Eight days after infection, we could detect DMelSV-derived 21 nt-long small RNAs in carcasses and ovaries, although less abundantly in the later (fig. 1). Sense and antisense viral siRNAs were found in equivalent proportions (fig 1) and mapped all along the viral genome (supplementary fig. S1, Supplementary Material online). The size distribution of small RNA reads aligning to DMelSV indicated that they corresponded to 21 nt-long small RNAs in the majority, that is, siRNAs (fig. 1), as described previously by others (Petit et al. 2016).

Fig. 1

DMelSV infection in D. simulans. (A) Proportion of flies which recovered after CO2 anesthesia through time. In these experimental conditions, 100% of the flies were successfully infected with DMelSV. (B) DMelSV detection in ovaries of infected flies using RT-PCR. (C) Size distributions (nt) for DMelSV-mapping small RNAs. Raw read counts. Sense and antisense alignments are counted above and below the x axes, respectively. Double bars correspond to both sequencing replicates. As expected, the majority of these small RNAs are 21 nt long.

Detailed analyses of the small RNA-seq data revealed viral sequences other than DMelSV (supplementary fig. S2, Supplementary Material online). Similar to our previous work (Roy et al. 2020), Nora Virus (NoraV) reads were detected in mock as well as DMelSV-infected conditions, and supposed not to interfere with the immune answer to the experimental infection. We also identified Drosophila A virus (DAV) small RNA reads, although only in DMelSV-infected D. simulans tissues. DAV read numbers were almost two orders of magnitude larger than DMelSV read numbers. This indicates that Makindu flies are chronically infected with NoraV, but not with DAV. Small RNA-seq of DMelSV viral stock (i.e., the solution used to experimentally infect D. simulans flies) did not reveal DAV sequences. In addition, we could not detect DAV RNAs using RT-PCR in Makindu flies nor in the DMelSV viral stock (data not shown). It has also to be noticed that NoraV and DAV were hardly detectable from RNA-seq data. This tends to suggest that the flies used to prepare DMelSV samples were not chronically infected with DAV. DAV detection such as the one observed in our Makindu samples seems to be widespread in fly samples (Palmer et al. 2018), and we also noticed such occurrence in our previous experiments using SINV (Roy et al. 2020). Overall, even if we cannot unambigously trace back the origin of these DAV small RNAs, the above elements let us assume that they do not interfere significantly with DMelSV infection.

TE Modulation in Somatic Tissues and Ovaries upon DMelSV Infection

We extracted RNAs and small RNAs at 8 dpi to perform RNA-seq and small RNA-seq. Small RNAs displayed the expected size distribution, that is, 21–30 nt in length (supplementary fig. S3, Supplementary Material online), as well as the 1U enrichment expected for piRNAs and siRNAs (Brennecke et al. 2007; Wang et al. 2014) (supplementary fig. S4, Supplementary Material online). Detailed sequence composition of these small RNA repertoires are provided as supplementary table S5, Supplementary Material online.

In somatic tissues, we could detect a clear increase in TE transcript amounts upon infection (24% increase, mean log2FC = 0.31, Wilcoxon paired test P value = 1.0e-17, fig. 2). TE-derived small RNAs were also affected: TE-derived 23–30 nt RNAs decreased (47% decrease, mean log2FC = −0.92, Wilcoxon paired test P value = 2.1e-35, fig. 2) as well as TE-derived 21 nt RNAs (28% decrease, mean log2FC = −0.47, Wilcoxon paired test P value = 1.2e-21, fig. 2). All but three TE families had decreased TE-derived 23–30 nt small RNA levels. The three exceptions are Tc3, Copia-1-Dmel, and ISBU2 families, which appear to display very low levels of piRNAs, as well as siRNAs and transcripts; therefore, this discrepancy should not be considered as relevant. TE families displaying an increase in TE transcript amount upon infection were not significantly enriched in TE families displaying decreased TE-derived 21 nt small RNAs (Fisher exact test P value = 0.10).

Fig. 2

TE modulation upon DMelSV infection, expressed as log2 of the ratio (counts in infected condition/counts in mock condition) for each annotated TE family. Each dot is a TE family. (A) TE transcript modulation upon infection in carcasses. (B) TE-derived 23–30 nt small RNA modulation upon infection in carcasses. (C) TE-derived 21 nt small RNA modulation upon infection in carcasses. (D) TE transcript modulation upon infection in ovaries. (E) TE-derived 23–30 nt small RNA modulation upon infection in ovaries. (F) TE-derived 21 nt small RNA modulation upon infection in ovaries. Small RNA amounts were normalized relative to miRNAs. P values were obtained using Wilcoxon paired tests.

In ovaries, TE transcript amount changes upon infection are the weakest (12% decrease, mean log2FoldChange [log2FC] = −0.18, Wilcoxon paired test P value = 0.005, fig. 2), whereas TE-derived 23–30 nt RNAs clearly increased (19% increase, mean log2FC = 0.26, Wilcoxon paired test P value = 5.2e-37, fig. 2) as well as TE-derived 21 nt RNAs (18% increase, mean log2FC = 0.24, Wilcoxon paired test P value = 4.4e-23, fig. 2). Surprisingly, TE-derived small RNA modulations were in opposite directions in ovaries compared with somatic tissues.

Host Immune Response against DMelSV Infection

We further analyzed host response upon DMelSV infection focusing on a subset of genes that are characteristic of known immune pathways in Drosophila, that is, the RNAi, Toll, Imd, and Jak/STAT pathways (supplementary table S6, Supplementary Material online). In somatic tissues, we could detect significant increases in expression for genes of the Toll and Jak/STAT pathways (fig. 3). Statistically significant responses were observed in carcasses for PGRP-SA (log2FC = 1.91), drosomycin (log2FC = 1.50), GNBP-1 (log2FC = 0.86), Spatzle (log2FC = 0.76), Toll (log2FC = 0.57), and dorsal (log2FC = 0.39) (Toll pathway), and tep1 (log2FC = 1.91), vir-1 (log2FC = 1.46), totA (log2FC = 0.55), and STAT (log2FC = 0.49) (Jak/STAT pathway). However, the corresponding log2-fold changes are all lower than 2, indicating moderate overexpression. On the contrary, the Imd pathway and genes involved in RNAi did not display significant transcription increase upon infection. In ovaries, none of the genes from the RNAi (piRNA and siRNA), Toll, Imd, nor Jak/STAT pathways displayed significant differential expression upon infection at the 0.05 threshold for adjusted P values.

Fig. 3

Host immune response against DMelSV infection. Log2-fold changes of transcript amounts upon infection as estimated using DESeq2. Significant adjusted P values at the 0.05 threshold are indicated by stars.

Meta-analysis of Antiviral Host Responses and TE Modulation

In order to go deeper into the understanding of the mechanisms that are triggered upon viral infection and impact TEs, we performed a meta-analysis of available RNA-seq data (table 1). Apart from our previously published data (Roy et al. 2020), available data sets were not prepared from isolated ovarian samples. Therefore, the following analysis will only deal with whole flies or carcasses. We previously analyzed RNA-seq data produced from Sindbis virus (SINV) infection in the same Makindu D. simulans strain and in the w control D. melanogaster strain (Roy et al. 2020). These experiments revealed that TE transcript amounts decreased in somatic tissues upon SINV infection in both strains. We also added data that correspond to DCV infection in D. melanogaster y strain (Merkling et al. 2015), which induces an increase in TE transcript amounts (Roy et al. 2020). In addition, we included data that correspond to Kallithea virus (KV) infection in D. melanogaster OregonR (Palmer, Medd, et al. 2018). We only kept data corresponding to females. In this last data set, our TE analysis procedure revealed an average decrease of TE transcript amounts (14% decrease, mean log2FC = −0.21, Wilcoxon paired test P value = 9.3e-9, supplementary fig. S7, Supplementary Material online).

Table 1

Characteristics of the Data Sets under Study

Data Set	Present Study	Roy et al. (2020)	Roy et al. (2020)	Merkling et al. (2015)	Palmer et al. (2018)
Mean TE log2FC	0.31	−0.25	−0.32	0.28	−0.21
Strain	Makindu	Makindu	w1118	y1	OregonR
Species	D. simulans	D. simulans	D. melanogaster	D. melanogaster	D. melanogaster
Virus	DMelSV	SINV	SINV	DCV	KV
Viral genome	ssRNA(−)	ssRNA(+)	ssRNA(+)	ssRNA(+)	dsDNA
Natural host	D. melanogaster	Mosquito	Mosquito	D. melanogaster	D. melanogaster
Infection route	Intrathoracic injections	Intrathoracic injections	Intrathoracic injections	Intrathoracic injections	Abdominal injections
Dpi	8	6	6	1	3
Tissue	Carcasses	Carcasses	Carcasses	Whole flies	Whole flies
Sex	Females	Females	Females	Females	Females

We are aware that experimental conditions vary across these studies. Therefore, rather than focusing on absolute transcript levels in each condition, we only considered log2-fold changes, which reflect changes induced by infection, regardless of baseline levels. We compiled these data and analyzed genes involved in immune pathways performing a principal component analysis (PCA) of log2-fold changes between infected and mock conditions for each experimental system (fig. 4 and supplementary fig. 8, Supplementary Material online). The first axis corresponded to 35.2% of the total variance and clearly separated {Makindu SINV} and {w SINV} from one hand apart from {Makindu DMelSV} and {y DCV} from the other hand. Therefore, this first axis segregates RNA virus samples that displayed TE transcript decrease from one hand and TE transcript increase from the other hand. This suggests that the transcriptional response of immune genes upon an RNA virus infection may allow to predict the effect on TE transcript modulation. However, this does not hold true for KV, the only DNA virus of the analysis. The second axis of the PCA explained 25.3% of the total variance. The first two axes clearly separated genes involved in RNAi from the other pathways and were associated with {Makindu SINV} sample.

Fig. 4

PCA of immune gene transcriptional changes upon viral infections. PCA was performed on log2FC values. The first and second axes explain 35.2% and 25.3% of the total variance, respectively. Samples are labeled according to the fly strain and the virus used. Samples displaying an increase in TE transcript amounts upon infection are tagged in yellow; otherwise, they are tagged in gray.

Discussion

DMelSV Replicates in D. simulans, Triggers an RNAi Response, and Activates the Toll and Jak/STAT Pathways

Although DMelSV is frequently found in D. melanogaster samples, it had never been observed in D. simulans natural populations (Webster et al. 2015). However, our results show that this virus is able to replicate in D. simulans and to induce the same symptoms following CO2 anesthesia. This is consistent with the findings of Longdon et al. (2011), who previously measured the ability of sigma viruses to replicate in various dipteran hosts following injection, and found that DMelSV could also replicate in D. simulans despite this species not being its natural host. We were able to detect DMelSV RNAs in ovaries of flies infected by thoracic injections. In addition, we found DMelSV-derived 21-nt small RNAs, equally produced in both orientations, indicative of an induced, efficient RNAi response. However, our RNA-seq data revealed that other immune pathways were also triggered. Several genes of the Toll and Jak/STAT pathways displayed increased transcript amounts upon infection, suggesting that these pathways play a role in the host immune response. Using microarrays, a previous study could not detect any evidence for activation of the Toll pathway upon DMelSV infection (Carpenter et al. 2009), whereas another study using RT-qPCR reported nonstatistically significant increases in the expression of Toll, Relish, and vir-1 (Tsai et al. 2008). The greater sensitivity of the RNA-seq approach may explain why we were able to detect weak increases in transcript amounts for these genes in the present study. In addition, these differences may also result from differences due to the species: the above studies were performed in D. melanogaster whereas ours used D. simulans.

As described in our previous study, we observed that Makindu flies were chronically infected with NoraV. This chronic infection was proposed not to impact the outcome of experimental acute infections (Roy et al. 2020). In addition, we found a strong signal for DAV small RNAs, whereas DAV RNAs were never observed. Using small RNA-seq and RT-PCR, we did not detect DAV small RNAs and RNAs from the DMelSV viral stock. The absence of detectable DAV RNA in the infected flies may be related to the time course of DAV replication, its context with the siRNA pathway, or specifically in the Makindu strain. Although we cannot completely rule out the possibility of a double infection with DMelSV and DAV, the relatively frequent occurrence of DAV small RNAs in Drosophila samples (Palmer, Medd, et al. 2018; Roy et al. 2020) allows us to assume that the potential impacts on the outcome of DMelSV infection are not significant.

TE Transcript Amounts Are Protected against Modulation in Ovaries

Using SINV in Makindu and w strains, we previously found that viral infections did not induce modulations in ovarian TE transcript amounts (Roy et al. 2020). However, we could observe that SINV did not enter ovaries, which could explain the absence of effect on ovarian TE transcripts. Here we chose to use DMelSV, which is known to be vertically transmitted through eggs (L’Héritier 1958), and the RNA of which was detected in ovarian samples in the present study. Nevertheless, we cannot exclude that viral replication is less intense in ovaries compared with somatic tissues. Surprisingly, in Makindu ovaries, DMelSV infection led to virtually no changes in TE transcript amounts, whereas TE-derived small RNA amounts increased. This result suggests that there is a strong pressure to maintain ovarian TE transcript homeostasis. Such a pattern is probably the consequence of strong selective pressures on the maintenance of germline cells' genomic integrity. Nevertheless, we have to note that our previous experiments using SINV in Makindu revealed that ovarian TE transcript homeostasis was associated with an increase in TE-derived siRNAs but a strong decrease in TE-derived piRNAs. On the contrary, no modulations of ovarian TE-derived small RNA repertoires were observed upon SINV infection in D. melanogaster w. These observations suggest that the dedicated pathways to TE homeostasis rely on TE-derived small RNA repertoire adjustments, and vary across host species, but the exact mechanisms remain to be established.

siRNA pathway genes can be thought of as being engaged in a Red Queen arms race due to the antagonistic evolutionary interactions between viruses and the host genomes (Obbard et al. 2006; Kolaczkowski et al. 2011). Similarly, piRNA pathway genes can be thought of as being engaged in a Red Queen relationship due to the antagonistic evolutionary interactions between TEs and the host genomes (Parhad et al. 2017; Palmer, Hadfield, et al. 2018) and the potential deleterious off-targets of genomic auto-immunity (Blumenstiel et al. 2016; Wang et al. 2020). Considering that these selective pressures lead to rapid, and a priori lineage independent, evolution of all partners, it is not unexpected to observe contrasted patterns between D. melanogaster and D. simulans. It is all the more expected in the ovaries, considering that it is the supposed place of the strongest interaction between the roles of both siRNA and piRNA pathways. This may add an additional layer of complex interactions acting on these different small RNA pathways. This reinforces the need for including inter and intraspecific variability in the study of such complex and rapid evolutionary processes.

In any case, the strong homeostasis of TE transcript amounts in ovaries suggests an optimal transposition rate in the oocyte (Dufourt et al. 2014), which we may imagine ensures an optimal mutational environment for the oocyte genome. This transcript homeostasis may rely on strong modulations of ovarian TE-derived small RNA amounts, which would lead to modulations of the TE-derived small RNA repertoires that are transmitted to the zygote. In the case of SINV, we observed a decrease in TE-derived small RNA amounts, which may lead to an increased mutation rate due to a release of TE activity. On the contrary, in the present study using DMelSV, we observed an increase in TE-derived small RNA amounts, which may ensure a greater genomic stability as a consequence of a stronger control of TEs in the zygote. For instance, we may envision that such a modulation reduces the instability resulting from hybrid crosses (Chambeyron et al. 2008; Akkouche et al. 2013). In any case, these contrasted patterns suggest that the study of a larger set of viruses and Drosophila strains will be needed to understand these aspects into more details (Palmer, Varghese, et al. 2018).

The Somatic Fate of TEs upon Viral Infections Depends on the Triggered Immune Pathways

The present study of DMelSV infection in D. simulans Makindu revealed that somatic TE transcript amounts were modulated upon infection, in agreement with our previous results using SINV (Roy et al. 2020). However, contrary to these previous observations, here we found an increase in somatic TE transcript amounts, which was associated with a decrease in TE-derived small RNAs. This suggests that somatic transposition rates may increase upon DMelSV infection. This is particularly interesting because somatic transposition was proposed to be involved in various physiological processes such as in brain functions and aging, especially in mammals (Faulkner and Garcia-Perez 2017) but also in Drosophila (Chang et al. 2019).

Although we cannot exclude that the abundant production of viral siRNAs limits the production of TE-derived siRNAs due to titration mechanisms and induces a global decrease in relative TE-derived siRNA amounts, the present results confirm that viral infections impact somatic TE transcript amounts in Drosophila, and suggest that the outcome of the infection on TEs likely depends on the considered virus and most probably on the main immune pathway that is triggered. In order to go deeper into the underlying mechanisms, we performed a meta-analysis of RNA-seq data produced upon infection using different viruses in different strains of D. melanogaster and D. simulans. We focused on a subset of genes known to be implicated in the various immune pathways characterized in flies.

It is true that many experimental conditions vary across samples (table 1). For instance, due to experimental constraints, samples were not produced at the same time (SINV: 6 dpi; DMelSV: 8 dpi; DCV: 1 dpi; KV: 3 dpi). In addition, the tissues are not exactly the same: In the DCV study, whole flies are considered, whereas we separated ovaries from the rest of the body (carcasses). Moreover, the viruses are different in their being natural pathogens of D. simulans or even of flies. They also differ regarding the organization of their genomes: single-stranded RNA genome of negative polarity for DMelSV, single-stranded RNA genome of positive polarity for SINV and DCV, and DNA genome for KV. However, despite all these differences, the meta-analysis allows to suggest future directions of research.

Our PCA analysis on immune genes log2FC displays a clear structuration of RNA virus samples according to the direction of the modulation of TE transcript amounts. The sample corresponding to KV, the only DNA virus from the analysis, does not follow this trend. This suggests that, although KV infection seems to trigger an RNAi response (Palmer et al. 2018), other mechanisms should be considered regarding the impacts on TE modulation. In the case of RNA viruses, what is particularly interesting is that the RNAi pathway appears to determine this dichotomy. This suggests that when RNAi is the first line of defense against the virus—this is particularly true for nonnatural pathogens, such as SINV in Drosophila—this immune response leads to an increase in TE-derived small RNAs and therefore a decrease in TE transcript amounts, in agreement with the model we previously proposed (Roy et al. 2020). On the contrary, when other immune pathways are on the front line, such as the Toll and Jak/STAT pathways here in the case of DMelSV, or when the virus encodes an RNAi suppressor (such as DCV), the host immune response leads to a global decrease in TE-derived small RNAs, and consequently an increase in TE transcript amounts. Very interestingly, DMelSV infections trigger the Toll pathway and induce TE transcript increase, whereas KV encodes a Toll pathway suppressor and leads to TE transcript decrease (Palmer et al. 2019). More work is needed to explore these avenues.

Conclusion

Based on the present study, we propose that ovarian TE transcript amounts are protected against modulation upon viral infection, even in the case of viruses that display a clear tropism to ovaries, such as DMelSV. The maintenance of ovarian TE transcript homeostasis is associated with significant modulations of TE-derived small RNAs. Determining the precise mechanisms will require further investigation.

In addition, using DMelSV this study also allows to confirm that somatic TE transcript amounts are modulated upon viral infection. If this was to lead to increased transposition rates, this would reinforce the involvement of viral infections in somatic genetic diversification. In addition, the analysis we perform here suggests that the direction of somatic TE transcript modulation depends on the immune pathways that are triggered (Palmer et al. 2018). This calls for the need of enlarging the panel of viruses under study in order to investigate the corresponding mechanisms.

Materials and Methods

Fly Strains and Husbandry

We used the Makindu D. simulans wild-type strain. Flies were reared on a standard corn medium at 25 °C. To eliminate the Wolbachia endosymbiotic bacteria, which is known to interfere with viral replication (Teixeira et al. 2008; Bhattacharya et al. 2017), flies were reared during two generations on standard agarose medium containing 0.25 mg/ml tetracycline hypochloride (Sigma-Aldrich). The absence of Wolbachia was confirmed using Wolbachia 16S primers (5′-TTG TAG CCT GCT ATG GTA TAA CT-3′ and 5′-GAA TAG GTA TGA TTT TCA TGT-3′), and Wolbachia wsp primers (5′-TGG TCC AAT AAG TGA TGA AGA AAC-3′ and 5′-AAA AAT TAA ACG CTA CTC CA-3′), as described previously (Riegler et al. 2005).

Virus Production and Titration

DMelSV stock was harvested from D. melanogaster chronically infected by Hap23 sigma virus (provided by Francis Jiggins, Cambridge University) (Carpenter et al. 2007). Briefly, DMelSV was extracted by grinding 100 infected fruit flies in PBS supplemented with 10% fetal calf serum (FCS). The supernatant was centrifuged 2 min at 10,000 × g and filtered at 0.45 µm. CO2 test was used to titrate DMelSV stock (Longdon et al. 2012), using 10-fold serial dilutions of stocks and intrathoracic injections (Nanoject II apparatus; Drummond Scientific) of naive D. melanogaster flies. At 14 dpi, infected flies were exposed to 100% CO2 and flies were subsequently deposited on ice for 15 min. At room temperature, surviving flies were counted 24 h after the CO2 test and stock viral titer was valued using the Reed and Muench’s method and expressed as log10 50% Drosophila infective doses (DID50) (Longdon et al. 2012).

Fly Infections and Validation of the Viral Infection

Flies were individually intrathoracically injected (Nanoject II apparatus; Drummond Scientific) with 207 nl containing 81 DMelSV infectious particles per fly (DID50) in PBS, 10% FCS. In parallel, control flies were mock infected with PBS, 10% FCS under the same conditions. Fly mortality at 1 dpi was attributed to damage due to the injection procedure, and these flies were excluded from further analyses. Flies were maintained at 25 °C and transferred to fresh medium every 2–3 days. Viral infection of individual flies (n = 10 per experiment) was validated using CO2 test. DMelSV injection experiments were performed independently, in duplicates.

RNA Extraction and RT-PCR

Total RNA from eight pairs of ovaries (biological duplicates) were extracted at 8, 10, 12, and 14 dpi using Trizol (Thermo Fisher Scientific) and RNeasy Mini Kit (Qiagen). Purified RNAs were treated using Turbo DNase (Ambion DNAfree kit). Reverse transcription was performed on 500 ng of total RNAs using Omniscript Reverse Transcriptase (Qiagen). Endpoint PCR was performed using the GoTaq G2 DNA Polymerase kit (Promega). The parameters for PCR cycling were as follows: one cycle of 94 °C for 2 min; 35 cycles of 94 °C for 30 s, 57 °C for 30 s and 72 °C for 30 s, followed by a terminal extension of 72 °C for 5 min. The primers used for DMelSV PCR amplify a fragment of 156 bp (DMelSV: 5′-ATG TAA CTC GGG TGT GAC AG-3′ and 5′-CCT TCG TTC ATC CTC CTG AG-3′). The positive control corresponds to the DMelSV stock. rp49 was used as RNA positive control (rp49: 5′-CGG ATC GAT ATG CTA AGC TGT-3′ and 5′-GCG CTT GTT CGA TCC GTA-3′) (Akkouche et al. 2012).

RNA-Seq

Female flies were carefully, manually dissected: pairs of ovaries were separated from the rest of the bodies, which will further be called “carcasses.” Ovaries had to remain intact for the samples to be kept. RNA extraction was performed at 8 dpi, which corresponds to the exponential phase of viral replication. Total RNA of 30 pairs of ovaries or 30 carcasses (biological duplicates) were extracted using Trizol (Thermo Fisher Scientific) coupled with the RNeasy Mini Kit (Qiagen). Purified RNAs (10 µg) were treated using TurboDNAse (Ambion DNAfree kit). RNA quality was validated using a Bioanalyzer (Agilent). Libraries were prepared after polyadenylated RNA selection. Purified RNAs were used for library preparation using the TruSeq RNA Sample Prep Kit v2. Enriched libraries were paired-end sequenced on a Hiseq 4000 apparatus (Illumina). Sequencing adaptors were removed using UrQt (Modolo and Lerat 2015) and trim Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Trimmed reads were aligned on reference genes retrieved from FlyBase (ftp://ftp.flybase.net/releases/FB2019_01/dsim_r2.02/fasta/dsim-all-gene-r2.02.fasta.gz) using TopHat2 (Kim et al. 2013). Gene count tables were generated using eXpress (Roberts et al. 2011). TE count tables were generated using the TEcount module of TEtools (Lerat et al. 2017), and the list of sequences used in Roy et al. (2020) and available at ftp://pbil.univ-lyon1.fr/pub/datasets/Roy2019/. Gene and TE count tables were concatenated to make the complete count table, which was further analyzed using the DESeq2 R package (version 1.18) (Love et al. 2014).

Tables compiling raw and normalized read counts for TE families are provided as supplementary material, Supplementary Material online, as well as DESeq2 results corresponding to TE families and immune genes.

Small RNA-Seq

Female flies were carefully, manually dissected: pairs of ovaries were separated from the rest of the bodies, which will further be called “carcasses.” Ovaries had to remain intact for the samples to be kept. Small RNA extraction was performed at 8 dpi, which corresponds to the exponential phase of viral replication. Small RNAs were isolated from 50 pairs of ovaries or 50 carcasses using a HiTrap Q HP anion exchange chromatography column (GE Healthcare Life Sciences), as described previously (Grentzinger et al. 2014). Pairs of ovaries or carcasses were crushed on HQ 0.1 solution (Hepes 20 mM, KOH 20 mM, EDTA 0.2 mM, MgCl2 1.5 mM, glycerol autoclaved 10%, Ack 0.1 M, Pefabloc SC 0.1 mg/ml DTT 1 mM) and the supernantant was added on the HiTrap column, followed by a high Ack concentration solution (HQ 4 M). Small RNAs were purified using phenol/chloroform isoamyl alcohol extraction (Sigma-Aldrich). The same experimental procedure was performed on two replicates of 150 µl of DMelSV viral stock, which corresponded approximately to 30 Hap23 infected flies each. Size selection of small RNAs (18–50 bp) was performed on gel at the GenomEast platform. Purified small RNAs were used for library preparation (SQ00/SIL-04-SR) and sequenced on a Hiseq 4000 apparatus (Illumina). We first removed adapter sequences using cutadatp (Martin 2011) -a TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACTTA -m 1. At this step, we drew the size distribution of read length, which almost exclusively ranged between 21 and 30 nt, that is, the expected size for siRNAs, piRNAs, and miRNAs; this made us confident that we did not extract RNA degradation products (supplementary fig. S3, Supplementary Material online). Using PRINSEQ lite version 0.20.4 (Schmieder and Edwards 2011), we filtered reads of size 23 to 30 nt and considered these as piRNAs, and reads of size 21 nt and considered these as siRNAs. Per base sequence contents were then computed using FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc/) and their analysis revealed the expected 1U enrichment (supplementary fig. S4, Supplementary Material online). Using a modified version of the TEcount module of TEtools (Lerat et al. 2017; Fablet et al. 2019), we mapped the reads against a list of TE sequences (available at ftp://pbil.univ-lyon1.fr/pub/datasets/Roy2019/) and virus sequences (sequences from Obbard laboratory, available at http://obbard.bio.ed.ac.uk/data/Updated_Drosophila_Viruses.fas.gz), either in the sense or antisense directions. We then computed the sums of sense and antisense counts. We deepened the description of sequences corresponding to these small RNA repertoires and counted the number of reads mapping against a multifasta corresponding to the concatenated following files from FlyBase: dsim-all-gene-r2.02 (after masking using Repeat Masker), dmel-all-miscRNA-r6.16.fasta, dmel-all-tRNA-r6.16.fasta, and dsim-all-miRNA-r2.02.fasta. Alignments were performed using bowtie --best (Langmead 2010).

Size distribution profiles of virus-mapping reads were performed using the following reference sequences for DMelSV, DAV, and NoraV: Genbank GQ375258, NC_012958.1, and NC_007919.3 accession numbers, respectively. Alignments were performed using bowtie --best (Langmead 2010). The mapping of DMelSV-derived small RNAs was performed using the SAMtools suite (Li et al. 2009) and the bamCoverage tool from deepTools2 (Ramírez et al. 2016).

Read count numbers were normalized according to the numbers of 19–39 nt reads aligning on miRNA sequences (miRNA sequences retrieved from FlyBase: dsim-all-miRNA-r2.02.fasta.gz). Tables compiling raw and normalized read counts for TE families are provided as supplementary material, Supplementary Material online.

Statistical Analyses

We tested whether the modulations of TE counts upon infection were different from zero using Wilcoxon paired tests comparing TE normalized counts in the infected versus mock conditions, for the considered 237 TE families. The null hypothesis corresponds to an equal number of dots above and below the 0 horizontal line. Differential expression of genes was tested using the DESeq2 R package (Love et al. 2014), and significance was assessed based on adjusted P values lower than 0.05.

Meta-analysis of RNA-Seq Data

We retrieved previously published RNA-seq data accessible under numbers PRJNA540249 (Roy et al. 2020), GSE56013 (Merkling et al. 2015), and PRJEB21366 (Palmer, Medd, et al. 2018). Roy et al.’s and Merkling et al.’s data were processed as described previously (Roy et al. 2020). Regarding Palmer et al.’s data, we first removed adapters using Trimmomatic (Bolger et al. 2014). To build the gene count table, we aligned reads against D. melanogaster reference genes (ftp://ftp.flybase.net/releases/FB2020_04/dmel_r6.35/fasta/dmel-all-gene-r6.35.fasta.gz) using HISAT2 (Kim et al. 2019), SAMtools (Li et al. 2009), and eXpress (Roberts et al. 2011). The TE count table was obtained following the same procedure as above. Gene and TE count tables were then concatenated to perform all subsequent analyses.

For each data set independently, we computed normalized counts using DESeq2 (Love et al. 2014), and manually calculated log2 fold changes between mean normalized counts in infected condition and mean normalized counts in mock condition, for genes involved in immune pathways (see supplementary table S6, Supplementary Material online for details). All log2FC values were then gathered in the same table, which was used to perform a PCA (ade4 R package, Dray and Dufour 2007).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.