Effects ofWolbachiaon Transposable Element Expression Vary BetweenDrosophila melanogasterHost Genotypes

GBE Effects of Wolbachia on Transposable Element Expression Vary Between Drosophila melanogaster Host Genotypes Ana T. Eugénio 1 , Marta S.P. Marialva 1 , and Patrícia Beldade 1,2, * 1 Instituto Gulbenkian de Ciência, Oeiras, Portugal 2 cE3c (Center for Ecology, Evolution and Environmental Changes) and CHANGE (Global Change and Sustainability Institute), Faculty of Sciences, University of Lisbon, Lisbon, Portugal *Corresponding author: E-mail: [email protected]. Accepted: 19 February 2023 Abstract Transposable elements (TEs) are repetitive DNA sequences capable of changing position in host genomes, thereby causing mutations. TE insertions typically have deleterious effects but they can also be beneficial. Increasing evidence of the contribution of TEs to adaptive evolution further raises interest in understanding what factors impact TE activity. Based on previous studies associating the bacterial endosymbiont Wolbachia with changes in the abundance of piRNAs, a mechanism for TE repression, and to transposition of specific TEs, we hypothesized that Wolbachia infection would interfere with TE activity. We tested this hypothesis by studying the expression of 14 TEs in a panel of 25 Drosophila melanogaster host genotypes, naturally infected with Wolbachia and annotated for TE insertions. The host genotypes differed significantly in Wolbachia titers inside individual flies, with broad-sense heritability around 20%, and in the number of TE insertions, which depended greatly on TE identity. By removing Wolbachia from the target host genotypes, we generated a panel of 25 pairs of Wolbachia-positive and Wolbachia-negative lines in which we quantified transcription levels for our target TEs. We found variation in TE expression that was dependent on Wolbachia status, TE identity, and host genotype. Comparing between pairs of Wolbachia-positive and Wolbachia-negative flies, we found that Wolbachia removal affected TE expression in 21.1% of the TE-genotype combinations tested, with up to 2.3 times differences in the median level of transcript. Our data show that Wolbachia can impact TE activity in host genomes, underscoring the importance this endosymbiont can have in the generation of genetic novelty in hosts. Key words: transposable element activity, Wolbachia infection, Drosophila melanogaster, DGRP, genetic variation. Significance Mobilization of transposable elements (TEs) generates mutations that can contribute to adaptive evolution, making it all the more relevant to understand what factors affect TE activity. We show that infection with a common endosymbiotic bacterium, Wolbachia, affects TE activity in Drosophila melanogaster hosts, in a manner that varies depending on TE identity and on host genotype. © The Author(s) 2023. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 1 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023 Eugénio et al. GBE Introduction Transposable elements (TEs) are repetitive DNA sequences capable of changing position independently in the host genome (Bourque et al. 2018; Mérel et al. 2020), and make up a significant fraction of many eukaryotic genomes (Guio and González 2019). They are divided into two major classes, depending on whether the mechanism of transposition does (for retrotransposons) or does not (for DNA transposons) involve an RNA intermediate that is reverse transcribed before integrating back into the host genome (Bourque et al. 2018). TE insertions can cause mutations, which typically have deleterious effects because they disrupt proper gene function in a variety of manners (McFaddenaf and Knowlesb 1997; Hedges and Deininger 2007; Belancio et al. 2008; Ayarpadikannan and Kim 2014). Consequently, host organisms have evolved mechanisms to control and repress TE activity, including the piRNA pathway in animals (Tóth et al. 2016). On the other hand, an increasing number of studies have been providing compelling examples of TE insertions with positive effects on host fitness, contributing to adaptation (González and Petrov 2009; González et al. 2010; van’t Hof et al. 2016), stress resistance (Guio et al. 2014; Pereira and Ryan 2019), and the origin of novel traits (Emera and Wagner 2012; Bennetzen and Wang 2014; Santos et al. 2014; Trizzino et al. 2017). Moreover, TEs might also contribute to reproductive isolation, as in the case of TE-mediated hybrid incompatibility (Petrov et al. 1995; Serrato-Capuchina and Matute 2018). TE contribution to adaptive evolution and diversification raises interest in understanding what factors impact TE activity. TE activity differs between TEs (Venner et al. 2009; Mérel et al. 2020) and between host genotypes (Anderson et al. 2019; Signor 2020; Wang et al. 2022). Furthermore, studies on different organisms have shown that TE activity can be affected by external environmental factors, including temperature (Chen et al. 2018), radiation (Newman et al. 2014), heavy metals (Habibi et al. 2014), starvation (Rep et al. 2005), and various other stressors (Miousse et al. 2015). Not much is known about how these factors can affect the molecular mechanisms responsible for TE regulation, including the piRNA pathway. On the other hand, Wolbachia, a common endosymbiotic bacterium, has been shown to affect the abundance of some piRNAs (in Aedes aegypti mosquitoes; Mayoral et al. 2014) and the rate of transposition of the retrotransposon gypsy (in Drosophila melanogaster; Touret et al. 2014). Moreover, the invasion of the DNA transposon P-element in populations of Drosophila reportedly co-occurred with a replacement of Wolbachia strain infecting those flies (Riegler et al. 2005). However, there has been no systematic analysis of the effects of Wolbachia on the activity of different TEs in different host genotypes. Wolbachia is a maternally inherited endosymbiont that is prevalent in invertebrates, including insects, arachnids, and nematodes (Werren et al. 2008; Kaur et al. 2021). Multiple studies have documented Wolbachia prevalence (Clark et al. 2005; Riegler et al. 2005; Weeks et al. 2007) and load (López-Madrigal and Duarte 2019; Liu and Li 2021) in natural and laboratory populations of Drosophila hosts. Associated with its mode of transmission, Wolbachia can have important effects on host reproduction, being responsible for phenomena such as cytoplasmic incompatibility, feminization, and male killing (Werren et al. 2008; Kaur et al. 2021). Wolbachia can also affect other aspects of host biology, including resistance to viral infection (Teixeira et al. 2008), gut microbiome composition (Simhadri et al. 2017), thermal preference (Truitt et al. 2019), sleep behavior (Bi et al. 2018), and fecundity and lifespan (Serga et al. 2021). At the molecular level, Wolbachia is known to affect host gene expression (Baião et al. 2019; Biwot et al. 2020), and meiotic recombination rate (Singh 2019), as well as the aforementioned TE-related properties (Riegler et al. 2005; Mayoral et al. 2014; Touret et al. 2014). Here, we test the impact of Wolbachia on TE expression by using host lines where Wolbachia is present versus where it was removed. Specifically, we use flies from the Drosophila melanogaster Genetic Reference Panel (DGRP), a panel of isogenic lines derived from a natural population, whose genomes have been fully sequenced and annotated for TE insertions (Mackay et al. 2012; Rahman et al. 2015). We selected 25 DGRP lines that were naturally infected with Wolbachia for which we estimated Wolbachia loads in individual flies and recorded the number of TE insertions for 14 TEs, representing different families. We found differences in Wolbachia loads and in number of TE insertions between genotypes, as well as an association between the two. We then generated a Wolbachia-free counterpart for each of the 25 target genotypes and used our panel of 25 paired Wolbachia-positive and Wolbachia-negative lines to quantify transcription levels of the 14 target TEs. We found variation in TE expression depending on host genotype, TE identity, and Wolbachia status. Whether Wolbachia removal led to increased or decreased TE expression appeared to be more of a property of host genotype than of TE identity. Results and Discussion To investigate the effect of Wolbachia infection on TE expression, we focused on 25 D. melanogaster genotypes, for which we documented differences in Wolbachia loads and in number of insertions of 14 target TEs (fig. 1). We focused specifically on what were called “novel insertions” (Mackay et al. 2012; Rahman et al 2015), which correspond 2Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023 Effects of Wolbachia on Transposable Element Expression GBE to TE insertions found in the DGRP genomes but not in Release 6 of the D. melanogaster’s reference genome (counts cf. Rahman et al. 2015). We then generated a corresponding panel of 25 lines from which Wolbachia was cleared, and compared expression level of our target TEs in adult females between the pairs of Wolbachia-positive (Wolb+) and Wolbachia-negative (Wolb−) flies (fig. 2). Host Genotypes Differ in Wolbachia Loads and in Number of TE Insertions We randomly chose 25 of the 85 DGRP lines known to be infected with the wMel strain of Wolbachia (Mackay et al. 2012; Richardson et al. 2012). For each of these lines, we measured Wolbachia loads in five individual adult females 10 days post-eclosion, the same sex and age used to measure TE expression. For this, we used quantitative real-time polymerase chain reaction (qPCR) with primers for one Wolbachia-specific gene (wsp), to estimate number of bacterial cells, and for one host-specific gene (actin), to assess number of host cells. Across the ∼125 flies assayed individually, Wolbachia loads varied between a minimum of 3.5 and a maximum of 51 Wolbachia cells per host cell. Only six individuals, of different genotypes, had >20 Wolbachia per host cell. These estimates of Wolbachia density fall along the same order of magnitude as those found through sequencing of the DGRP lines (0.9–17.1 copies per host cell; Richardson et al. 2012), or through qPCR of whole bodies (Bénard et al. 2021; Chrostek et al. 2021) and gonadal tissues (Correa and Ballard 2012) of other D. melanogaster genotypes, as well as for other Wolbachia strains (Chrostek and Teixeira 2015). We found differences in Wolbachia loads between host genotypes, with median values ranging from 5 to 15 copies of Wolbachia per host cell (fig. 1A), and estimated FIG. 1.—Characterization of our 25 target host lines in relation to Wolbachia load in individual flies (A) and to the number of novel TE insertions in their genomes (B). The 25 genotypes are organized along the y-axis in order of the median value of Wolbachia load. (A) Wolbachia load relative to number of host cells (x-axis). Each blue dot is a biological replicate and represents one single female. Wolbachia load is significantly different across genotypes (ANOVA; F 24,5508 = 4.5e + 27, P < 2e−16). (B) Heatmap representing the predicted number of novel insertions for our 14 target TEs. TEs are organized in the x-axis by median number of novel insertions across genotypes. The scale of gray, from white to dark gray, represents, respectively, from the lowest to highest number of novel insertions. The TE pogo in genotype RAL-21 is out of the scale, with 216 novel insertions annotated. There was no information in Rahman et al. (2015) for RAL-855. The number of novel insertions differed significantly between genotypes (ANCOVA; genotype: F 22,4974 = 1.3e + 28, P < 2e−16, TEs: F 13,4974 = 1.7e + 29, P < 2e−16, and Wolbachia titers: F 1,4974 = 2.2e + 27, P < 2e−16). Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 3 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023 Eugénio et al. GBE FIG. 2.—Expression levels of 14 TEs in adult female flies of 25 genotypes with versus without Wolbachia. TEs are ordered from left to right by median (and average, for tied medians) number of novel insertions. Statistical significance for expression differences between Wolb+ and Wolb− is shown as * for P < 0.05, ** for P < 0.01, and *** for P < 0.001 (ANOVA, see Material and Methods). (A) Expression of the 14 TEs in genotypes RAL-21, RAL-181, RAL-712, and RAL-737. The same plots for all other genotypes can be found in supplementary figure S2, Supplementary Material online. (B) Expression of different TEs across various genotypes, illustrating cases where expression levels are statistically different between Wolbachia status. Each dot in plots (A) and (B) represents a biological replicate, corresponding to a pool of ten female flies. (C) Heatmap representing differences in expression level for the 14 target TEs between Wolb + and Wolb− flies of all 25 different genotypes. Genotypes in the y-axis are ordered by Wolbachia load (as in fig. 1B). Cells are displayed in a gradient of color, representing effect size (color intensity) and whether expression is higher in Wolb+ relative to Wolb− (blue shades; top half of the gradient legend right to the heatmap) or the other way around (pink shades; bottom half of the gradient legend right to the heatmap). Underlined asterisks represent significant differences after Benjamini–Hochberg correction for multiple comparisons. (D) Volcano plot representing the effect size (x-axis) and P-value (testing for log2 foldchange of TE expression differences between Wolb+ and Wolb−; y-axis). Dots relative to largest effect sizes (blastopia in RAL-321 and blood in RAL-595), and to nonsignificant effect size below −1, corresponding to 100% difference in TE expression after Wolbachia removal (blastopia in RAL-440 and Quasimodo in RAL-440), are labelled. The dashed grey line represents the threshold of statistical significance for differences in TE expression levels (P = 0.05). Effect sizes in (C) and (D) correspond to the difference of the log2 fold-change in TE expression between Wolb+ and Wolb− paired genotypes. 4Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023 Effects of Wolbachia on Transposable Element Expression GBE broad-sense heritability (H 2 ) of 0.22 (among-line variance = 10.3, within-line variance = 36.0). Although little is known about what host loci harbor natural allelic variation contributing to variation in Wolbachia loads, we do know that loads vary with environmental factors, including temperature (Wiwatanaratanabutr and Kittayapong 2009), host diet (Ponton et al. 2015; Serbus et al. 2015), and viral infection (Kaur et al. 2020). First, we validated in silico predictions of insertions (Mackay et al. 2012) by PCR with primers for the fly genomic sequence flanking 132 predicted novel insertions in 11 genotypes (supplementary table S1, Supplementary Material online). The amplicons from each of the insertions were sized (agarose gel) and sequenced to confirm the presence, length, and identity of the inserted DNA (supplementary table S1, Supplementary Material online). For 100% of the predicted insertion locations we tested, we confirmed the presence of a TE insertion, and, in most cases, we also confirmed that the size and the sequence of the inserted DNA corresponded to the predicted TE identity (supplementary fig. S1A, Supplementary Material online). For 113 (85.6%) of the insertions tested, the inserted TE corresponded to the most likely expected identity (cf. the predictions made from the whole-genome sequence data), and for 16, it corresponded to the second most likely TE (Mackay et al. 2012). We observed that 67 insertions (50.8%) had the size corresponding to the expected full length of that TE, 44 (33.3%) were smaller, and 21 (15.9%) were larger (supplementary fig. S1A and table S1, Supplementary Material online). With predictions of novel insertions validated, we used data from the TIDAL-FLY v1.0 tool of Rahman et al. (2015) to gather information about the number of novel insertions for each of our 14 target TEs in 24 of our 25 study genotypes (there were no data for genotype RAL-855). We found significant differences in the number of novel insertions between genotypes and Wolbachia titers. The retrotransposons Cr1a, gypsy5, and Idefix had the lowest number of predicted novel insertions (with zero for the majority of the lines), whereas the DNA transposons 1360, pogo, and P-element generally had the highest number of novel insertions, in accordance with other studies describing DNA transposons as most active (Bourque et al. 2018). For most individual TEs, the estimated number of novel insertions varied between 0 and 44, with the exception of the TE pogo, predicted to have 216 novel insertions in the line RAL-21 (fig. 1B; Rahman et al. 2015). Note that the in silico predictions of the number of TE insertions are likely to be underestimates of the actual number of insertions. The DGRP lines were originally sequenced using a combination of Illumina and 454 sequencing technologies (Mackay et al. 2012), which generate short-reads and, as such, are not ideal for detecting TE insertions (Fiston-Lavier et al. 2015; Goerner-Potvin and Bourque 2018; Panda and Slotkin 2020; Rech et al. 2022). Moreover, new insertions may also have occurred after sequencing. We confirmed experimentally multiple “false negatives” in TE insertion predictions for the DGRPs. By running PCRs with TE-specific primers and DNA from eight DGRP genotypes predicted in Mackay et al. (2012) to have no insertions (novel or shared) of particular TEs. In all 17 cases tested, we verified the presence of those TEs (supplementary fig. S1B, Supplementary Material online). However, even if predictions are underestimates of actual number of insertions, the effects should be similar/random across TEs and genotypes with equivalent sequence coverage depth. TE Transcription Level Varies With Wolbachia Status in a Host Genotype-dependent Manner Using qPCR with TE-specific primers and a reference host gene, we quantified the expression of our 14 target TEs in eight replicate pools of ten 10-day-old females each, for each of the 25 Wolb+ and Wolb− pairs of genotypes (Cq data in supplementary table S2, Supplementary Material online). Expression levels differed significantly (analysis of variance, ANOVA), between TEs (F 13,4796 = 4.75, P = 2.97e−08), genotypes (F 24,4796 = 26.63, P < 2.2e−16), and with Wolbachia status (F 1,4796 = 28.79, P = 8.5e−08), with all interactions being significant (P < 0.0001 in all cases; fig. 2). Given the significant effect of Wolbachia status on TE expression, we then specifically compared expression of each of the 14 TEs in each of the 25 host genotypes with versus without Wolbachia. We found statistically significant differences in TE expression between Wolb+ and Wolb− lines for a total of 74 of the 350 (21.1%) genotype-TE combinations tested (fig. 2; supplementary fig. S2 and table S3, Supplementary Material online). We observed distinct scenarios depending on TE and genotype: higher expression in Wolb− flies for 47 in 350 cases (13.4%) and higher expression in Wolb+ flies for 26 in 350 cases (7.4%). For any given TE, the effect of Wolbachia on expression was not the same across genotypes, and, for any given genotype, the effect of Wolbachia on TE expression was not the same across TEs. However, for some genotypes, we observed some consistency in the effects of Wolbachia on TE expression. For genotypes RAL-142 and RAL-181, when statistically significantly different between Wolb+ and Wolb− flies, TE expression was always higher in Wolb+ (blue shades in fig. 2C) relative to Wolb− flies. Conversely, for genotypes RAL-712, RAL-21, and RAL-321, when significantly different, TE expression was always higher in Wolb− (pink shades in fig. 2C) relative to Wolb+ flies. Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 5 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023 Eugénio et al. GBE Effect Size of Wolbachia Removal on TE Expression Levels For each TE in each Wolb+/Wolb− genotype pair, the raw effect size of Wolbachia removal was calculated by subtracting the median log2 TE expression (normalized to reference gene) in the Wolb− flies from the median log2 TE expression in the Wolb+ flies (fig. 2C and D; supplementary table S5, Supplementary Material online). There were more cases in which Wolbachia removal resulted in an increase in TE expression than the reverse (i.e., more cases with significantly higher TE expression in Wolb− relative to Wolb+ flies; pink shades in fig. 2Cand pink dots in 2D). Also, the size of the effect of Wolbachia removal on TE expression tended to be larger when this removal resulted in increased expression relative to when it resulted in decreased expression (i.e., further from zero in the x-axis of fig. 2D). The median effect size for statistically significantly higher expression in Wolb− relative to Wolb+ (i.e., pink dots in fig. 2D) was approximately −1.0, corresponding to almost 100% increase in expression upon Wolbachia removal. On the other hand, the median effect size for statistically significant lower expression in Wolb− relative to Wolb+ (i.e., blue dots in fig. 2D) was around +0.5, corresponding to around 71% reduction in expression upon Wolbachia removal. In 24 of the 350 TE-genotype combinations tested, Wolbachia removal resulted in an increase in TE expression >100% increase (effect size <−1.0 in fig. 2D). Of these 24 cases, only two were not statistically significant (blastopia and Quasimodo, both in RAL-440; highlighted in fig. 2D). Interestingly, two genotypes, RAL-21 and RAL-712, stood out for having among the highest Wolbachia loads (fig. 1A) and a general increase in expression for most TEs after Wolbachia removal (fig. 2A and C), with median effect sizes across TEs of around −1.1 (107% increase) and around −0.9 (93% increase), respectively. In three of the 350 TE-genotype combinations tested, Wolbachia removal resulted in a >100% reduction (effect size >+1.0) in TE expression. Of these, blood in RAL-595 stood out with an effect size around 2.3, corresponding to a 246% reduction in expression after Wolbachia removal (fig. 2D). Conclusion In this study, we investigated the hypothesis of a relationship between Wolbachia infection and TE activity. Wolbachia is a prevalent endosymbiotic bacterium whose impact on TE mobilization was suggested by distinct lines of evidence, including: (1) effect on piRNA expression (Mayoral et al. 2014), (2) effect on rate of transposition of retrotransposon gypsy (Touret et al. 2014), and (3) Wolbachia strain-replacement co-occurring with invasion of DNA transposon P-element (Riegler et al. 2005). We tested whether the expression of 14 diverse TEs was different between Wolbachia-infected and Wolbachia-free D. melanogaster flies of 25 distinct genotypes differing in Wolbachia loads and in number of TE insertions. We focused on TE expression, which is often used as a proxy for TE activity (e.g., Becking et al. 2020; Torres et al. 2021), reasoning that higher expression creates more opportunities for insertions. TE transcript levels were quantified using qPCR, with effort put into carrying out and explaining in detail data structure, quality control, and analyses (see Materials and Methods). However, transcription is only one, albeit necessary, step in TE mobilization, and Wolbachia, or other factors, may impact TE integration post-transcriptionally. We found statistically significant differences in levels of TE transcript between Wolb+ and Wolb− flies in 21.1% of the 350 genotype-TE combinations analyzed, and a maximum effect size of 2.3 lower expression upon Wolbachia removal (supplementary table S5, Supplementary Material online). The observed effects of Wolbachia removal were not uniform for any given TE (i.e., one same TE could increase, decrease, or not change expression depending on genotype) nor for most genotypes (i.e., one same genotype could have TEs that increased, TEs that decreased, and TEs that did not change expression). However, some genotypes did stand out in having multiple TEs for which the direction Wolbachia effect on expression was the same. In particular, genotypes RAL-21 and RAL-712 showed both some of the highest Wolbachia loads and significant increase in expression for most TEs when Wolbachia was removed. Various factors can potentially lead to Wolbachia removal affecting TE expression. Wolbachia may affect piRNAs, as has been shown for a narrow set of piRNAs in A. aegypti mosquitoes (Mayoral et al. 2014), which are naturally devoid of Wolbachia. Wolbachia removal might also act as a stress factor for host genotypes, which might have adjusted to having Wolbachia. This study was performed using a subset of the DGRP lines, a panel of isogenic and fully sequenced genotypes that provide the possibility of looking at genotypic variation. Even though the genotypes are not naturally occurring, in that they were highly isogenized post-collection of a natural population, they represent naturally segregating allelic variants. Many studies showed differences between DGRP genotypes for various different types of traits (e.g., Magwire et al. 2012; Weber et al. 2012; Durlam et al. 2014; Ivanov et al. 2015; Howick and Lazzaro 2017; Lafuente et al. 2018; Mackay and Huang 2018). Our results highlight differences between genotypes in Wolbachia loads and number of TE insertions, as well as in the effect of Wolbachia removal on TE activity. The inter-genotype differences further emphasize the importance of analyzing multiple 6Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023 Effects of Wolbachia on Transposable Element Expression GBE genotypes to have a more complete understanding of any biological phenomena. Studies that only focus on a single or a few genotypes may miss or misrepresent general properties. Novel genetic variants created by TE mobilization can be, and are often, deleterious (McFaddenaf and Knowlesb 1997; Hedges and Deininger 2007; Belancio et al. 2008; Ayarpadikannan and Kim 2014). As such, high TE activity can put natural populations under stable conditions at a disadvantage. On the other hand, TE insertions can also be beneficial and, particularly in conditions of environmental perturbation, TE activity could contribute to novel genetic variants better adjusted to the changed conditions (e.g., Rey et al. 2016). The question of which and how intrinsic and extrinsic factors affect TE activity is a fundamentally interesting and largely unresolved question, especially for animal when comparing with plant TEs (Thieme et al. 2017). Our study shows that the maternally inherited Wolbachia endosymbiont, which is prevalent in insects and nematodes, affected TE expression in D. melanogaster. We expect future studies to provide insight about which and how different factors affect TE mobilization; including Wolbachia and other environmental factors in multiple hosts. Materials and Methods Confirming in silico Predictions of TE Insertions in DGRP Lines We looked to validate in silico predictions in terms of both potential false positives (focusing on specific insertions) and potential false negatives (focusing on particular TEs deemed as having no insertions in some genotypes). TE insertions in the DGRPs have been classified as “shared” versus “novel” depending on whether they were versus were not present in the reference genome, release 6 (Mackay et al. 2012; Rahman et al. 2015). Although nonreference insertions might not necessarily be novel, we kept the terminology from the original articles that documented TE insertions in the DGRPs, and which is used in various other studies referring to those data. First, we selected 132 of the predicted novel insertions in 12 of the DGRP lines and designed primers for the sequence flanking those insertions (supplementary table S1, Supplementary Material online). For each of the lines, we extracted gDNA from pools of ten males (homogenized using pestles), using DNeasy Blood and Tissue kit (Qiagen), following manufacturer’s instructions. We then used 4 ng of this gDNA in 15 μl long PCRs with 0.5 μM primers, 2% DMSO, 0.5 mM dNTPs mix, 0.21 μl of GoTaq enzyme (Promega). Thermocycler conditions included 2 min at 92 °C; 10 cycles of 92 °C for 10 s, 60 °C for 15 s, 68 °C for 10 min; 30 cycles of 92 °C for 15 s, 60 °C for 30 s, 68 °C for 10 min + 20 s cycle elongation for each successive cycle; 7 min at 68 °C. Amplicons were sized (1% agarose gel electrophoresis) and sequenced (ThermoFisher BigDye Terminator v1.1, or SUPREMErun™ from; same forward primers used for amplification) and these were NZYTech compared with the size and sequence of the canonical Drosophila transposons (Flybase version 9.42). Second, we tested the absence of specific TEs in genotypes annotated as having no insertions of that TE. We ran PCR with primers specific for each of seven TEs (blood, copia, gypsy5, H-element, jockey, opus, pogo; supplementary table S4, Supplementary Material online) and gDNA extracted from pools of ten adult females (extractions as described above) of eight genotypes (RAL-109, RAL-161, RAL-237, RAL-350, RAL-362, RAL-555, RAL-776, and RAL-808) predicted to not have one or more of those TEs (Mackay et al. 2012), confirming the presence or absence of insertion band in 1% agarose gel (supplementary fig. S1B, Supplementary Material online). With the gDNA from each of the target genotypes, we ran two types of positive controls: (1) with TE-specific primers with gDNA extracted from a genotype (RAL-321) predicted to have insertions of all six TEs, and (2) with primers for the Drosophila gene RPL32 (supplementary table S4, Supplementary Material online) present in every line. gDNA extracted as described above was used in 10 µl PCRs containing 0.4 ng gDNA, 0.25 U GoTaq (Promega), 1.5 mM MgCl 2 , and 0.5 μM of each primer. The thermal cycling protocol was: 10 min at 95 °C; 35 cycles of 95 °C for 30 s, 60 °C for 1 min, 72 °C for 30 s; 5 min at 72 °C. Fly Lines and Husbandry We randomly chose 25 DGRP lines described to be infected with Wolbachia (Mackay et al. 2012) and none were described to be infected with the endosymbiont Spiroplasma (Richardson et al. 2012). See the complete list of target DGRP genotypes in supplementary tables S2 and S3, Supplementary Material online. For each of the lines selected, we generated a Wolbachia-free version following procedures described in Teixeira et al. (2008) and Chrostek et al. (2013). In short, flies were first rid of Wolbachia by feeding on food supplemented with tetracycline antibiotic (0.05 mg/ml) for two generations. Their gut flora was then restored by placing sterilized eggs of Wolbachia-cleared flies (10 min in 50% bleach followed by washing in sterilized water) on food supplemented with a bacterial inoculum (150 µl of a mix prepared by mixing 2 ml of sterile water with 1 g of a month-old food filtered to remove eggs and larvae) of each respective untreated (Wolbachia-positive) line. Flies were Wolbachia-free and gut microbiota-homogenized for at least five generations before the experiments were initiated. Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 7 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023 Eugénio et al. GBE Flies were reared at 25 °C and 12 h:12 h light:dark cycle, in vials with cornmeal-agar food (45 g/l molasses, 75 g/l white sugar, 70 g/l corn flour, 20 g/l yeast extract, 10 g/l agar-agar, and 25 ml Nipagin at 10%) and similar density conditions. For our experiments, we transferred newly eclosed adult flies to vials in groups of ten females and six males. Females were sampled for extraction of DNA (for quantification of Wolbachia) or of RNA (for quantification of TE expression) at 10 days of age. Wolbachia Presence and Loads We used Wolbachia-specific primers against the Wolbachia surface protein gene (wsp; sequence from Teixeira et al. 2008) to confirm that the tetracycline-treated Wolb− lines were indeed Wolbachia free and to quantify Wolbachia loads in the untreated Wolb+ lines. We confirmed the absence of Wolbachia in each of the tetracycline-treated lines in 10 µl PCRs, containing 0.4 ng gDNA template, 0.25 U GoTaq (Promega), 1.5 mM MgCl 2 , and 0.5 μM of each primer (wsp). We used gDNA extracted (Qiagen’s DNeasy Blood and tissue kit, following manufacturer’s indications) from 3 pools of 10 females (mixed ages) from each of the 25 Wolb− lines, homogenized using Qiagen Tissue Lyser II (2 min at 23 s/f). As positive control, we extracted gDNA from the 25 Wolb+ lines (same protocol) and used those samples as template. Thermal cycle was 4 min at 95 °C; 35 cycles of 95 °C for 30 s, 60 °C for 1 min and 72 °C for 30 s; 5 min at 72 °C. PCR amplicons were checked by electrophoresis gel (1% agarose) and we confirmed the successful removal of Wolbachia (no amplicon) in all 25 target DGRP lines. We measured Wolbachia loads in 5 individual females (10 days post-eclosion) from each of our 25 target Wolb+ lines. Individual females were homogenized in Qiagen ATL Buffer in 96 well-plates with a sterile glass bead per well on a Tissue Lyser II (Qiagen) at 23 s/f for 2 min, before DNA was extracted using the Quick-DNA™ 96 kit (Zymo Research), following manufacturer’s instructions. DNA was eluted in 200 µl buffer AE from the kit and stored at −20 °C until qPCR, which was run using primers (supplementary table S4, Supplementary Material online) either for a Wolbachia gene (wsp; measuring Wolbachia load) or for a host gene (actin, proxy for number of host cells) as described below. RNA Extraction and cDNA Synthesis for TE Expression Quantification by qPCR To quantify TE expression, we extracted RNA from eight replicate pools of 10 co-housed, 10-day-old females, for each of the 25 pairs of Wolb+ and Wolb− genotypes (total of 50 lines). Whole bodies were homogenized in 400 µl TRIzol (Invitrogen) using a sterile glass bead in microcentrifuge tubes and a Tissue Lyser II (Qiagen) at 26 s/f for 1 min. Homogenates were stored at −80 °C until further processing. Once thawed, we added 80 µl of chloroform, centrifuged (12,000 × g for 15 min at 4 °C) and collected the supernatant aqueous phase containing the RNA (to avoid carrying fly tissues and fat to the RNA extraction step), and then we added 400 µl more TRIzol. Total RNA was then extracted using the Direct-zol™ 96 RNA Kit (Zymo Research), following manufacturer instructions. We used 4 µg of RNA to synthesize cDNA with NZY First-Strand cDNA Synthesis Kit (NZYTech), following manufacturer’s instructions. cDNA was then diluted 1:10 in sterile water (Sigma) to be used as template in qPCR with primers (supplementary table S4, Supplementary Material online) against each of the 14 target TEs (412, 1360, blastopia, blood, copia, Cr1a, gypsy5, Idefix, Juan, mdg1, opus, Quasimodo, P-element, pogo) or against one reference gene, EF1, chosen from a number of candidates (18S, Act5c, actin, EF1, ELF2, Gapdh1, Mnf, Rpl32, Rps20, TBP, tubulin) using Normfinder (Andersen et al. 2004) and selecting a gene with Cq values similar to that of the TEs being tested (qPCR reagents and thermocycle as described below). qPCR With Standard Curves We used qPCR to measure both Wolbachia titers (gDNA template and primers for one Wolbachia-specific gene, wsp, and one host-specific reference gene, actin) and TE expression (cDNA template and primers for each of the 14 target TEs and one reference gene). Template preparation and primers were described above. Our qPCR studies followed MIQE guidelines (Taylor et al. 2010), including technical and biological replication, ensuring template quality, careful selection of reference genes, and correction for differences in primer efficiency. Moreover, all samples being directly compared were ran together and using the same batch of reagents. For each biological replicate sample, we ran two technical replicate reactions in an QuantStudio™ 7 Flex Real-Time PCR System (Applied Biosystems™). We used 4 µl of genomic template, 0.5 µl of each primer (0.2 µM) and 5 µl of SYBR Green I® (Bio Rad), and the following thermal cycling conditions: 2 min at 50 °C; 10 min at 95 °C; 40 cycles of 95 °C for 30 s, 60 °C for 1 min and 72 °C for 30 s. We discarded biological replicates for which the standard deviation between Cq values of the two technical replicates was >0.5, and calculated the mean Cq value between technical replicates for each of all other biological replicates. Processing of Cq data of biological replicates is detailed below. For each gene and each TE, we also obtained standard curves relating amount of template and Cq values. These were obtained by using as template a 1:10 serial dilution (8 dilutions) of a cleaned (Macherey-Nagel’s NucleoSpin 8Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023 Effects of Wolbachia on Transposable Element Expression GBE Gel and PCR Clean-up) and quantified (Invitrogen’s Qubit™) PCR product (obtained by PCR on gDNA extracted from flies from the standard line Oregon R). We used the equations for the linear regression of log quantity of starting template (x-axis) and Cq value (y-axis) to: (1) do absolute quantification of wsp and actin, as there in no obvious calibrator sample for analysis of Wolbachia loads, and (2) calculate primer efficiency required for the relative quantification of TE expression using the Pfaffl method (2001). Processing qPCR Cq Data to Quantify Wolbachia and TE Expression For the absolute quantification of Wolbachia loads, we used the mean Cq values of each biological sample and the standard curves for wsp and actin to determine the quantity of each of the genes in the sample used as template: quantity = 10 ((Cq−b)/m) , where b is the intercept and m is the slope of the linear regression equation. We estimated the quantity of both wsp and actin in each sample and then calculated the ratio between the two (quantity of wsp/quantity of actin) as a measurement of Wolbachia load in relation to host cells. For the relative quantification of TE expression, we used the Pfaffl method: expression ratio = E (TE) ΔCq(TE) /E (EF1) ΔCq(EF1) . E is the amplification efficiency for each primer pair and is calculated based on the equation of the linear regression of the respective standard curve: E = 10 −1/slope (primer efficiencies in supplementary table S4, Supplementary Material online). For all TE × genotype samples, ΔCq refers to the difference in Cq values between a calibrator sample (average of same-genotype Wolb+ samples) and each sample for that genotype (Wolb− and Wolb+). Statistical Analysis All statistical analyses were performed in R (version 4.3.1), using Rstudio (version 2022.07.2). We estimated broad sense heritability (H 2 ) for Wolbachia loads as H 2 = σ 2 A/(σ 2 A + σ 2 W), where σ 2 A is the amonglines variance and σ 2 W is the within-line variance. Variance components were extracted using the VCA R package (Schuetzenmeister and Dufey 2020). We tested for differences between the 25 target DGRP genotypes in: (1) Wolbachia loads, using ANOVA with genotype as fixed factor: aov(Wolbachia load ∼ genotype) in R syntax, and (2) the number of novel TE insertions, using ANCOVA with Wolbachia load as covariate, and genotype and TE as fixed factors: aov(novel TE insertions ∼ mean Wolbachia load + genotype * TE) in R syntax. To account for variation in TE expression with Wolbachia status (Wolb+/Wolb−), we used ANOVA with TE, genotype, and Wolbachia status as fixed factors: aov(log2 TE expression normalized to reference gene expression ∼ TE * genotype * Wolbachia status) in R syntax. Then, for each TE in each paired Wolb+/Wolbgenotype, we compared TE expression between Wolb+ and Wolb− flies using ANOVA with Wolbachia status as fixed factor: aov(log2 TE expression normalized to reference gene expression of TE expression ∼ Wolbachia status) in R syntax. We plotted the residuals of the models and found that their distributions were sufficiently close to normal to justify parametric tests. However, we also applied a nonparametric test (Kruskal–Wallis test), which gave mostly equivalent results (supplementary table S3, Supplementary Material online). Supplementary Material Supplementary data are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/). Acknowledgments The authors thank Renato Alves for helping with bioinformatic tools, Daniel Sobral for aligning next generation sequencing reads against TE canonical sequences, Jorge Carneiro for crucial help with the statistical analysis, the Fly Facility at Instituto Gulbenkian de Ciência for fly maintenance, and Élio Sucena and Isabel Gordo for serving as thesis committee members on ATE’s PhD work. This work was supported by the Portuguese science funding agency, Fundação para a Ciência e Tecnologia (FCT): PhD fellowships for AT Eugénio (SFRH/BD/115535/2016, COVID/BD/ 151707/2021), PhD fellowships for MSP Marialva (SFRH/ BD/51882/2012), and research grants to P Beldade (PTDC/BIA-EVF/0017/2014, PTDC/BIA-EVL/0321/2021). Data Availability Data submitted in supplementary table S2, Supplementary Material online. Literature Cited Andersen CL, Jensen JL, Orntoft TF. 2004. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 64:5245–5250. Anderson SN, et al. 2019. Dynamic patterns of transcript abundance of transposable element families in maize. G3 (Bethesda) 9: 3673–3682. Ayarpadikannan S, Kim HS. 2014. The impact of transposable elements in genome evolution and genetic instability and their implications in various diseases. Genomics Inform. 12:98–104. Baião GC, Schneider DI, Miller WJ, Klasson L. 2019. The effect of Wolbachia on gene expression in Drosophila paulistorum and its implications for symbiont-induced host speciation. BMC Genomics. 20:465. Becking T, Gilbert C, Cordaux R. 2020. Impact of transposable elements on genome size variation between two closely related crustacean species. Annal Biochem. 600:113770. Genome Biol. Evol. 15(3) https://doi.org/10.1093/gbe/evad036 Advance Access publication 28 February 2023 9 Downloaded from https://academic.oup.com/gbe/article/15/3/evad036/7059559 by Universidade de Lisboa user on 03 July 2023