INVESTIGATION Dynamic Interactions Between the Genome and an Endogenous Retrovirus: Tirant in Drosophila simulans Wild-Type Strains Marie Fablet,*,1 Angelo Jacquet,* Rita Rebollo,† Annabelle Haudry,* Carine Rey,‡ Judit Salces-Ortiz,*,2 Prajakta Bajad,§ Nelly Burlet,* Michael F. Jantsch,§ Maria Pilar García Guerreiro,** and Cristina Vieira* *Laboratoire de Biométrie et Biologie Evolutive, Université de Lyon; Université Lyon 1; CNRS; UMR 5558, Villeurbanne 69622, France, †BF2i, Univ Lyon, INSA-Lyon, INRA, UMR 0203, Villeurbanne 69621, France, ‡Laboratoire de Biologie et Modélisation de la Cellule, UnivLyon, ENS de Lyon, Univ Claude Bernard, CNRS UMR 5239, INSERM U1210, Lyon 69007, France, §Center for Anatomy and Cell Biology, Division of Cell- and Developmental Biology, Medical University of Vienna, Vienna 1090, Austria, and **Grup de Genòmica, Bioinformàtica i Biologia Evolutiva, Departament de Genètica i Microbiologia, Universitat Autònoma de Barcelona, Bellaterra-Barcelona 08193, Spain ORCID ID: 0000-0002-7819-4541 (M.F.) ABSTRACT All genomes contain repeated sequences that are known as transposable elements (TEs). KEYWORDS Among these are endogenous retroviruses (ERVs), which are sequences similar to retroviruses and are Transposable transmitted across generations from parent to progeny. These sequences are controlled in genomes element through epigenetic mechanisms. At the center of the epigenetic control of TEs are small interfering RNAs of Retrotransposon the piRNA class, which trigger heterochromatinization of TE sequences. The tirant ERV of Drosophila Drosophila simulans displays intra-specific variability in copy numbers, insertion sites, and transcription levels, providing epigenetic us with a well-suited model to study the dynamic relationship between a TE family and the host genome control through epigenetic mechanisms. We show that tirant transcript amounts and piRNA amounts are positively piRNA correlated in ovaries in normal conditions, unlike what was previously described following divergent chromatin crosses. In addition, we describe tirant insertion polymorphism in the genomes of three D. simulans wild- type strains, which reveals a limited number of insertions that may be associated with gene transcript level changes through heterochromatin spreading and have phenotypic impacts. Taken together, our results participate in the understanding of the equilibrium between the host genome and its TEs. All genomes contain repeated sequences that are known as transpos- mutations. A few decades ago, TEs were only acknowledged as parasitic able elements (TEs). These are sequences that can move and multiply “junk DNA”; however, we now have evidence that host-TE relation- along the chromosome arms, generating most of the time deleterious ships may range from parasitism to mutualism, as do all symbiotic interactions (Biémont and Vieira 2006; Jangam et al. 2017). TEs may be of various structures (Wicker et al. 2007), among which are endog- Copyright © 2019 Fablet et al. doi: https://doi.org/10.1534/g3.118.200789 enous retroviruses (ERVs). ERVs are sequences similar to retroviruses Manuscript received October 10, 2018; accepted for publication January 14, 2019; and are transmitted across generations from parent to progeny. Due to published Early Online January 18, 2019. the potential harmful activity of TEs, mechanisms that allow to control This is an open-access article distributed under the terms of the Creative them have been favored by natural selection, and consist mainly Commons Attribution 4.0 International License (http://creativecommons.org/ of epigenetic processes (Slotkin and Martienssen 2007; Siomi et al. licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 2011). In Drosophila, knowledge of such processes is increasingly Supplemental material available at Figshare: https://doi.org/10.25387/ accumulating. At the center of the epigenetic control of TEs are small g3.7578380. interfering RNAs of the piRNA class, which trigger heterochromatini- 1 Correspondence: Laboratoire de Biométrie et Biologie Evolutive, Université de zation of TE sequences through trimethylation of histone 3 lysine Lyon; Université Lyon 1; CNRS; UMR 5558, Villeurbanne 69622, 43 bv 11 novembre 1918, France. E-mail: marie.fablet@univ-lyon1.fr 9 (H3K9me3) (Sienski et al. 2012; Le Thomas et al. 2013). 2 Current address: Institute of Evolutionary Biology (CSIC-Universitat Pompeu TE impacts on genomes have varied nature and extents, from dele- Fabra), Barcelona 08003, Spain. terious to adaptive (reviewed in Cordaux and Batzer 2009; Casacuberta Volume 9 | March 2019 | 855 and González 2013; Fablet et al. 2017). For instance, TE insertions may Using a modified version of the TEcount module of TEtools (Lerat cause diseases (Hancks and Kazazian 2016), whereas in some cases, TEs et al. 2017), we mapped the reads against a list of TE sequences may be recruited by the host genome and perform endogenous functions. either in the sense or antisense directions. The list of TE sequences The most famous example is the syncytin gene in mammals, which is an is made of the TE insertions retrieved from D. simulans sequenced ancient TE gene and is now essential for placenta formation (Mi et al. genome (Lerat et al. 2017). Read count numbers were normalized to 2000). TEs may also be involved in adaptation. As an example, the Accord miRNA read count numbers (miRNA sequences retrieved from TE inserted into the 59 end of the Cyp6g1 gene of Drosophila mela- FlyBase: dsim-all-miRNA-r2.02.fasta.gz). nogaster is found in some strains and provides them with increased We looked for ping-pong signatures using signature.py with the resistance to insecticides (Daborn et al. 2002). options min_size = 23 and max_size = 30 (Antoniewski 2014). We Despite their potential harmful impacts, TEs may reach high pro- used as input files SAM alignments obtained from cleaned small portions of genome sequences (Biémont and Vieira 2006; Tenaillon RNA reads against our C-subfamily complete (8.5 kb) reference et al. 2010; Elliott and Gregory 2015). Therefore, it is of fundamental tirant sequence: accession number AC0054444, positions 50,203 to interest to precisely understand TE dynamics within genomes, espe- 58,729, extracted from the D. melanogaster genome, as already de- cially at the transcriptional level, which is the preliminary stage before scribed in Fablet et al. (2006). We used bowtie (–best) (Langmead mobilization. In this study, we propose to use the tirant element of et al. 2009). We analyzed base composition at each position for Drosophila simulans to tackle this important issue. piRNAs aligned against tirant using SAMStat (Lassmann et al. 2011). Tirant belongs to the ERV class of retrotransposons (Terzian et al. To visualize read mapping along tirant sequence, we sampled 2001). It is made of three open reading frames (ORFs), namely gag, pol, 40,000,000 reads from each small RNA-seq sample using the fastq- and env, which provide all proteins necessary to fulfill the retroviral sample program from fastq-tools 0.8 (https://homes.cs.washington. cycle, and is bordered by two long terminal repeats (LTRs), which edu/dcjones/fastq-tools/), and then we selected 23-30 nt-long reads include regulatory signals. Two subfamilies can be distinguished based as mentioned above. We used bowtie (–best) (Langmead et al. 2009) to on their sequence similarity and the structure of a minisatellite located map these small RNA-seq reads against our C-subfamily reference in the 59 untranslated region (UTR) (Fablet et al. 2006). Subfamily C is tirant sequence (see above). We removed the 39 LTR, so that reads responsible for tirant activity while subfamily S is heterochromatic, corresponding to LTRs should map at a unique location. We filtered low-copy-number, and silent (Fablet et al. 2006, 2009). Tirant was first out alignments with a mapq score below 20 using samtools (Li et al. discovered in D. melanogaster (Garrell and Modolell 1990) and found 2009). Visualization of the alignments was performed using BamView to occur in related species (Fablet et al. 2007). Its number of copies and (Carver et al. 2010). insertion sites, together with expression profiles vary across D. simulans wild-type strains (Vieira et al. 1999; Fablet et al. 2006; Akkouche et al. Cell culture and reporter gene assays 2012), providing a powerful system to investigate the impacts of We used reporter gene assays to test the sense and antisense promoter ERVs on the host genome, particularly regarding neighboring gene potentials of the LTRs of both subfamilies. LTR sequences were PCR expression and epigenetic signatures. amplified from the Mayotte strain (see Supplementary Material S1 The analyses of RNA-seq and small RNA-seq data indicate that for primer sequences). Primers were designed to include KpnI and tirant transcript amounts and piRNA amounts are positively correlated HindIII restriction sites, to allow subsequent directional cloning into in ovaries in normal conditions, while they were previously found to the pGL4.10 plasmid (Promega, luciferase reporter gene) using T4 be negatively correlated in some divergent crosses leading to tirant DNA ligase (NEB). S2 Drosophila cells were transfected using Cell- deregulation. These results allow us to propose that the shape of the fectin (Invitrogen) and lyzed for analysis after 24 h of incubation correlation between the amounts of TE transcripts and the correspond- (6 replicates per construct). Luminescence was then estimated using ing piRNAs are indicative of genome stability. In addition, we describe the Dual-Luciferase Reporter Assay System (Promega). tirant insertion polymorphism in the genomes of three D. simulans wild-type strains, which reveals a limited number of insertions that may Genome sequencing and assemblies be associated with gene transcript level changes through heterochro- Genomic DNA samples were treated with RNAse A at 37° for 30 min matin spreading and have phenotypic impacts. Taken together, our and purified using Phenol:Chloroform extraction. The purified geno- results participate in the understanding of the equilibrium between mic DNA samples were fragmented using Covaris S220 sonicator to the host genome and its TEs. get fragments of size 300-800 bp. Sequencing libraries were prepared from these fragmented DNA samples using NEBNext UltraTM DNA MATERIAL AND METHODS Library Prep Kit (Illumina, # E7370L). All three libraries were sequenced in one lane of an Illumina HiSeqV4 with 125 bp paired-end reads at Drosophila stocks the Next Generation Sequencing unit of the Vienna Biocenter Core Drosophila strains were maintained in the laboratory at 24° as small- Facilities (VBCF http://vbcf.ac.at). mass cultures. Chicharo, Mayotte, Makindu and Zimbabwe are wild- Approximately 100 million read pairs (24 Gb) from three accessions type strains, anciently sampled from the field; w501 is the major strain (Mayotte, Makindu and Chicharo) were obtained and raw reads were of the D. simulans sequenced genome. filtered using the unsupervised quality trimming program UrQt (Modolo and Lerat 2015). Contigs were then generated by the Ray Small RNA analysis assembler v2.3.1 (Boisvert et al. 2010) using kmer length varying from We used small RNA-seq data obtained from ovaries from Akkouche 25 to 45. Best assemblies were selected by maximizing both NG50 and et al. (2013) and Lerat et al. (2017), corresponding to the Chicharo, the number of contigs larger than 10 kb. For each accession, the Mayotte and Makindu strains. We first removed adapter sequences assembled genome size was over 100% of the expected genome size using cutadatp (Martin 2011) -a CTGTAGGCACCATCAA. Using (146.7 Mb). The final contig NG50 were 5.38 kb (k = 43), 11.89 kb (k = PRINSEQ lite version 0.20.4 (Schmieder and Edwards 2011), we 39) and 7.36 kb (k = 35), with maximum scaffold length of 130 kb, filtered reads of size 23 to 30 nt and considered these as piRNAs. 495 kb and 495 kb for Mayotte, Makindu and Chicharo, respectively. 856 | M. Fablet et al. In order to assess completeness of genome assemblies, we searched of the strains, apytram will be able to catch it. Then, to analyze Hs6st for conserved genes across arthropods using BUSCO (Simão et al. and tkv transcript structures in our samples, we aligned output tran- 2015). Among the set of 2,675 genes, 2,254 (84%), 2,403 (90%) and scripts on D. simulans genome (ASM75419v3.41) using exonerate (Slater 2,269 (85%) were recovered as complete and single copy in the and Birney 2005). Finally, we manually checked the hit distribution assemblies of Mayotte, Makindu and Chicharo, respectively. against D. simulans genome to search for the presence of tirant in tran- script sequences. We used TopHat2 (Kim et al. 2013) and ggsashimi Genome and transcriptome analyses (Garrido-Martín et al. 2018) to view transcript structures using sashimi We ran RepeatMasker (Smit et al. 2013) on the above assemblies plots. (default parameters, -species Drosophila) and retrieved “Tirant” accessions in the output file in order to identify tirant insertion sites. ChIP assays Then, we extracted the flanking sequences for the insertions found Biological duplicates were produced. For each replicate, seventy above and located them on droSim1 assembly using the BLAT tool ovary pairs were dissected from three to five day-old females and (Kent 2002) of UCSC Genome Browser (http://genome.ucsc.edu). flash frozen to maintain chromatin integrity. Chromatin extraction We checked that shared tirant insertions had the same breakpoints and sonication were carried out as described in Akkouche et al. in all strains. (2013) with the following modifications: 30 sonication cycles RNA-seq data from ovaries and the corresponding computed read (30 sec on and 30 sec off on high) were performed to obtain chro- counts for all TE families, including tirant, were obtained from Lerat matin fragments comprised between 200 and 500 bp. The chroma- et al. (2017) for five strains (Chicharo, Makindu, Mayotte, Zimbabwe, tin fragments were then aliquoted and diluted into six IP samples of and w501). Differential expression analysis was performed using 500 mL along with a 10% Input sample. IP samples were then immu- DESeq2 (Love et al. 2014), which model internally corrects for library noprecipitated following the manufactor’s protocol (Magna ChIP, size, and provides normalization. Read counts obtained for genes on Millipore) with 3 mg of antibodies (H3K4me3: Abcam # ab8580; one hand and TEs on the other hand (Lerat et al. 2017) were concat- H3K9me3: Actif Motif # 39161). enated to make the complete read count table that was analyzed at Quantitative PCR (qPCR) was used to estimate IP enrichment as once using DESeq2. previously described (Fablet et al. 2009; Akkouche et al. 2013) and The tkv gene had no assigned D. simulans ortholog in the FlyBase performed as technical triplicates. For each of the studied tirant gene list that we used for RNA-seq analysis at that time (Lerat et al. 2017). insertions, we used qPCR primers located at the insertion site and In 2018, FBgn0194057 has been assigned as D. simulans tkv ortholog 2 kb away, in both directions. At the insertion site, we designed (OrthoDB v9.1). Thus, we retrieved the FBgn0194057 sequence, per- primer pairs displaying one primer within tirant and the other formed the same procedure as for the other genes (Lerat et al. 2017), primer in the flanking sequence. For strains devoid of tirant inser- and added the corresponding read counts to our count table. tion, we used the same primer in the flanking region as mentioned To visualize read mapping along tirant sequence, we sampled above and designed a second primer in the flanking region beyond 90,000,000 reads from each DNA-seq sample and 25,000,000 reads the site of tirant insertion. These two primer pairs allow us to am- from each RNA-seq sample (one biological replicate from each strain) plify tirant insertion sites in strains where tirant is present as well as to ensure even abundance of reads. In addition, for RNA-seq data, we in strains were tirant is absent. Primer sequences are provided in only used reads longer than 20 nt (PRINSEQ filter, as described above). Supplementary Material S1. As controls, we used light for H3K9me3 We used the fastq-sample program from fastq-tools 0.8 (https://homes. enrichment (Fw: 59-GCT AGG CAA TGA CAA AGT CCT TTG cs.washington.edu/dcjones/fastq-tools/) and the fastq_sampler.py GG-39 and Rv: 59-GCA TTC GTC TGA AGT CGG CAG ATA script (Modolo 2018). Then, we used bowtie2 (–very-sensitive) G-39), and rpl32 for H3K4me3 enrichment (Fw: 59-CAG CTT (Langmead and Salzberg 2012) and mapped DNA-seq and RNA-seq CAA GAT GAC CAT C-39 and Rv: 59-GTT CGA TCC GTA ACC reads against our C-subfamily reference tirant sequence (see above). GAT GT-39). Raw data from qPCRs are provided in Supplementary We removed the 39 LTR, so that reads corresponding to LTRs should Material S2. map at a unique location. We filtered out alignments with a mapq score below 20 using samtools (Li et al. 2009). Visualization of the Long PCR alignments was performed using BamView (Carver et al. 2010). We PCR amplified tirant insertion at the Hs6st locus using DreamTaq Transcript structure analysis: we used apytram (apytram v1.1 (Fermentas), with a 7 min elongation time and a 68° hybridization. (Rey et al. 2017)) to assemble independently Hs6st and tkv in Chi- Primer sequences are 59- GCAACACTGACAGCAACTACA -39 and charo, Makindu and Mayotte RNA-seq samples (replicates #1). 59- TCCTTGCTAGCTACATGGAAC -39. Apytram is an implementation of the Target Restricted Assembly Method (Johnson et al. 2013). This method allows to focus on the assembly of only one gene of a transcriptome and not the whole Data Availability transcriptome. It is based on an iterative process: at the first itera- Genomic raw reads are available under the SRA accession number tion, a reference gene from another species is used as bait sequences SRP128969. Genome assemblies, modified versions of TEcount.py and to fish reads in RNA-seq data using BLAST (Camacho et al. 2009). the list of TE sequences are available at ftp://pbil.univ-lyon1.fr/pub/ Paired reads are used to enlarge this batch of reads. Reads are then datasets/Fablet2018. SupplMat S1 contains primer sequences. Sup- de novo assembled using Trinity (Grabherr et al. 2011), and a new plMat S2 contains ChIP qPCR results. SupplMat S3 is an alignment iteration begins with the reconstructed sequences as baits. Here, we of reference C and S sequences. Alignments of piRNAs against tirant used the Hs6st and tkv genes from D. melanogaster as references. reference sequence are provided Figure S1. Figure S2 displays graphics This iterative process allows to reconstruct transcripts step by of piRNA base composition. Gels obtained from PCR products are step from homologous regions between D. melanogaster and shown Figure S3. Sashimi plots for tkv and Hs6st splicing are provided D. simulans, and thus non homologous regions can be reconstructed Figure S4. Supplemental material available at Figshare: https://doi.org/ in D. simulans transcripts. In addition, if ever tirant is exonized in one 10.25387/g3.7578380. Volume 9 March 2019 | Tirant ERV Dynamics in Drosophila | 857 RESULTS being the strongest promoter and LTR S sense being the weakest. In addition, we find that both subfamily LTRs display weak antisense Tirant is differentially expressed across wild- promoter activity (Figure 3) (Wilcoxon tests, C antisense vs. S sense: type strains p-value = 0.002, S antisense vs. S sense: p-value = 0.002, C antisense Using RNA-seq data obtained from ovary samples, we were able to vs. S antisense: p-value = 0.485). Therefore, the results suggest that, quantify tirant transcript levels in D. simulans strains (Figure 1). As although subfamily S is not significantly involved in the production expected (Fablet et al. 2006, 2009), subfamily S transcript levels of tirant transcripts, and albeit the weak activity observed in vitro, were found at background level in all strains, and so were sub- it may rather be involved in the production of antisense piRNA family C transcript levels in w501, Chicharo and Zimbabwe strains. precursors. Tirant produces significant amounts of transcripts only in Makindu and Mayotte strains (Figure 1A). The mapping of tirant reads against Tirant displays insertional polymorphism the reference sequence shows enrichment in the env region (Figure 1B). The influence of a TE on the host genome depends on its insertion sites. This pattern could reflect either the existence of a master copy deleted To investigate tirant’s impacts, we produced genome assemblies for Chi- for gag and pol, or the prevalence of splicing events of the gag and charo, Makindu and Mayotte, and used RepeatMasker (Smit et al. 2013) pol ORFs, as described by Marsano et al. (2000) in D. melanogaster. to retrieve tirant insertion sites. We excluded scaffolds exclusively made The higher abundance of env reads in the Mayotte strain is in agree- of tirant. Nested tirant insertions into other TE sequences were also ment with the production of the tirant Env protein in this strain, excluded since we are not able to locate them precisely on the genome. which we previously observed using immunostaining (Akkouche The final number of retrieved tirant insertions is congruent with our et al. 2012). knowledge of tirant copy numbers (Fablet et al. 2006, 2009; Akkouche et al. 2012), indicating that we did not missed a significant number of loci. Tirant displays specific sense and antisense piRNAs Overall, we localized precisely 15 insertions, the majority of which are We specifically analyzed small RNAs for the two strains that show unique to one strain (Table 1). Among these, six are located in introns significant tirant transcript levels (i.e., Makindu and Mayotte), and and two upstream transcription start sites (TSSs), which are the inser- one strain presenting background levels (i.e., Chicharo). We found tions the most likely to have an impact on genes. Insertion #10 is sense and antisense tirant-specific piRNAs, corresponding to both located 242 bp upstream of CanA1 TSS, and insertion #11 533 bp the C and S subfamilies. Sequence similarity between subfamilies C upstream of Lamp1 TSS. We identified five additional insertions, at- and S varies along tirant sequence: it is the lowest in the UTR 59 tributed to the U chromosome (Unknown chromosome), and which region (65% (Fablet et al. 2006)) and displays higher values within likely correspond to heterochromatic regions. Among these, one is coding sequences (76–85% global similarity rate (Fablet et al. shared across strains. 2009)). Although we cannot exclude that some reads may be mis- It is to be noted that, following this procedure, the tirant insertion into attributed to the other subfamily, it should not happen for the the tkv gene (insertion #12) is only found in the Makindu genome, while majority. Except for a peak in the LTR region, tirant piRNAs evenly we have experimental evidence that it also exists in the Mayotte strain. map along tirant reference sequence in all three strains (Figure S1). This stresses on the non exhaustive nature of this approach, which Tirant piRNAs displayed comparable proportions of antisense strongly relies on scaffold assemblies. Indeed, we used a PCR approach piRNAs to total piRNAs in all three strains and tirant subfamilies, to confirm tirant insertion polymorphism for insertions #3 and #12, ranging from 0.798 to 0.869 (Figure 2A). In addition, we looked for which are studied in more details below. (See Supplementary Material S1 ping-pong signatures, which are 10 nt overlaps between sense and for primer sequences). For insertion #3, PCR profiles were congruent antisense piRNAs (Brennecke et al. 2007). We found clear, similar with tirant being present only in Mayotte and absent in Makindu and enrichments for 10 nt overlaps in all three strains (Figure 2C), in- Chicharo (Figure S3A). For insertion #12, PCR profiles indicated dicative of similar, functional, secondary piRNA pathways against tirant presence in Makindu and Mayotte, and tirant absence in tirant in all three strains. In addition, these piRNAs showed Chicharo (Figure S3B). However, we could also detect a very faint expected enrichment in 1U and 10A (Figure S2) (Saito et al. 2006; band indicating tirant absence in Makindu. This suggests that this Brennecke et al. 2007). insertion may be heterozygous or even absent in few individuals. Moreover, our data allowed us to reveal a positive correlation be- Accessing the precise structure of each tirant insertion is difficult tween tirant piRNA read counts and tirant transcript read counts without long reads. Therefore, we mapped genomic reads against our (Pearson correlation, r = 0.963, p-value = 0.002, Figure 2B). The cor- tirant reference sequence to have an idea of the representativeness of relation remains strong using only sense (Pearson correlation, r = 0.965, each region within the genomes of our strains. We found that reads map p-value = 0.002), or antisense (Pearson correlation, r = 0.946, p-value = evenly along the reference sequence for Chicharo and Makindu, sug- 0.004) piRNAs. gesting the presence of full-length copies (Figure 4). Mayotte displays the highest number of reads after normalization, in agreement with a higher Tirant LTRs display antisense promoter activity number of tirant copies. We note a slight depression in the center of the As new TE insertions may behave as dual-strand piRNA clusters pol and env ORFs (Figure 4), suggesting the existence of internally de- (Mohn et al. 2014; Shpiz et al. 2014), we wondered whether tirant leted copies. The drop of mapping coverage within the 59 UTR region is copies could be able to produce antisense piRNA precursors. To due to our mapq filter (see Material and Methods), which prevents investigate such hypothesis, we tested the C and S LTRs for their multi-mapping against repeated regions, such as the minisatellite that sense and antisense promoter abilities. We known from a previous is described within tirant 59 UTR (Fablet et al. 2006). work that sequence similarity between C and S LTRs is 65% on average, while sequence similarity within each subfamily is higher Mild impact of tirant on neighboring genes than 98% (Fablet et al. 2006). It appears that the different constructs For those insertions located near or within genes (within introns), we behave significantly differently regarding gene expression promo- tested whether tirant could have an effect on gene expression through tion (Kruskal Wallis test, df = 3, p-value = 2.1024), with LTR C sense modifications of chromatin structure. To do so, we used data produced 858 | M. Fablet et al. Figure 1 Tirant trancripts analyzed from RNA-seq data. A. Tirant transcript normalized read counts. B. Mapping of tirant reads against the reference tirant copy (subfamily C). The upper part is tirant structure, to scale (in bp). The lower part is read coverage along tirant sequence, obtained from samples of 25,000,000 paired-end reads. We removed LTR 39 to get rid of multi-mapping issues. from ovaries, and we focused on the insertions of tirant that fulfill the Since tkv mutations are known to alter wing vein phenotypes following criteria: i) the gene into which tirant is inserted displays a (Lindsley and Zimm 1992), we examined wing morphologies in our high enough number of reads in our RNA-seq dataset (. 100 reads), ii) flies (we did not separate sexes). We found seven flies out of 97 in the its expression pattern is contrasted across strains with and without the Mayotte strain which displayed thick L5 veins (Figure 6C, 6D), tirant insertion (Lerat et al. 2017). We found only two of such inser- while we did not find any out of 100 Makindu flies nor 100 Chicharo tions (Figure 5). The two insertions of interest are located within the flies (Fisher exact test, p-value = 0.014). It is to be noted that Hs6st and tkv genes, respectively (insertions #3 –only in Mayotte– and Mayotte is the strain that shows the strongest reduction in tkv #12 –in Mayotte and Makindu– in Table 1, respectively). In both cases, transcript levels. gene transcript levels are lower for strains where tirant is present com- As it was recently shown that piRNAs may also regulate TEs pared to strains devoid of tirant at the considered locus (DESeq2 padj through splicing (Teixeira et al. 2017), we analyzed transcript struc- values: Hs6st: Chicharo vs. Mayotte: 0.0157, Makindu vs. Mayotte: 0.0916; tures in Hs6st and tkv. Hs6st displays virtually no transcript isoform tkv: Chicharo vs. Makindu: 0.0617, Chicharo vs. Mayotte: 2.1028). More- variation in the reference D. melanogaster genome. Unsurprisingly, over, FlyBase D. simulans genome annotation indicates the presence we found the described Hs6st transcript to be produced in all three of additional genes within Hs6st (FBgn0190828, FBgn0270486, strains (Figure 6E). In contrast, four tkv isoforms are described in FBgn0190826, and FBgn0269631) and tkv (FBgn0270637, FBgn0194717, D. melanogaster. Isoforms B and C are not found to be produced FBgn0270509, and FBgn0194718) introns. However, virtually none of the in our samples, and isoform A displays the same reference structure reads of our data sets map against these sequences. in all three strains. Variability was observed regarding isoform D: We investigated chromatin structure to look for epigenetic modifica- the first exon is 1,371 bp longer in the Makindu strain compared to tions that could spread on neighboring regions. We performed ChIP the others (Figure 6F). Sashimi plots are provided Figure S4. In- experiments on H3K4me3 and H3K9me3. The former is characteristic of terestingly, this variable splice site is located less than 800 bp from open chromatin conformation, particularly found at promoters; the later is tirant insertion in Makindu and Mayotte, and thus in the region the heterochromatic mark known to be involved in TE transcriptional where we find chromatin structure to be modified in association silencing. We used qPCR primers located at the insertion site and 2 kb with tirant insertion. It is thus tempting to speculate that tirant- away, in both directions (Figure 6A, 6B). Primers at the insertion sites induced chromatin modifications may affect splice site definition. In were the same as those mentioned above to detect tirant presence (Sup- addition, we note that we do not find evidence of tirant exonization plementary Material S1). No obvious effect of tirant presence on either neither in Hs6st nor in tkv. chromatin mark at the Hs6st gene is observed: Mayotte enrichment levels (tirant present) lie in between Chicharo and Makindu levels (tirant DISCUSSION absent). On the contrary to Hs6st insertion, it appears that tirant neigh- borhood is moderately enriched in H3K9me3 at the tkv insertion site The repertoire of tirant-specific piRNAs compared to strains without tirant (fold change = 1.43 ; t-test, p-value = Our in-depth analysis of tirant-specific small RNAs revealed the presence 0.098). In addition, these histone post-translation modifications spread at of both sense and antisense piRNAs. While the primary piRNA pathway least 2 kb away downstream tirant (fold change = 1.41 ;t-test, p-value = exclusively produces antisense piRNAs, the secondary piRNA pathway 0.002). Regarding H3K4me3 enrichments, we could not detect differ- leads to the production of both sense and antisense piRNAs, character- ences between strains with and without tirant at the tkv locus. ized by the so-called ping-pong signature (Brennecke et al. 2007). Volume 9 March 2019 | Tirant ERV Dynamics in Drosophila | 859 Figure 2 Tirant-specific piRNAs. A. Relative proportions of sense and antisense piRNAs. B. Posi- tive correlation between tirant tran- script normalized read counts and tirant piRNA read counts in wild- type strains. C. Ping-pong signa- tures: there is a clear enrichment in 10 nt overlaps. Thus, the detection of sense tirant piRNAs indicates that the secondary However, we cannot exclude the involvement of the primary pathway as pathway is involved in tirant control, as already proposed by previous well, especially considering the high proportions of antisense tirant- experimental analyses (Akkouche et al. 2013). This is reinforced by the specific piRNAs. identification of ping-pong signatures, which are similar across strains. So far, we do not know which piRNA clusters are involved in the production of tirant-specific piRNAs in D. simulans strains; however, our observation of antisense promoter properties for both subfamilies, albeit weak, indicates that some copies may have the ability to behave as dual-strand piRNA clusters, as proposed by two research studies (Mohn et al. 2014; Shpiz et al. 2014). While the S subfamily is not transcribed and apparently not involved in tirant activity, this potential role in piRNA precursor production may also explain why it is found to be conserved across strains (Fablet et al. 2006). However, although we provide in vitro evidence that some tirant LTR sequences may behave as antisense promoters, we do not know whether the antisense tirant piRNAs that we detect are indeed produced from these promoters or from promoters outside of the element itself. Moreover, while we observed significant differences in tirant piRNA amounts across strains, we could not detect differences in sense vs. antisense piRNA ratios neither in ping-pong signatures, indicating that tirant control by the piRNA pathway may take place the same way and with comparable efficiencies in the three considered strains. Figure 3 Tirant LTR promoter activity. LTRs of C and S subfamilies in Dynamics of tirant control by piRNAs sense and antisense orientations were cloned upstream of a luciferase Our survey in wild-type strains revealed a positive correlation between reporter gene and tested for promoter activity in S2 cells. the amounts of tirant transcripts and tirant piRNAs in ovaries in 860 | M. Fablet et al. n Table 1 Tirant insertion sites as identified using RepeatMasker on genome assemblies. +: insertion found in the genome assembly of the strain. Due to the fragmented nature of our assemblies, we are not able to interpret “absence of detection” as “absence of tirant”. These data do not distinguish tirant subfamilies. TSS: Transcription Start Site. U: Unkown chromosome. ‡ LTR: we detected one LTR (417 bp) at the extremity of a scaffold, which means that we do not know whether the insertion is a solo LTR or longer. Insertion #10 is located 242 bp upstream of CanA1 TSS, and insertion #11 533 bp upstream of Lamp1 TSS. D. simulans orthologs as assigned in OrthoDB v9.1 (retrieved via FlyBase in november 2018) are mentioned in brackets and italics. Insertion # Chr Annotation Insertion size (bp) Mayotte Makindu Chicharo 1 3R Intron CG34383 (FBgn0190438) 41 + + + 2 3L 102 + 3 3R Intron Hs6st (FBgn0041661) $ LTR + 4 2L Intron vri (FBgn0194720) $ LTR + 5 3R Intron CG11873 (NA) $ LTR + 6 3R . 95 + 7 2L 446 + + 8 3h 490 + 9 2h . 589 + 10 3R Near CanA1 TSS (FBgn0192983) . 96 + 11 2L Near Lamp1 TSS (FBgn0195671) 704 + + + 12 2L Intron tkv (FBgn0194057) $ LTR + 13 2L . 100 + 14 2R . 95 + 15 2R Intron CG17684 (FBgn0268563) 8014 + 16 U . 99 + + + 17 U .98 + + 18 U 115 + 19 U . 93 + 20 U 438 + normal conditions. This is in agreement with what we found using all We propose that the nature of the correlation between TE transcripts TE families in the same strains (Lerat et al. 2017). This is also congruent and piRNAs reflects genome stability. In normal conditions – which with data obtained by others from D. melanogaster laboratory strains in corresponds to a balanced relationship between the genome and its normal conditions (Kelleher and Barbash 2013), and in accordance TEs, the correlation is positive. When the equilibrium is broken – with the accepted model claiming that TE transcripts fuel the ping- and so far we do not know which sensor could detect such event –, pong loop (Senti et al. 2015). the correlation becomes negative. Such patterns were also observed in Nevertheless, these results contrast with what was found in other the case of the I non-LTR retrotransposon in D. melanogaster strains studies. Notably, in a previous work, in order to understand tirant (Chambeyron et al. 2008). Indeed, Chambeyron et al. described that control, we performed crosses between Makindu and Chicharo, strain JA has no I transcripts nor I-specific piRNAs while strain HT2 because they display contrasted tirant copy numbers and activities displays high amounts of both I transcripts and piRNAs (Chambeyron (Akkouche et al. 2013). When the mother comes from Makindu et al. 2008). This corresponds to the positive correlation we propose at and the father comes from Chicharo, tirant is properly controlled equilibrium. On the other hand, when dysgenic crosses are performed in the progeny. This direction of cross is called “RT” (Regulated and disrupt genome / I stable relationship, an opposite correlation is Tirant). However, when the reciprocal cross is performed, tirant observed: SF females –which suffer from sterility– accumulate I tran- shows a strong accumulation of transcripts in the somatic follicle scripts due to a lack of I piRNAs; on the contrary, RSF females –in cells of the progeny’s ovaries (Akkouche et al. 2013). This direc- which fertility is restored– lower I transcript levels due to a high pro- tion of cross is called “NRT” (Non-Regulated Tirant). This disrup- duction of I piRNAs. Thus, the results obtained from tirant in the tion of tirant control is associated with a lack of tirant-specific present work may reflect a behavior common to all TE families. piRNAs, while these piRNAs remain relatively abundant in the RT direction of cross (Akkouche et al. 2013). It results in a negative Tirant’s soft influence on neighboring genes correlation between tirant transcript amounts and tirant-specific In most organisms, TE control is mainly achieved by epigenetic piRNA amounts. mechanisms (Slotkin and Martienssen 2007; Siomi et al. 2011). In Figure 4 Mapping of genomic reads against tirant reference sequence. The upper part represents tirant structure, to scale (in bp). The lower part depicts read coverage along tirant sequence (subfamily C), obtained from samples of 90,000,000 paired-end reads. We removed LTR 39 to get rid of multi-mapping issues. Volume 9 March 2019 | Tirant ERV Dynamics in Drosophila | 861 Figure 5 RNA-seq results for genes at or near tirant insertions. Normalized read counts. Please refer to Table 1 for tirant insertion details. Error bars are standard deviations on biological duplicates. When tirant insertion is polymorphic, we provide adjusted p-values comparing strains with and without tirant insertion, as calculated by DESeq2. Gene CG17684 (insertion #15) has no D. simulans ortholog in the list of genes used for RNA-seq data mapping. D. melanogaster, these essentially consist in piRNAs and H3K9me3 tirant insertion, and could not find any. This suggests that the (Saito et al. 2006; Aravin et al. 2007; Sienski et al. 2012; Le Thomas effects observed on chromatin structure and gene expression et al. 2013). We know that this heterochromatic mark may spread to may be attributed to the insertion of tirant. Moreover, we observed flanking regions up to a 5 kb distance (Rebollo et al. 2011). Rebollo et al. a thick veins phenotype in 7% of Mayotte flies while we could not used mouse strains and did not find that heterochromatin spreading observe it in Makindu nor Chicharo flies. Mayotte is the strain recurrently affected gene activity, probably due to the fact that such displaying the strongest reduction in tkv transcript levels (Mayotte TE insertions would be too deleterious (Rebollo et al. 2011). In a vs. Chicharo: log2FC = -0.80). Although we cannot exclude other Drosophila system, Sienski and colleagues also observed such spread- differences in genetic backgrounds, we propose that the tirant ing, which in addition did impact neighboring genes in mutants in cell insertion into tkv in Mayotte may be responsible for this thick cultures (Sienski et al. 2012). In the present study, we used the tirant veins phenotype, through tkv expression reduction induced by system to investigate the spreading of TE chromatin marks in a natural heterochromatin formation. Nevertheless, variation in tkv tran- context. We found a large level of insertion polymorphism across the script levels is observed –probably associated with micro-environmental strains we investigated (Table 1). We identified two cases, correspond- variability–, which may explain why most flies do not show the ing to insertions into introns of the Hs6st and tkv genes, associated with phenotype. changes in gene transcript levels. While we could not detect any effect In addition, as it was recently shown that piRNAs may also of tirant insertion on chromatin structure in the case of Hs6st, tirant regulate TEs through splicing (Teixeira et al. 2017), we analyzed insertion into the tkv gene was associated with a moderate enrichment transcript structures in Hs6st and tkv. We found no effect of tirant in H3K9me3 marks at the insertion site and up to 2 kb away from the insertion in Hs6st transcript structures. This observation parallels insertion site, compared to the strain devoid of tirant. In addition, we the absence of tirant effect on Hs6st chromatin conformation. In noted that the tirant effect that we detect in Makindu may be under- contrast, we observed variation in tkv splice sites nearby tirant estimated since there may be heterozygosity or absence of the insertion insertion in the Makindu strain: the first exon of isoform D is in some individuals. We searched for other TEs in the genome 1,371 bp longer than expected. It is tempting to speculate that assemblies of all strains in a 10 kb window around this particular tirant-induced chromatin modifications disrupt splice site definition; 862 | M. Fablet et al. Figure 6 Analysis of two tirant insertions within introns. A. and B. Histone mark enrichments at two tirant insertion sites. P-values in the upper left corners are produced from t-tests between strains with the tirant insertion (yellow) and strains without tirant (uncolored). Left panel: A. tirant insertion #3 into Hs6st, right panel B.: tirant insertion #12 into tkv. Genes are drawn to scale in the middle of each panel; thin lines: introns, thick boxes: exons. H3K4me3 and H3K9me3 enrichments are shown in green and red, respectively, and quantified relative to control genes (see Material and Methods section). They are displayed in boxes corresponding to the location of the qPCR amplicons: either at the tirant insertion site or 2 kb away. Biological duplicates were produced for ChIP experiments. “+” indicates the mean between biological replicates. Normalized transcript levels for the consid- ered genes are provided in the upper right corners (as already shown Figure 5). Error bars are standard deviations. Strains in which tirant is present at the considered locus are highlighted in yellow. Statistical significance for RNA-seq results was assessed by padj values provided by the DESeq2 analysis. Padj values: Hs6st: Chicharo vs. Mayotte: 0.0157, Makindu vs. Mayotte: 0.0916, tkv: Chicharo vs. Makindu: 0.0617, Chicharo vs. Mayotte: 2.1028. C. and D. Vein phenotype in the Mayotte strain. C. We observed thick L5 veins (red arrows) in seven out of 97 flies in the Mayotte strain. D. Wild-type vein phenotype observed in 90 out of 97 flies in Mayotte (the same phenotype is also observed in all examined flies in Makindu and Chicharo). E. and F. Transcript structures. Reference structures from D. melanogaster, as retrieved from the UCSC Genome Browser, are depicted in dark green, thin lines: introns, thick boxes: exons. Reconstructed exons are shown in light blue. E. Hs6st transcript structures. All three strains display the same transcript structure as described in D. melanogaster. F. tkv transcript structures. We observe a Makindu-specific transcript of isoform D. however, we do not find the same pattern in the Mayotte strain, Tirant insertions have either no or moderate effects on nearby gene which also displays tirant insertion. Further investigation is expression. The insertion into tkv may have an impact on chroma- needed at a larger scale to determine the impacts of TEs on splice tin structure and gene expression whereas the tirant insertion into site definition. Hs6st does not. Our results are congruent with the recently published Volume 9 March 2019 | Tirant ERV Dynamics in Drosophila | 863 study of Lee and Karpen, although regarding a different histone LITERATURE CITED mark. Based on a D. melanogaster genome-wide study, they report Akkouche, A., T. Grentzinger, M. Fablet, C. Armenise, N. Burlet et al., the spreading of TE-associated H3K9me2 marks on regions flank- 2013 Maternally Deposited Germline PiRNAs Silence the Tirant Ret- ing the insertions and an effect on the expression level of flanking rotransposon in Somatic Cells. EMBO Rep. 14: 458–464. https://doi.org/ genes only for about half of the TE insertions (Lee and Karpen 10.1038/embor.2013.38 Akkouche, A., R. Rebollo, N. Burlet, C. Esnault, S. Martinez et al., 2017). So far, we have not identified which properties of a TE 2012 Tirant, a Newly Discovered Active Endogenous Retrovirus in insertion make it affect nearby genes or not. We know from a Drosophila Simulans. J. Virol. 86: 3675–3681. https://doi.org/10.1128/ previous work that tirant is full length at the tkv locus in the JVI.07146-11 Makindu strain (Fablet et al. 2009). We performed a long PCR Antoniewski, C., 2014 Computing SiRNA and PiRNA Overlap Signa- at the Hs6st locus and observed a band congruent with a full- tures. Methods Mol. Biol. 1173: 135–146. https://doi.org/10.1007/978- length insertion as well (Figure S3C). In addition, both tirant 1-4939-0931-5_12 sequences are oriented the same way as the gene. We may note Aravin, A. A., G. J. Hannon, and J. Brennecke, 2007 The Piwi-PiRNA that the Hs6st locus displays insertion polymorphism for helitron Pathway Provides an Adaptive Defense in the Transposon Arms Race. N1_Dbi, which is found approximately 1 kb away from tirant in Science 318: 761–764. https://doi.org/10.1126/science.1146484 Chicharo and Mayotte but not in Makindu. Apart from this heli- Biémont, C., and C. Vieira, 2006 Genetics: Junk DNA as an Evolutionary Force. Nature 443: 521–524. https://doi.org/10.1038/443521a tron and tirant, no other TE sequences are found within a 4 kb Boisvert, S., F. Laviolette, and J. Corbeil, 2010 Ray: Simultaneous Assembly window around tirant insertion site. To sum up, we still do not of Reads from a Mix of High-Throughput Sequencing Technologies. know why tirant affects tkv but not Hs6st. J. Comput. Mol. Cell Biol. 17: 1519–1533. Tirant effects are of small size: on average, a 30% reduction in Brennecke, J., A. A. Aravin, A. Stark, M. Dus, M. Kellis et al., 2007 Discrete tkv transcript amounts and an increase of 33–42% in H3K9me3 Small RNA-Generating Loci as Master Regulators of Transposon Activity enrichment. If it were not a small effect, we speculate that the in Drosophila. Cell 128: 1089–1103. https://doi.org/10.1016/j. insertion would be too deleterious and removed from the genome cell.2007.01.043 by natural selection. In addition, cases were also reported in which Camacho, C., G. Coulouris, V. Avagyan, N. Ma, J. Papadopoulos et al., TE insertions had only a modest impact on neighboring genes in 2009 BLAST+: Architecture and Applications. BMC Bioinformatics 10: standard conditions, but significant effects under stress conditions 421. https://doi.org/10.1186/1471-2105-10-421 Carver, T., U. Böhme, T. D. Otto, J. Parkhill, and M. Berriman, (Naito et al. 2009). 2010 BamView: Viewing Mapped Read Alignment Data in the Context of the Reference Sequence. Bioinformatics 26: 676–677. https://doi.org/ Conclusion 10.1093/bioinformatics/btq010 In the present study, we use natural variability regarding the tirant ERV Casacuberta, E., and J. González, 2013 The Impact of Transposable Ele- of D. simulans to identify potential impacts of TEs on the host ments in Environmental Adaptation. Mol. Ecol. 22: 1503–1517. https:// genome. We show that, to a moderate extent, a tirant insertion doi.org/10.1111/mec.12170 may spread heterochromatin to flanking regions where genes lie, Chambeyron, S., A. Popkova, G. Payen-Groschêne, C. Brun, D. Laouini et al., and be associated with transcript level reduction for the gene, with 2008 PiRNA-Mediated Nuclear Accumulation of Retrotransposon Transcripts in the Drosophila Female Germline. Proc. Natl. Acad. Sci. potential phenotypic impacts. Our results also illustrate that TE USA 105: 14964–14969. https://doi.org/10.1073/pnas.0805943105 effects on the genome may not be as pervasive as recently proposed Cordaux, R., and M. A. Batzer, 2009 The Impact of Retrotransposons on (Lee and Karpen 2017) since we detected them only for one inser- Human Genome Evolution. Nat. Rev. Genet. 10: 691–703. https://doi.org/ tion and not for the others. In addition, our data suggest the exis- 10.1038/nrg2640 tence of a dynamic relationship between TE transcripts and the Daborn, P. J., J. L. Yen, M. R. Bogwitz, G. Le Goff, E. Feil et al., 2002 A piRNAs that control them. This is illustrated by a positive correla- Single P450 Allele Associated with Insecticide Resistance in Drosophila. tion between these two variables in normal conditions, which is Science 297: 2253–2256. https://doi.org/10.1126/science.1074170 reverted to a negative correlation in the case of dysgenic crosses. Elliott, T. A., and T. R. Gregory, 2015 What’s in a Genome? The C–Value We speculate this may happen when asymmetry is too high between Enigma and the Evolution of Eukaryotic Genome Content. Philos. Trans. ping-pong partners abundances, –ie maternally transmitted piR- R. Soc. Lond. B Biol. Sci. 370: 20140331. https://doi.org/10.1098/ rstb.2014.0331 NAs on one hand and TE transcripts on the other hand, for the Fablet, M., E. Lerat, R. Rebollo, B. Horard, N. Burlet et al., 2009 Genomic piRNA pathway to be efficient. Environment Influences the Dynamics of the Tirant LTR Retrotransposon Taken together, our results participate in the understanding of the in Drosophila. FASEB J. Off. Publ. Fed. Am. Soc. Exp. Biol. 23: 1482–1489. equilibrium between the host genome and its TEs. This study opens the Fablet, M., J. F. McDonald, C. Biémont, and C. Vieira, 2006 Ongoing Loss way to the investigation of the transition from the equilibrium disruption of the Tirant Transposable Element in Natural Populations of Drosophila to a restored balance between the genome and its TEs, which are Simulans. Gene 375: 54–62. https://doi.org/10.1016/j.gene.2006.02.008 fundamental aspects of the understanding of TE biology and genome Fablet, M., J. Salces-Ortiz, B. Fraga Menezes, M. Roy, and C. Vieira, stability. 2017 Self and nonself from a genomic perspective: transposable ele- ments, pp. 111–128 in Evolutionary Biology: Self/Nonself Evolution, Spe- ACKNOWLEDGMENTS cies and Complex Traits Evolution, Methods and Concepts, edited by This work was performed using the computing facilities of the CC Springer, Cham. LBBE/PRABI, the luminometer of the IUT Lyon 1 Biology Fablet, M., S. Souames, C. Biémont, and C. Vieira, 2007 Evolutionary Pathways of the Tirant LTR Retrotransposon in the Drosophila Mela- Department, and the Vienna Biocenter Campus Facility. nogaster Subgroup of Species. J. Mol. Evol. 64: 438–447. https://doi.org/ We thank Justine Picarle, Sonia Martinez, Nabila Berabez and Anaïs 10.1007/s00239-006-0108-9 Alquezar for technical assistance. We thank Philippe Veber, Matthieu Garrell, J., and J. Modolell, 1990 The Drosophila extramacrochaetae locus, Boulesteix, and attendees of the XXth French National Congress on an antagonist of proneural genes that, like these genes, encodes a helix- Transposable Elements for fruitful discussions. We thank Labex loop-helix protein. Cell 61: 39–48. https://doi.org/10.1016/0092- ECOFECT. 8674(90)90213-X 864 | M. Fablet et al. Garrido-Martín, D., E. Palumbo, R. Guigó, and A. Breschi, 2018 ggsashimi: Modolo, L., 2018 l-modolo/fastq_sampler: everything is working. Zenodo. Sashimi plot revised for browser- and annotation-independent splicing https://doi.org/10.5281/zenodo.1227899 visualization. PLOS Comput. Biol. 14: e1006360. https://doi.org/10.1371/ Mohn, F., G. Sienski, D. Handler, and J. Brennecke, 2014 The Rhino- journal.pcbi.1006360 Deadlock-Cutoff Complex Licenses Noncanonical Transcription of Dual- Grabherr, M. G., B. J. Haas, M. Yassour, J. Z. Levin, D. A. Thompson et al., Strand PiRNA Clusters in Drosophila. Cell 157: 1364–1379. https://doi. 2011 Full-Length Transcriptome Assembly from RNA-Seq Data with- org/10.1016/j.cell.2014.04.031 out a Reference Genome. Nat. Biotechnol. 29: 644–652. https://doi.org/ Naito, K., F. Zhang, T. Tsukiyama, H. Saito, C. N. Hancock et al., 10.1038/nbt.1883 2009 Unexpected Consequences of a Sudden and Massive Transposon Hancks, D. C., and H. H. Kazazian, 2016 Roles for Retrotransposon In- Amplification on Rice Gene Expression. Nature 461: 1130–1134. https:// sertions in Human Disease. Mob. DNA 7: 9. https://doi.org/10.1186/ doi.org/10.1038/nature08479 s13100-016-0065-9 Rebollo, R., M. M. Karimi, M. Bilenky, L. Gagnier, K. Miceli-Royer et al., Jangam, D., C. Feschotte, and E. Betrán, 2017 Transposable Element Do- 2011 Retrotransposon-Induced Heterochromatin Spreading in the mestication As an Adaptation to Evolutionary Conflicts. Trends Genet. Mouse Revealed by Insertional Polymorphisms. PLoS Genet. 7: e1002301. TIG 33: 817–831. https://doi.org/10.1016/j.tig.2017.07.011 https://doi.org/10.1371/journal.pgen.1002301 Johnson, K. P., K. K. O. Walden, and H. M. Robertson, 2013 Next-Gen- Rey, C., T. Lorin, M. Sémon, and B. Boussau, 2017 apytram v1.1. Zenodo. eration Phylogenomics Using a Target Restricted Assembly Method. Mol. https://doi.org/10.5281/zenodo.804416 Phylogenet. Evol. 66: 417–422. https://doi.org/10.1016/j. Saito, K., K. M. Nishida, T. Mori, Y. Kawamura, K. Miyoshi et al., ympev.2012.09.007 2006 Specific Association of Piwi with RasiRNAs Derived from Retro- Kelleher, E. S., and D. A. Barbash, 2013 Analysis of PiRNA-Mediated Si- transposon and Heterochromatic Regions in the Drosophila Genome. lencing of Active TEs in Drosophila Melanogaster Suggests Limits on the Genes Dev. 20: 2214–2222. https://doi.org/10.1101/gad.1454806 Evolution of Host Genome Defense. Mol. Biol. Evol. 30: 1816–1829. Schmieder, R., and R. Edwards, 2011 Quality Control and Preprocessing of https://doi.org/10.1093/molbev/mst081 Metagenomic Datasets. Bioinformatics 27: 863–864. https://doi.org/ Kent, W. J., 2002 BLAT–the BLAST-like Alignment Tool. Genome Res. 12: 10.1093/bioinformatics/btr026 656–664. https://doi.org/10.1101/gr.229202 Senti, K.-A., D. Jurczak, R. Sachidanandam, and J. Brennecke, 2015 piRNA- Kim, D., G. Pertea, C. Trapnell, H. Pimentel, R. Kelley et al., 2013 TopHat2: guided slicing of transposon transcripts enforces their transcriptional accurate alignment of transcriptomes in the presence of insertions, de- silencing via specifying the nuclear piRNA repertoire. Genes Dev. 29: letions and gene fusions. Genome Biol. 14: R36. https://doi.org/10.1186/ 1747–1762. https://doi.org/10.1101/gad.267252.115 gb-2013-14-4-r36 Shpiz, S., S. Ryazansky, I. Olovnikov, Y. Abramov, and A. Kalmykova, Langmead, B., and S. L. Salzberg, 2012 Fast Gapped-Read Alignment with 2014 Euchromatic Transposon Insertions Trigger Production of Novel Bowtie 2. Nat. Methods 9: 357–359. https://doi.org/10.1038/nmeth.1923 Pi- and Endo-SiRNAs at the Target Sites in the Drosophila Germline. Langmead, B., C. Trapnell, M. Pop, and S. L. Salzberg, 2009 Ultrafast and PLoS Genet. 10: e1004138. https://doi.org/10.1371/journal.pgen.1004138 Memory-Efficient Alignment of Short DNA Sequences to the Human Sienski, G., D. Dönertas, and J. Brennecke, 2012 Transcriptional Silencing of Genome. Genome Biol. 10: R25. https://doi.org/10.1186/gb-2009-10-3-r25 Transposons by Piwi and Maelstrom and Its Impact on Chromatin State and Lassmann, T., Y. Hayashizaki, and C. O. Daub, 2011 SAMStat: Monitoring Gene Expression. Cell 151: 964–980. https://doi.org/10.1016/j.cell.2012.10.040 Biases in next Generation Sequencing Data. Bioinformatics 27: 130–131. Simão, F. A., R. M. Waterhouse, P. Ioannidis, E. V. Kriventseva, and https://doi.org/10.1093/bioinformatics/btq614 E. M. Zdobnov, 2015 BUSCO: Assessing Genome Assembly and An- Le Thomas, A., A. K. Rogers, A. Webster, G. K. Marinov, S. E. Liao et al., notation Completeness with Single-Copy Orthologs. Bioinformatics 31: 2013 Piwi Induces PiRNA-Guided Transcriptional Silencing and Es- 3210–3212. https://doi.org/10.1093/bioinformatics/btv351 tablishment of a Repressive Chromatin State. Genes Dev. 27: 390–399. Siomi, M. C., K. Sato, D. Pezic, and A. A. Aravin, 2011 PIWI-Interacting https://doi.org/10.1101/gad.209841.112 Small RNAs: The Vanguard of Genome Defence. Nat. Rev. Mol. Cell Biol. Lee, Y. C. G., and G. H. Karpen, 2017 Pervasive epigenetic effects of Drosophila 12: 246–258. https://doi.org/10.1038/nrm3089 Euchromatic transposable elements impact their evolution. eLife 6: e25762. Slater, G. S., and E. Birney, 2005 Automated generation of heuristics for Lerat, E., M. Fablet, L. Modolo, H. Lopez-Maestre, and C. Vieira, biological sequence comparison. BMC Bioinformatics 6: 31. https://doi. 2017 TEtools Facilitates Big Data Expression Analysis of Transposable org/10.1186/1471-2105-6-31 Elements and Reveals an Antagonism between Their Activity and That of Slotkin, R. K., and R. Martienssen, 2007 Transposable Elements and the PiRNA Genes. Nucleic Acids Res. 45: e17. Epigenetic Regulation of the Genome. Nat. Rev. Genet. 8: 272–285. Li, H., B. Handsaker, A. Wysoker, T. Fennell, J. Ruan et al., 2009 The https://doi.org/10.1038/nrg2072 Sequence Alignment/Map Format and SAMtools. Bioinformatics 25: Smit, A. F., R. Hubley, and P. Green, 2013 RepeatMasker Open-4.0. http:// 2078–2079. https://doi.org/10.1093/bioinformatics/btp352 www.repeatmasker.org. Lindsley, D. L., and G. G. Zimm, 1992 The Genome of Drosophila Mela- Teixeira, F. K., M. Okuniewska, C. D. Malone, R.-X. Coux, D. C. Rio et al., nogaster, Academic Press, Inc., San Diego 2017 PiRNA-Mediated Regulation of Transposon Alternative Splicing Love, M. I., W. Huber, and S. Anders, 2014 Moderated Estimation of Fold in the Soma and Germ Line. Nature 552: 268–272. Change and Dispersion for RNA-Seq Data with DESeq2. Genome Biol. Tenaillon, M. I., J. D. Hollister, and B. S. Gaut, 2010 A Triptych of the 15: 550. https://doi.org/10.1186/s13059-014-0550-8 Evolution of Plant Transposable Elements. Trends Plant Sci. 15: 471–478. Marsano, R. M., R. Moschetti, C. Caggese, C. Lanave, P. Barsanti et al., https://doi.org/10.1016/j.tplants.2010.05.003 2000 The Complete Tirant Transposable Element in Drosophila Mela- Terzian, C., A. Pélisson, and A. Bucheton, 2001 Evolution and Phylogeny of nogaster Shows a Structural Relationship with Retrovirus-like Retrotrans- Insect Endogenous Retroviruses. BMC Evol. Biol. 1: 3. https://doi.org/ posons. Gene 247: 87–95. https://doi.org/10.1016/S0378-1119(00)00115-3 10.1186/1471-2148-1-3 Martin, M., 2011 Cutadapt Removes Adapter Sequences from High- Vieira, C., D. Lepetit, S. Dumont, and C. Biémont, 1999 Wake up of Throughput Sequencing Reads. EMBnet.journal 17: 10–12. https://doi.org/ Transposable Elements Following Drosophila Simulans Worldwide Col- 10.14806/ej.17.1.200 onization. Mol. Biol. Evol. 16: 1251–1255. https://doi.org/10.1093/ox- Mi, S., X. Lee, X. Li, G. M. Veldman, H. Finnerty et al., 2000 Syncytin Is a fordjournals.molbev.a026215 Captive Retroviral Envelope Protein Involved in Human Placental Mor- Wicker, T., F. Sabot, A. Hua-Van, J. L. Bennetzen, P. Capy et al., 2007 A phogenesis. Nature 403: 785–789. https://doi.org/10.1038/35001608 Unified Classification System for Eukaryotic Transposable Elements. Nat. Modolo, L., and E. Lerat, 2015 UrQt: An Efficient Software for the Un- Rev. Genet. 8: 973–982. https://doi.org/10.1038/nrg2165 supervised Quality Trimming of NGS Data. BMC Bioinformatics 16: 137. https://doi.org/10.1186/s12859-015-0546-8 Communicating Editor: J. Birchler Volume 9 March 2019 | Tirant ERV Dynamics in Drosophila | 865