Species-specific chromatin landscape determines how transposable elements shape genome evolution

  1. Yuheng Huang
  2. Harsh Shukla
  3. Yuh Chwen G Lee  Is a corresponding author
  1. Department of Ecology and Evolutionary Biology, University of California, Irvine, United States

Abstract

Transposable elements (TEs) are selfish genetic parasites that increase their copy number at the expense of host fitness. The ‘success’, or genome-wide abundance, of TEs differs widely between species. Deciphering the causes for this large variety in TE abundance has remained a central question in evolutionary genomics. We previously proposed that species-specific TE abundance could be driven by the inadvertent consequences of host-direct epigenetic silencing of TEs—the spreading of repressive epigenetic marks from silenced TEs into adjacent sequences. Here, we compared this TE-mediated local enrichment of repressive marks, or ‘the epigenetic effect of TEs’, in six species in the Drosophila melanogaster subgroup to dissect step-by-step the role of such effect in determining genomic TE abundance. We found that TE-mediated local enrichment of repressive marks is prevalent and substantially varies across and even within species. While this TE-mediated effect alters the epigenetic states of adjacent genes, we surprisingly discovered that the transcription of neighboring genes could reciprocally impact this spreading. Importantly, our multi-species analysis provides the power and appropriate phylogenetic resolution to connect species-specific host chromatin regulation, TE-mediated epigenetic effects, the strength of natural selection against TEs, and genomic TE abundance unique to individual species. Our findings point toward the importance of host chromatin landscapes in shaping genome evolution through the epigenetic effects of a selfish genetic parasite.

Editor's evaluation

Transposable elements are genomic parasites and the fraction of the genome that is made up of such elements varies greatly between species, and models suggest that this must reflect the balance between the rate at which they multiply, and the rate at which selection purges them from the genome. Precisely how selection acts against transposable element insertions is not clear. This paper provides evidence that the strength of selection depends on the extent to which epigenetic silencing spreads to nearby genes – although the mechanism is obscure, as gene expression is not affected. This is a very interesting hypothesis that deserves more attention, and the paper is an excellent example of trying to combine population genetics models with a mechanistic understanding of the process modeled.

https://doi.org/10.7554/eLife.81567.sa0

eLife digest

All the instructions required for life are encoded in the set of DNA present in a cell. It therefore seems natural to think that every bit of this genetic information should serve the organism. And yet most species carry parasitic ‘transposable’ sequences, or transposons, whose only purpose is to multiply and insert themselves at other positions in the genome.

It is possible for cells to suppress these selfish elements. Chemical marks can be deposited onto the DNA to temporarily ‘silence’ transposons and prevent them from being able to move and replicate. However, this sometimes comes at a cost: the repressive chemical modifications can spread to nearby genes that are essential for the organism and perturb their function.

Strangely, the prevalence of transposons varies widely across the tree of life. These sequences form the majority of the genome of certain species – in fact, they represent about half of the human genetic information. But their abundance is much lower in other organisms, forming a measly 6% of the genome of puffer fish for instance. Even amongst fruit fly species, the prevalence of transposable elements can range between 2% and 25%. What explains such differences?

Huang et al. set out to examine this question through the lens of transposon silencing, systematically comparing how this process impacts nearby regions in six species of fruit flies. This revealed variations in the strength of the side effects associated with transposon silencing, resulting in different levels of perturbation on neighbouring genes. A stronger impact was associated with the species having fewer transposons in its genome, suggesting that an evolutionary pressure is at work to keep the abundance of transposons at a low level in these species. Further analyses showed that the genes which determine how silencing marks are distributed may also be responsible for the variations in the impact of transposon silencing. They could therefore be the ones driving differences in the abundance of transposons between species.

Overall, this work sheds light on the complex mechanisms shaping the evolution of genomes, and it may help to better understand how transposons are linked to processes such as aging and cancer.

Introduction

Transposable elements (TEs) are widespread genetic parasites that copy and insert themselves across host genomes. The presence and movement of TEs could impair host genome functions. TEs disrupt genes and functional elements (Finnegan, 1992), introduce ectopic regulatory sequences (Chuong et al., 2017), and trigger highly deleterious chromosomal rearrangements through nonhomologous recombination (Langley et al., 1988; Montgomery et al., 1991). Nevertheless, the ability to self-replicate has allowed TEs to successfully occupy nearly all eukaryotic genomes surveyed (reviewed in Wells and Feschotte, 2020). Within a eukaryotic genome, TEs are prevalent in both gene-poor, repeat-rich heterochromatic and gene-rich euchromatic regions (e.g., Kaminker et al., 2002; Bergman et al., 2006). TEs in the heterochromatic genome are oftentimes fragmented, losing their ability to replicate (e.g., Hoskins et al., 2007; Hoskins et al., 2015). On the contrary, euchromatic TEs, which could intersperse with functional elements, commonly retain the potential to replicate, making them active players for not only their own evolutionary dynamics, but also the function and evolution of the euchromatic genome. Because of that and the technical challenges associated with identifying TEs in the heterochromatic regions (Salzberg and Yorke, 2005; Treangen and Salzberg, 2011), studies have been largely focused on the evolution of TEs in the euchromatic genome. Intriguingly, the abundance of TEs in the euchromatic genome substantially varies across the phylogenetic tree (Huang et al., 2012; Elliott and Gregory, 2015; Wells and Feschotte, 2020). For instance, in the assembled vertebrate genomes, which mostly consist of euchromatic regions, the proportion occupied by TEs ranges from only 6% in pufferfish (Volff et al., 2003) to more than 65% in salamander (Nowoshilow et al., 2018). Even within the same genus, genomic TE abundance differs widely (e.g., 2.5–25% of assembled genome sequences in Drosophila; Clark et al., 2007; Rius et al., 2016). Deciphering the role of this prevalent parasite in shaping genome evolution has remained a central question in genomics (Kazazian, 2004; Feschotte and Pritham, 2007; Arkhipova and Kumar, 2018); however, the ultimate causes of such dramatic divergence in TE abundance in the euchromatic genome remain unclear.

Theoretical analyses proposed that, in panmictic host populations with unrestricted recombination, TE abundance is determined by how quickly TEs replicate and how fast they are removed from the populations by natural selection against their harmful fitness effects (Charlesworth and Charlesworth, 1983, reviewed in Lee and Langley, 2010). Under this model, divergent genome-wide TE abundance could be driven by between-species differences in the strength of selection against TEs. Currently available evolutionary models that address this possibility have focused on population genetic parameters that influence the efficacy of selection removing TEs, such as mating systems (Wright and Schoen, 1999; Dolgin and Charlesworth, 2006; Boutin et al., 2012) and effective population size (Lynch and Conery, 2003). Yet, empirical support for such a hypothesis has been mixed (Dolgin et al., 2008; Lockton and Gaut, 2010; de la Chaux et al., 2012; Arunkumar et al., 2015; Agren et al., 2014; Mérel et al., 2021; Oggenfuss et al., 2021). On the other hand, between-species differences in the magnitude of harmful effects exerted by TEs, and accordingly the strength of selection against TEs, could also determine genomic TE abundance, a plausible hypothesis that is yet to have empirical investigations.

A new avenue for exploring how these genetic parasites shape the function and evolution of eukaryotic genomes was opened by the recently discovered host-directed silencing of TEs and the associated ‘inadvertent’ deleterious epigenetic effects (reviewed in Choi and Lee, 2020). To counteract the selfish increase of TEs in host genomes, eukaryotic hosts have evolved small RNA-mediated mechanisms to transcriptionally silence TEs (reviewed in Slotkin and Martienssen, 2007; Czech et al., 2018; Deniz et al., 2019). Host protein complexes are guided by small RNAs to TEs with complementary sequences, which is followed by the recruitment of methyltransferases that modify DNA or histone tails at TE sequences (Qi et al., 2006; Aravin et al., 2008; Wang and Elgin, 2011; Sienski et al., 2012; Le Thomas et al., 2013). Such a process results in the enrichment of DNA methylation or di- and tri-methylation on lysine 9 of histone H3 (H3K9me2/3), both repressive epigenetic modifications that are typically found in heterochromatic regions and associated with repressed gene expression (reviewed in Pikaard and Mittelsten Scheid, 2014; Allis and Jenuwein, 2016). This repressed transcription of TEs results in reduced RNA intermediates (for RNA-based TEs) and proteins (e.g., transposase and reverse transcriptase) necessary for TE replication, effectively slowing the selfish propagation of TEs.

While such epigenetic silencing of TEs should benefit their hosts, studies in various model species have found that repressive marks enriched at silenced TEs ‘spread’ beyond TE boundaries, leading to local enrichment of such marks at TE-adjacent sequences across the euchromatic genomes (i.e., Mus, Drosophila, Arabidopsis, and Oryza; Rebollo et al., 2012; Sienski et al., 2012; Pezic et al., 2014; Lee, 2015; Quadrana et al., 2016; Choi and Purugganan, 2018, reviewed in Choi and Lee, 2020). Furthermore, TEs with such effects were observed to have lower population frequencies (Hollister and Gaut, 2009; Lee, 2015; Lee and Karpen, 2017), suggesting that selection acts to remove them. These discoveries highlight the potential importance of TE-triggered epigenetic effects in shaping genome evolution. Interestingly, the strength of TE-mediated local enrichment of repressive marks substantially varies between distantly related taxa (reviewed in Choi and Lee, 2020). Investigations on pairs of closely related species further revealed that this ‘epigenetic effect of TE’ differs and is stronger in the species with fewer euchromatic TEs (Arabidopsis thaliana vs. Arabidopsis lyrata and Drosophila melanogaster vs. Drosophila simulans; Hollister et al., 2011; Lee and Karpen, 2017). These observations spurred our previous hypothesis that, across species, different TE-mediated enrichment of repressive marks could result in varying functional consequences and thus differences in the strength of selection against TEs, eventually contributing to divergent TE abundance in the euchromatic genome (Lee and Karpen, 2017). We further postulated that this difference in TE-mediated epigenetic effects could have resulted from species-specific genetic modulation of the repressive chromatin landscape (Lee and Karpen, 2017), which was shown to determine the spreading of repressive epigenetic marks from constitutive heterochromatin into adjacent euchromatic sequences (reviewed in Girton and Johansen, 2008; Elgin and Reuter, 2013).

To fully examine the hypothesis that varying host chromatin landscape drives between-species differences of TEs through epigenetic mechanisms, one needs to connect species-specific regulation of chromatin landscape, TE-mediated enrichment of repressive marks, the associated functional consequence and resultant selection against TEs, and genomic TE abundance. However, former analyses that compared TE-mediated epigenetic effects between species have limited sampling (two species) and thus lack sufficient statistical power for robust inference (Hollister et al., 2011; Lee and Karpen, 2017). Also, support for key links of the hypothesis is lacking. For instance, selection against TE-mediated epigenetic effects was expected to result from the associated reducing effects on the expression of neighboring genes. Yet, investigations in multiple taxa reported weak or no associations between the epigenetic effects of TEs and neighboring gene expression (Quadrana et al., 2016; Stuart et al., 2016; Lee and Karpen, 2017; Choi and Purugganan, 2018, reviewed in Kelleher et al., 2020; Choi and Lee, 2020). These inconclusive observations cast doubt on the possibility that TE-mediated epigenetic effects impair host fitness by silencing neighboring genes and whether this particular deleterious consequence indeed shapes genome evolution. Furthermore, previous comparisons of population frequencies between TEs with and without epigenetic effects, an approach used to infer the strength of natural selection removing TEs, could not exclude the confounding influence of other harmful effects of TEs on their population frequencies (e.g., Hollister and Gaut, 2009; Lee, 2015). Accordingly, those analyses could not unequivocally support selection against TE-mediated epigenetic effects. Multi-species studies that span an appropriate evolutionary distance and connect the missing links in the proposed hypothesis would be needed to test the predicted importance of TE-mediated epigenetic effects in determining between-species differences of TEs.

In this study, we investigated the prevalence of TE-mediated local enrichment of repressive epigenetic marks, or ‘TE-mediated epigenetic effects’, in the euchromatic genome of six species in the D. melanogaster subgroup (diverged around 10 MYR; Obbard et al., 2012, Figure 1A). These species are from the two well-studied species complexes (melanogaster and yakuba complexes), providing good phylogenetic resolution to address the role of TE-mediated epigenetic effects in genome evolution. While TE insertions in all species studied result in robust local enrichment of repressive epigenetic marks, the strength of such effects varies substantially within genomes, among species, and between species complexes. Our larger sample size allowed us to re-examine the still debated question about the impacts of TE-mediated enrichment of repressive marks on neighboring gene expression, which surprisingly revealed their complex interactions. Importantly, our multi-species analysis provides the power to test the predicted associations between TE abundance in the euchromatic genome, TE-mediated epigenetic effects, selection against such effects, and host chromatin environment, while uncovering the evolutionary causes for the wide variety of TE abundance between species.

Figure 1 with 8 supplements see all
Enrichment of H3K9me2 around euchromatic transposable elements (TEs) across species.

(A) Phylogenetic relationship among species included in this study. Species in the Drosophila melanogaster complex are in pink, while those in Drosophila yakuba complex are in blue. Numbers after each species denote the number of euchromatic TEs called by RepeatModeler2, before assignment into TE families and merging of adjacent copies (see Materials and methods). (B) Genome-wide average H3K9me2 HMD levels around euchromatic TEs with LOESS smoothing (span = 15%) in studied genomes; 95% confidence interval around smooth is shown as shaded areas. (C) Genome-wide average H3K9me2 HMD levels at homologous sequences in the presence (orange) and absence (gray) of euchromatic TEs in species that have data for two strains (Drosophila simulans and D. yakuba). The average H3K9me2 HMD level was smoothed with LOESS (span = 15%) with 95% confidence intervals around smooth shown as the shaded areas. (D) An Integrated-Genome-Viewer view showing the local enrichment of H3K9me2 around a hobo TE, and the two estimates (magnitude and extent of H3K9me2 enrichment) estimated for quantifying the epigenetic effects of individual TE.

Results

TE-mediated local enrichment of repressive marks is prevalent across studied species

We investigated whether previously reported local enrichment of heterochromatic marks around euchromatic TEs in model Drosophila species (Lee, 2015; Lee and Karpen, 2017) is prevalent across species in the D. melanogaster subgroup. H3K9me2 and H3K9me3 are histone modifications that are highly enriched in the constitutive heterochromatin in D. melanogaster (Riddle et al., 2011; Kharchenko et al., 2011) and are generally considered ‘heterochromatic marks’. Previously, it was shown that TEs lead to a local enrichment of both of these two histone modifications in the D. melanogaster euchromatic genome (Lee and Karpen, 2017), and we chose to focus on one of them (H3K9me2) in this study. We performed spike-in controlled chromatin immunoprecipitation (ChIP)-seq targeting H3K9me2 using 16–18 hr embryos (see Materials and methods). We estimated histone modification density (HMD), which is the ratio of fragment coverage in ChIP samples to that in matching input samples, standardized by the ratio of spike-in fragments (Lam et al., 2019) (see Materials and methods). TE insertions in strains used for the ChIP-seq experiment were annotated by running Repeatmodeler2 (Flynn et al., 2020) on genomes assembled using long-read PacBio sequencing. The high continuity of these genomes enables a more comprehensive identification of TEs than previous studies based on short-read sequencing data (Khost et al., 2017; Chakraborty et al., 2019). Because our analysis mainly focuses on the evolutionary dynamics of euchromatic TEs and their role in driving the evolution of the euchromatic genome, we excluded TEs in or near heterochromatic regions from our analysis (see Materials and methods).

Across all six species analyzed, we observed significant enrichment of H3K9me2 HMD around TEs, and this enrichment decreases to the background level within 10 kb (Figure 1B). To exclude the possibility that this enrichment of H3K9me2 is due to TE preferentially inserted into regions that are already enriched with heterochromatic marks (e.g., Dimitri and Junakovic, 1999), we collected H3K9me2 epigenomic data for two genomes of D. simulans and Drosophila yakuba and compared the enrichment of H3K9me2 for homologous sequences with and without a TE insertion. The H3K9me2 enrichment is observed in the vicinity of TEs in the genome where they are present, but not at homologous sequences in the other genome (Figure 1C). This observation expanded previous single-species observations (Lee, 2015; Lee and Karpen, 2017) and supported that the local enrichment of heterochromatic marks in the euchromatic regions is induced by TEs in multiple species.

Strength of TE-mediated epigenetic effects depends jointly on TE attributes and host genetic background

In order to investigate biological factors associated with the strength of TE-mediated epigenetic effects, we quantified the ‘magnitude’ and ‘extent’ of local enrichment of H3K9me2 for individual TEs (Figure 1D). We estimated the magnitude of TE-mediated epigenetic effects as the H3K9me2 enrichment in the 1 kb window immediately adjacent to a TE insertion. To identify the extent of TE-mediated enrichment of H3K9me2, we scanned from TE to locate the farthest 1 kb window in which H3K9me2 enrichment level is above that of the local background (see Materials and methods for details). For D. simulans and D. yakuba, we also estimated the magnitude and extent of TE-induced H3K9me2 enrichment by comparing H3K9me2 enrichment at homologous sequences between strains with and without focal TEs (two-genome estimates; see Materials and methods). Estimates based on one genome or two genomes strongly correlate (Spearman rank correlation coefficient ρ=0.64–0.85 [magnitude] and 0.40–0.64 [extent], p<10–10 for all tests, Figure 1—figure supplement 1). Because an important aspect of our analyses is the comparison of TE-induced H3K9me2 across species, we reported analyses based on single-genome estimates henceforth.

The assembly of constitutive heterochromatin has been proposed to depend on the concentration of heterochromatic enzymes and structural proteins, which should be the highest at the nucleation site and gradually decrease, leading to the cis ‘spreading’ of repressive marks (Locke et al., 1988). Under this model, the magnitude of the enrichment for repressive marks at the nucleation site should determine the extent of the enrichment for repressive marks. Yet, predictions of this model were found to be inconsistent with several empirical observations, including the discontinuous ‘spreading’ of repressive marks from constitutive heterochromatin (Belyaeva and Zhimulev, 1991; Talbert and Henikoff, 2000) and the dependency of the extent of such effect on factors other than heterochromatin mass (Sabl and Henikoff, 1996). If similar molecular mechanisms are also applicable to epigenetically silenced TE in the euchromatic genome, the magnitude and the extent of TE-mediated local enrichment of H3K9me2 would not perfectly correlate and should capture different aspects of such TE-mediated effects. While we have no a priori predictions for the relative importance of these two indexes in the questions that this study aims to address, we anticipate that the extent of TE-mediated H3K9me2 enrichment should more likely be influenced by local genomic context than the magnitude of such effect. Nevertheless, we found that the magnitude and extent of TE-mediated H3K9me2 enrichment strongly correlated within genomes (Spearman rank correlation coefficient ρ=0.58–0.75, p<10–16, Figure 1—figure supplement 2).

TE length was postulated to be an important factor determining the strength of TE-mediated epigenetic effects because silenced TEs that are longer in length are expected to represent larger heterochromatin mass (Lee, 2015). Consistent with the prediction, we observed significant, though weak, positive correlations between TE length and the strength of TEs’ epigenetic effects within most genomes studied (Spearman rank correlation coefficient ρ=0.12–0.22, p<0.05; Figure 1—figure supplement 3). It is worth noting that previous analysis on non-reference D. melanogaster strains was unable to test this prediction, due to the inability to assemble internal sequences of TEs with short-read resequencing data (Lee and Karpen, 2017). For the same reason, previous analysis could not study whether the epigenetic effects differed between full-length and truncated TEs that have different potential to be transcribed. This distinction between TE insertions could be important because TE-mediated local enrichment of repressive marks was mainly observed with transcriptionally active TEs (Pezic et al., 2014, but see Sentmanat and Elgin, 2012). Consistently, we found that full-length TEs exert significantly larger magnitude and extent of H3K9me2 enrichment than truncated TEs in several genomes (Mann-Whitney U test, p<0.05 for D. melanogaster, Drosophila mauritiana, Drosophila santomea (magnitude) and for D. melanogaster and D. santomea (extent); Figure 1—figure supplement 4).

We next compared TEs of different classes, which are classifications based on the transposition mechanisms of TEs (Wicker et al., 2007) and previously observed to associate with varying strength of epigenetic effects within D. melanogaster genomes (Lee, 2015; Lee and Karpen, 2017). While different classes of TEs showed similar levels of H3K9me2 enrichment, there is a very strong species-complex effect in which TEs in genomes of yakuba complex showed much larger magnitude of epigenetic effects than those in genomes of melanogaster complex (Figure 2A, an average 1.7-fold larger). On the other hand, the extent of H3K9me2 enrichment is more variable and does not show a similar trend (Figure 2A). Interestingly, the extent and magnitude of H3K9me2 spreading vary substantially between TEs of different families within a class and between TEs of the same family (Figure 2B for D. simulans and see Figure 2—figure supplement 1 for other melanogaster complex species and Figure 2—figure supplement 2 for yakuba complex species). Moreover, the rank order of the extent and magnitude of TE-induced H3K9me2 enrichment of TE families varies between species and even between strains of the same species. These observations strongly suggest that the strength of TE-mediated epigenetic effects depends on both TE family attributes and host genetic background. It is worth noting that the percentage of TEs with a family assigned is higher in the melanogaster complex than that in the yakuba complex (Figure 2—figure supplement 3), and our analysis likely missed TE families that are highly divergent in and/or unique to the yakuba complex species.

Figure 2 with 4 supplements see all
Variation in the epigenetic effects of transposable elements (TEs) within genomes.

(A) The mean magnitude (left) and extent (right) of TE-induced H3K9me2 enrichment for different types of TEs in eight genomes from six species are shown. Different colors represent TEs of different classes, including Terminal Inverted Repeat (TIR), Long Terminal Repeat (LTR), and non-Long Terminal Repeat (non-LTR, also known as LINE) insertions. (B) The magnitude (left) and extent (right) of TE-induced H3K9me2 enrichment for different TE families in the two strains of Drosophila simulans. Only TE families with at least five identified copies in a genome were included. See Figure 2—figure supplements 2 and 3 for other genomes. (C) The median magnitude (left) and extent (right) of TE-induced H3K9me2 enrichment for six TE families with at least five copies in all genomes studied.

To further study the effects of host genome-by-family interaction across species, we compared the epigenetic effects of TEs from families that are found in all strains and have at least five copies (Figure 2C). Interestingly, while some TE families show universally strong H3K9me2 enrichment (e.g., 1360 for the magnitude) or weak H3K9me2 enrichment (e.g., FB for the magnitude) across species, some families’ epigenetic effects clearly depend on host genotype (e.g., H element shows larger magnitude and extent of H3K9me2 enrichment in yakuba complex species than in melanogaster complex species). Using a linear regression model, we found significant family-by-strain interaction effects for both the magnitude (ANOVA F-value: 3.4, df = 35, p=1.4 × 10–10) and the extent (ANOVA F-value: 4.4, df = 35, p<5.8 × 10–16) of H3K9me2 enrichment. While the significant host genome-by-family interaction on the magnitude of H3K9me2 enrichment could be strongly driven by Cr1a family, all TE families investigated showed strong host genome-by-family interactions for the extent of H3K9me2 enrichment (Figure 2C).

In addition to the intrinsic biological attributes of individual TEs, the insertion locations of TEs may also influence the magnitude and extent of their epigenetic effects, especially given that the chromatin environment is quite different between genic and non-genic sequences (Filion et al., 2010; Kharchenko et al., 2011). We categorized TEs according to their insertion locations relative to genes (intergenic, intronic, and exonic), but did not find a significant difference between them in terms of the magnitude of H3K9me2 enrichment except in a few cases (intergenic > exonic (Drosophila teissieri and D. santomea); exonic > intergenic ~ intronic (D. simulans strain 1); Mann-Whitney U test, p<0.05; Figure 2—figure supplement 4). On the other hand, for the extent of TE-mediated H3K9me2 enrichment, intergenic TEs exert larger such effects than other TEs in multiple genomes (Mann-Whitney U test, p<0.05 for D. simulans strains 1 and 2, D. santomea, and D. yakuba strain 2; Figure 2—figure supplement 4). Together, our observations revealed that the epigenetic effects of TEs jointly depend on the class, family identity, length, and insertion locations of TEs as well as the host genetic background (also see below).

Complex relationship between TE-induced enrichment of H3K9me2 and neighboring gene expression

Euchromatic TEs are interspersed with actively transcribing genes, and TE-mediated spreading of H3K9me2 was previously observed to extend into neighboring genes (Lee, 2015; Lee and Karpen, 2017). The enrichment of H3K9me2, a repressive histone modification, is generally associated with suppressed gene expression (Kouzarides, 2007). Accordingly, TE-mediated epigenetic effects were expected to lower neighboring gene expression. Yet, previous studies in several model species found limited effects of TE-induced enrichment of repressive marks on adjacent gene expression (reviewed in Kelleher et al., 2020; Choi and Lee, 2020). These observations left an important yet unsolved question about the functional importance of TE-mediated epigenetic effects.

To revisit this question with expanded data, we studied whether TE-mediated epigenetic effects associate with H3K9me2 enrichment or the expression level of their nearest adjacent genes. Because TEs inside genes could alter gene expression through other mechanisms and potentially confound the analysis (e.g., the disruption of functional sequences, see Kelleher et al., 2020; Choi and Lee, 2020), we only included intergenic TEs for analyses in this section. Similar to previous observations made in D. melanogaster (Lee, 2015; Lee and Karpen, 2017), H3K9me2 enrichment level at gene body positively correlates with the magnitude and extent of TE-induced H3K9me2 enrichment in all species analyzed (Spearman rank correlation coefficient ρ=0.16–0.33 [magnitude] and 0.18–0.31 [extent], p<0.05 for all tests, Figure 3A). Importantly, for the magnitude of TE-mediated epigenetic effects, this association is much stronger for genes close to TEs than those far from TEs (regression analysis: genic H3K9me2~TE epigenetic effects + close/far from genes + interaction, p-value for interaction term <0.05 for all genomes except for D. mauritiana [p=0.06] and D. teissieri [p=0.09]; Figure 3—figure supplement 1). For the extent of TE-mediated epigenetic effects, we observed similar, but weaker, distance-dependent effects (regression analysis, p-value for interaction term <0.001 for D. mauritiana, D. simulans strain 1, D. teissieri, and p=0.06 for D. melanogaster; Figure 3—figure supplement 2). Regression analysis without binning genes into those close or distant to TEs reached a similar conclusion that the associations between TE-mediated epigenetic effects and genic H3K9me2 enrichment depend on gene-TE distance (genic H3K9me2 ~ TE epigenetic effects + distance to TE + interaction; p-value for interaction term <0.05 for all tests, except for the extent of the effect for D. yakuba strain 1, p=0.05).

Figure 3 with 13 supplements see all
Associations between transposable element (TE)-mediated H3K9me2 enrichment and the epigenetic states and expression of neighboring genes.

(A) Spearman rank correlation coefficients (ρ) between the magnitude (filled bars) and extent (open bars) of TE-mediated H3K9me2 enrichment and genic H3K9me2 enrichment level (top) and gene expression rank (higher rank means lower expression; bottom) of nearby genes. Most of the Spearman rank correlation coefficients are significantly different from 0 for comparisons of genic H3K9me2 enrichment level (top), but not for comparisons of gene expression rank (bottom). (B) Spearman rank correlation coefficients (ρ) between the magnitude (top) and extent (bottom) of TE-mediated H3K9me2 enrichment and the H3K9me2 enrichment level (left) and gene expression rank (right) of nearby genes for TEs 5’ (dark blue/orange) and 3’ (light blue/light orange) to genes. For the genic H3K9me2 enrichment, all Spearman rank correlation coefficients are significantly different from 0 except for one test. For gene expression rank, few correlations are significantly different from 0. (C) z-Scores for comparing the H3K9me2 enrichment (left) and expression rank (right) of homologous genic alleles whose nearby TEs with (blue/orange) or without (gray) epigenetic effects (as defined as the magnitude of H3K9me2 enrichment >1; see Figure 3—figure supplement 7 for categorizing TEs with the extent of H3K9me2 enrichment, which gives similar results). A positive z-score means the allele with TE has higher H3K9me2 enrichment or larger expression rank (i.e., lower expression level) than the homologous allele without TE in another strain. (D) A cartoon describing the ‘genic side’ and ‘intergenic side’ extent of H3K9me2 enrichment mediated by TEs is shown on the left. z-Scores for comparing the extent of TE-mediated H3K9me2 enrichment on the intergenic side and on the genic side for TEs close to (green) and far (gray) from genes whose expression is at least 10 RPKM. A positive z-score means that the extent of TE-mediated H3K9me2 enrichment is more restricted on the genic side than the intergenic side. In several genomes, z-scores for TEs close to genes are significantly different from 0 and/or larger than those for TEs distant to genes. mel: Drosophila melanogaster, mau: Drosophila mauritiana, sim1: Drosophila simulans strain 1, sim2: D. simulans strain 2, tei: Drosophila teissieri, san: Drosophila santomea, yak1: Drosophila yakuba strain 1, yak2: D. yakuba strain 2. ***p<0.001, **p<0.01, *p<0.05 for Spearman rank correlation tests (A, B) and Mann-Whitney U tests (C, D).

Surprisingly, we found nearly absent correlations between TE-mediated H3K9me2 enrichment and the expression rank (ranking from the highest expressed genes) of adjacent genes in all genomes, except in one incidence (Spearman rank correlation tests, p>0.05 for all tests except for the extent of the effect of D. yakuba strain 2, Spearman rank correlation coefficient ρ=0.17, p<0.01; Figure 3A). Furthermore, the associations between gene expression and the epigenetic effects of TEs do not differ between close and distant gene-TE pairs for most genomes (gene expression ~ TE epigenetic effects + close/far from genes + interaction, p-value for interaction term >0.05 for all genomes except for D. simulans strain 1 (magnitude) and D. santomea (extent); Figure 3—figure supplements 3 and 4 for the magnitude and extent of TE-mediated epigenetic effects, respectively) nor depend on gene-TE distance except in one genome (gene expression ~ TE epigenetic effects + distance to TE + interaction; p-value for interaction term >0.05 for all genomes except for D. santomea, p<0.05 for both magnitude and extent of the effects). These observations substantially differ from our observed prevalent associations between TE-mediated effects and genic epigenetic states across species (Figure 3A, Figure 3—figure supplements 1 and 2).

TEs upstream to genes (i.e., 5’ to genes) are expected to have greater potential to influence the promoters than those 3’ to genes. If TE-mediated epigenetic effects could indeed lower gene expression, the associations should be more likely observed with TEs 5’ to genes. Analyzing TEs 5’ and 3’ to genes together could thus potentially obscure the signal. To test this possibility, we first investigated whether TEs 5’ and 3’ to genes have differential impacts on the epigenetic states of genes. We still observed significant positive associations between the magnitude and extent of TE-mediated H3K9me2 enrichment and genic H3K9me2 enrichment for both TEs 5’ and 3’ to genes (Spearman rank correlation coefficient ρ=0.13–0.38 [TEs 5’ to genes] and 0.16–0.43 [TEs 3’ to genes], p<0.05 for all tests except for the magnitude of the effect for TEs 5’ to genes, D. yakuba strain 2; Figure 3B). Also, there is no consistent trend on whether the correlation is stronger for TEs 5’ or 3’ to genes across genomes (Figure 3B). The dependency of this association on TE-gene distance also generally holds for both TEs 5’ and 3’ to genes (p-value for interaction term <0.05 for more than half of the genomes, Figure 3—figure supplement 5). On the other hand, we still observed nearly absent associations between TE-mediated epigenetic effects and gene expression rank even when analyzing TEs 5’ and 3’ to genes separately, and the few significant associations are all for TEs 5’ to genes (Spearman rank correlation test, p>0.05 for all tests except for TEs 5’ to genes, D. yakuba strain 1 [both magnitude and extent] and D. yakuba strain 2 [extent] Figure 3B). Again, the associations between TE-mediated epigenetic effects and gene expression do not depend on TE-gene distance either for TEs 5’ and 3’ to genes for most of the cases (regression analysis interaction term, p>0.05 for all genomes and for both TEs 5’ and 3’ to genes, except for the magnitude for TEs 5’ to genes [D. simulans strain 1] or for the extent for TE 3’ to genes [D. mauritiana], Figure 3—figure supplement 6).

The above analyses may have limited power because variation in gene expression levels within a genome could be due to intrinsic gene properties, instead of TE-mediated effects. To directly test the effects of TE-mediated H3K9me2 enrichment on nearby gene expression, we compared the expression of homologous alleles with and without adjacent TEs in D. simulans and D. yakuba, two species that we have epigenomic and transcriptomic data for two strains. We estimated z-scores, which compare H3K9me2 enrichment and expression rank between homologous alleles (see Materials and methods). A positive z-score indicates that the allele adjacent to TE insertions has higher H3K9me2 enrichment or expression rank (i.e., lower expression) than the homologous allele without a nearby TE. We again observed TEs with stronger epigenetic effects associated with higher z-scores of genic H3K9me2 enrichment, which confirms their impacts on genic epigenetic states (Mann-Whitney U test, p<0.05 for all comparisons; Figure 3C and Figure 3—figure supplement 7 for categorizing genes according to the magnitude and extent of TE-induced H3K9me2 enrichment, respectively). By analyzing TEs 5’ and 3’ to genes, we observed a consistent trend that z-scores for genes near TEs with epigenetic effects are larger, even though the comparisons are only significant for a subset of comparisons, probably due to the reduced sample size (Figure 3—figure supplement 8). On the contrary, there is no association between TE-mediated epigenetic effects and z-scores of gene expression rank (Mann-Whitney U test, p>0.05 for all comparisons; Figure 3C and Figure 3—figure supplement 7 for categorizing genes according to the magnitude and extent of TE-induced H3K9me2 enrichment, respectively; Figure 3—figure supplement 8 for looking at TEs 5’ and 3’ to genes separately). Directly comparing the expression of homologous alleles with and without TEs reached similar conclusions (Figure 3—figure supplement 9). Overall, our results suggest that, across genomes, TE-mediated epigenetic effects lead to robust enrichment of heterochromatic marks at neighboring genes, which, contrary to expectation, does not have a predominantly negative impact on the expression of neighboring genes. It is worth noting that our analyses focus on identifying the genome-wide average pattern; it is still plausible that the local enrichment of H3K9me2 induced by individual TEs occasionally lowers nearby gene expression (e.g., genes with positive z-scores in Figure 3C or genes with TE-mediated regulation; Ninova et al., 2020). In fact, many early examples of TE-mediated epigenetic effects were discovered by the phenotypic consequences of reduced neighboring gene expression (reviewed in Choi and Lee, 2020).

The extent of local enrichment of repressive epigenetic marks for a handful of TE insertions was previously suggested to be influenced by the expression of neighboring genes in a mouse cell line (Rebollo et al., 2012). This observation suggests the possibility that the extent of TEs’ epigenetic effects is, in return, restrained by the expression of neighboring genes. If true, this phenomenon may explain our observed lack of negative associations between TE-mediated epigenetic effects and neighboring gene expression. To investigate this possibility on a genome-wide scale, we compared the H3K9me2 enrichment for a TE on the side that faces a gene with appreciable expression in 16–18 hr embryo (at least 10 RPKM; genic side) and the other side that does not face the gene (intergenic side; Figure 3D). We predict that the transcriptional effects of genes on TE-mediated enrichment of H3K9me2 should be the most prominent on the ‘genic side’ of a TE. Also, the differences between the two sides of a TE should be larger for TEs closer to genes.

To test these predictions, we estimated the normalized difference between TE-mediated H3K9me2 enrichment on the intergenic and genic sides by calculating a z-score. A positive z-score would mean that the magnitude or extent of TE-mediated H3K9me2 enrichment is more restricted on the genic side than on the intergenic side. Consistent with these predictions, there is a general trend that the z-score for the extent of spread is positive for TEs close to highly expressed genes (RPKM >10), although the comparisons are statistically significant only for a subset of the genomes (Mann-Whitney U test, p<0.001 for D. simulans strain 1 and D. teissieri, Figure 3D). On the other hand, z-scores for the extent of spread are not significantly different from 0 for TEs far from highly expressed genes (Mann-Whitney U test, p>0.05 for all genomes, Figure 3D) and are significantly larger for TEs close to highly expressed genes than for those far from highly expressed genes (Mann-Whitney U test, p<0.05 for D. simulans strain 1, D. teissieri, and D. santomea; Figure 3D). We also found significant negative associations between the z-score of a TE and its distance to the nearest gene (Spearman rank correlation coefficient ρ=–0.17 [D. mauritiana], –0.13 [D. simulans strain 1], –0.11 [D. teissieri], and –0.10 [D. santomea], p<0.01 for all tests), further supporting that the restricted extent of TE-mediated H3K9me2 enrichment on the genic side depends on TE-gene distance. Curiously, comparisons based on the magnitude of TE-induced H3K9me2 enrichment did not find differences between the genic side and the intergenic side, nor between TEs that are close to or far from highly expressed genes (Mann-Whitney U test and Spearman rank correlation test, p>0.05 for all comparisons, Figure 3—figure supplement 10; also see Discussion).

If the restricted extent of TE-mediated H3K9me2 enrichment on the genic side indeed results from the transcription activities of neighboring genes, this effect should mainly be observed with TEs near highly, but not lowly, expressed genes and would be stronger for TEs 5’ to genes (i.e., near promoters) than those 3’ to genes. Consistent with these predictions, z-score for TEs close to lowly expressed genes (RPKM <10) are not significantly different from 0, nor do they differ between TEs close to or far from genes (Mann-Whitney U test, p>0.05 for both comparisons in all genomes, Figure 3—figure supplement 11). Interestingly, TEs 5’ to highly expressed genes are more likely to have z-scores significantly larger than 0 (for four out of eight genomes) than those of TEs 3’ to genes (for one out of eight genomes), suggesting a more restricted extent of H3K9me2 enrichment on the genic side for TEs near promoters (Mann-Whitney U test, p<0.05 for D. simulans strain 1, D. teissieri, D. santomea, and D. yakuba strain 2 [TEs 5’ to genes] and for D. teissieri [TEs 3’ to genes], Figure 3—figure supplement 12). Similarly, the negative associations between the z-score of a TE and its distance to the nearest gene are significant for TEs 5’ to genes in four genomes (Spearman rank correlation coefficient ρ=–0.24 [D. mauritiana], –0.18 [D. simulans strain 1], –0.11 [D. teissieri], and –0.17 [D. santomea], p<0.05 for all tests) but not for those 3’ to genes (Spearman rank correlation tests, p>0.05 for all tests). Such observation is consistent with the idea that the distance-dependent effect of gene transcription on the extent of TE-mediated epigenetic effects only holds for TEs near promoters. Overall, our findings reveal a complex relationship between TE-mediated epigenetic effects and neighboring gene expression, and strongly suggest that the extent of TE-mediated H3K9me2 enrichment is influenced by the expression of adjacent genes, especially on the side of a TE that faces the promoter of a nearby, highly expressed gene.

TEs with epigenetic effects are selected against across species

TEs exerting epigenetic effects were previously suggested to experience stronger purifying selection than other TEs (reviewed in Choi and Lee, 2020). If TE-mediated epigenetic effects indeed impair host fitness and are thus selected against, such effects could play important roles in shaping the evolution of both TEs and their host genomes. Yet, previous analyses testing the presence of selection against TEs with epigenetic effects were restricted to few model species. More importantly, the confounding effects of other deleterious mechanisms of TEs (e.g., ectopic recombination; see below) on TE population frequencies could not be ruled out in those studies (Hollister and Gaut, 2009; Lee, 2015; Lee and Karpen, 2017).

To investigate the fitness impacts of TE-mediated H3K9me2 enrichment, we estimated the frequencies of TEs that are included in our analysis in four species, using previously published populations of genomes (D. melanogaster [Lack et al., 2015], D. simulans and D. yakuba [Rogers et al., 2014], and D. mauritiana [Garrigan et al., 2014]). Population frequencies of TEs are strongly influenced by the strength of natural selection against their deleterious effects (reviewed in Charlesworth and Langley, 1989; Lee and Langley, 2010; Barrón et al., 2014), with low-frequency TEs generally expected to be more deleterious than high-frequency TEs. We categorized TEs in our focused genome into high-frequency and low-frequency according to whether they are (high-frequency) or are not (i.e., singletons; low-frequency) identified in the population samples. We found that the magnitude of TE-induced H3K9me2 enrichment is significantly higher for low-frequency TEs than for high-frequency TEs in all four species (Mann-Whitney U test, p<0.05 for all but D. simulans strain 1, Figure 4), which extends previous investigations in D. melanogaster. Intriguingly, the difference is much weaker when comparing the extent of TE-mediated H3K9me2 enrichment, and the comparison is significant for only one genome (Mann-Whitney U test, p<0.05 for strain 2 of D. simulans, but >0.05 for all other genomes; Figure 4, see Discussion).

Associations between the population frequencies and epigenetic effects of transposable elements (TEs).

The magnitude (top) and extent (bottom) of TE-mediated enrichment of H3K9me2 for low-frequency TEs (usually considered as strongly selected, green) and high-frequency TEs (gray) in Drosophila melanogaster, Drosophila mauritiana, Drosophila simulans, and Drosophila yakuba. Low-frequency TEs are those that are only found in the focused genome, while high-frequency TEs are identified in both the focused genome and the population of genomes. Mann-Whitney U test, **p<0.01, *p<0.05.

A potential confounding factor for our observed associations between TEs’ epigenetic effects and population frequencies is TE length. TE length was previously observed to negatively correlate with population frequencies of TEs (Petrov et al., 2003; Petrov et al., 2011), either because of the larger potential of long TEs in disrupting functional elements or their higher propensity to be involved in deleterious ectopic recombination (Petrov et al., 2003). Because the strength of TE-mediated epigenetic effects also positively correlated with TE length (Figure 1—figure supplement 3), our observed negative associations between TE frequencies and epigenetic effects could instead result from other harmful effects of TEs. To investigate this possibility, we performed logistic regression analysis to test the effects of TE-mediated H3K9me2 enrichment on population frequencies while accounting for the influence of TE length (population frequency ~ length + TE epigenetic effects). Because of the co-linearity between predictor variables in the regression model (TE length and epigenetic effects, Figure 1—figure supplement 3), this analysis is expected to have restricted statistical power. Nevertheless, regression coefficients for the magnitude of H3K9me2 enrichment are negative for all but one genome (Table 1), and the coefficient is significantly negative for D. melanogaster.

Table 1
Logistic regression coefficients for the effects of transposable element (TE) length and TE-mediated H3K9me2 enrichment (magnitude and extent) on the population frequencies of TEs.
Magnitude of H3K9me2 enrichmentExtent of H3K9me2 enrichment
TE lengthMagnitudeTE lengthExtent
D. melanogaster–4.95E-044.59E-01*–5.22E-048.63E-05
D. mauritiana–1.83E-032.74E-02–1.87E-031.55E-05
D. simulans (strain 1)–2.56E-041.07E-01–2.52E-049.00E-06
D. simulans (strain 2)–3.13E-032.50E-01–3.21E-032.41E-04
D. yakuba (strain 1)–3.63E-041.68E-01–3.73E-047.77E-05
D. yakuba (strain 2)–4.66E-041.27E-01–4.72E-042.87E-04
  1. Negative regression coefficients for TE-mediated epigenetic effects on TE population frequencies are in bold. *p<0.05.

Local meiotic recombination rate is another factor that could potentially confound our observed negative associations between TE population frequencies and epigenetic effects. If TEs with weaker epigenetic effects tend to locate in genomic regions with low recombination rate, lower probability of ectopic recombination (Langley et al., 1988), or reduced efficacy of selection (Hill and Robertson, 1966; Felsenstein, 1974; reviewed in Charlesworth and Langley, 1989; Lee and Langley, 2010; Barrón et al., 2014; Kent et al., 2017) may instead drive their higher population frequencies. However, we found no associations between TE-mediated epigenetic effects and local recombination rate in D. melanogaster, the species with a high-resolution recombination map (Comeron et al., 2012) (Spearman rank correlation tests, p=0.59 [magnitude] and 0.98 [extent]). Such observation is consistent with a previous analysis using different D. melanogaster strains (Lee and Karpen, 2017) and suggests that the difference in meiotic recombination rate unlikely confounded our analysis of TE population frequencies.

The average deleterious effects of TE insertions were found to differ between TE families (reviewed in Charlesworth and Langley, 1989; Barrón et al., 2014), which could also confound our analysis, given that the strength of TE-mediated epigenetic effects varies along the same axis (Figure 2). Accordingly, we studied the associations between TE-induced H3K9me2 enrichment and TE population frequencies among copies of the same TE family within species, aiming to exclude the potential confounding effects of family identity on TE population frequencies. We performed logistic regression analysis for individual TE families that have at least 10 insertions while accounting for the effects of TE length (population frequency ~ TE epigenetic effects + length). With collinearity among predictor variables (see above) and the small number of TEs included for each TE family (fewer than 30 copies for 19 out of 23 TE families, with a median of 19 of TEs included in the regression analysis), these analyses are again underpowered. For most of the TE families tested, we observed negative regression coefficients for the magnitude of TE-mediated H3K9me2 enrichment (Table 2), which is more than expected by chance (18 out of 23 TE families tested, binomial test, p=0.0106). We observed similar results with the extent of TE-mediated H3K9me2 enrichment (17 out of 23 TE families tested, binomial test, p=0.0346, Table 2). Despite the lack of statistical power, H and Cr1a families showed significant negative regression coefficients for TE-mediated epigenetic effects in one of the D. yakuba strains. Overall, our results revealed that, after controlling for the effects of TE length and family identity, we still found negative associations between the strength of TE-mediated epigenetic effects and TE population frequencies. By excluding the potential impacts of confounding factors on TE population frequencies, our observation strongly supports that TE-mediated enrichment of repressive marks is disfavored by natural selection in multiple species.

Table 2
Logistic regression coefficients for the effects of transposable element (TE) length and TE-mediated H3K9me2 enrichment (magnitude and extent) on the population frequencies of TEs from different families.
Magnitude of H3K9me2 enrichmentExtent of H3K9me2 enrichment
TE familyTE lengthMagnitudeTE lengthExtent
D. melanogaster
BS–8.10E-032.65E+00–5.18E-035.16E+02
297–4.86E-041.43E+00–5.12E-041.10E-04
jockey2.17E-041.30E+002.50E-041.42E-03
pogo9.46E-041.04E+009.18E-049.47E-05
Doc2.65E-042.20E-012.96E-042.55E-04
hopper–1.53E-041.36E+00–1.63E-042.14E-04
D. mauritiana
HB–2.94E-032.10E-01–2.96E-033.23E-05
Bari–4.78E-036.58E-01–4.96E-033.65E-04
hopper6.66E-031.96E+008.70E-033.72E-04
D. simulans (strain 1)
H1.84E-029.70E-021.76E-023.08E-04
transib2–3.69E-031.26E+00–3.71E-033.08E-04
Tc1–1.19E-013.56E-01–1.16E-012.33E-04
1,360–3.84E-021.37E-01–4.02E-028.88E-05
diver22.16E-041.39E-01–2.54E-041.39E-03
Helena–3.99E-032.85E-01–3.62E-034.20E-04
HB–2.45E+004.26E+00–2.48E+001.61E-03
roo–4.01E-042.71E-01–5.67E-044.10E-04
D. simulans (strain 2)
H–4.97E-036.63E-01–4.94E-033.02E-05
Tc1–2.28E-038.68E-01–2.46E-037.59E-04
D. yakuba (strain 1)
H–1.80E-033.79E-01–1.50E-032.75E-04
Cr1a9.53E-056.66E-01–1.05E-043.35E-04
D. yakuba (strain 2)
H–2.29E-023.89E-01–5.48E-028.98E-03*
Cr1a1.32E-031.32E+00*4.42E-041.27E-04
  1. Negative regression coefficients for TE-mediated epigenetic effects on TE population frequencies are in bold. *p<0.05.

Epigenetic effects of TEs negatively associate with genomic TE abundance across species

The magnitude and extent of TE-mediated H3K9me2 enrichment vary not only within genomes and between individuals, but also between species (Figure 1). According to a previously proposed hypothesis (Lee and Karpen, 2017), this between-species difference in the strength of TE-mediated H3K9me2 enrichment could determine the average strength of selection removing TEs, leading to species-specific TE abundance in the euchromatic genome. To test this prediction, we investigated the associations between euchromatic TE numbers and the average magnitude and extent of TE-induced local enrichment of H3K9me2. Because the assignment of TEs into families is biased against TEs in yakuba complex species (Figure 2—figure supplement 3), all TEs, irrespective of whether we could assign their family identity, were included in this between-species analysis (see Materials and methods).

We found that TEs in species of the yakuba complex show a much larger magnitude of H3K9me2 enrichment than those in species of the melanogaster complex (Figure 5A, Mann-Whitney U test, p=0.029). After controlling for the strong impacts of species complex, we found significant negative associations between the number of euchromatic TEs and the average magnitude of TE-mediated H3K9me2 enrichment across genomes (Figure 5A, TE abundance ~ species complex + epigenetic effects; regression coefficient for the magnitude of H3K9me2 enrichment: –1603.5; ANOVA F-value: 9.05, df = 1, p=0.030). Exclusion of D. simulans strain 2, a sample that shows lower H3K9me2 enrichment consistency between replicates (see Materials and methods) gave consistent results (regression coefficient: –1602; ANOVA F-value: 5.26, df = 1, p=0.080). On the other hand, the extent of TE-mediated H3K9me2 spreading does not differ between species complexes (Mann-Whitney U test, p=0.89) nor associate with the number of euchromatic TEs (Figure 5A, regression coefficient for the extent of H3K9me2 enrichment [kb]: 925; ANOVA F-value: 0.03, df = 1, p=0.87; excluding D. simulans strain 2 – regression coefficient: 1028, F-value: 0.13, df = 1, p=0.78; see Discussion). This finding echoes our observations of the nearly absent associations between the extent of TE-mediated H3K9me2 enrichment and TE population frequencies (Figure 4). It is worth noting that the associations between the magnitude, but not the extent, of TE-mediated H3K9me2 enrichment and genomic TE abundance was also documented in a much smaller study previously (Lee and Karpen, 2017). To account for the phylogenetic non-independence among species, we used phylogenetic generalized least squares (PGLS, Grafen, 1989; Martins and Hansen, 1997) to repeat the regression analysis and found consistent results (regression coefficient for the magnitude of H3K9me2 enrichment: –1948; ANOVA F: 7.49, df = 1, p=0.034; regression coefficient for the extent of H3K9me2 enrichment [kb]: 599; ANOVA F: 0.14, df = 1, p=0.72). Overall, the negative associations between the magnitude of TE-mediated H3K9me2 enrichment and abundance of euchromatic TEs within species complex support our prediction that varying strength of TE-mediated epigenetic effects could drive between-species differences in genomic TE abundance on a short evolutionary time scale.

Transposable element (TE)-mediated epigenetic effects associate with genomic TE abundance and repressive chromatin landscape.

(A) The associations between the mean magnitude (left)/extent (right) of TE-mediated enrichment of H3K9me2 and the number of estimated euchromatic TEs across species. (B) An example gene (bin3) whose expression rank negatively correlates with the mean magnitude of TE-mediated H3K9me2 enrichment across species (left). The distributions of Spearman rank correlation coefficient (ρ) between the magnitude of TE-mediated H3K9me2 enrichment and the expression ranks significantly differ between Su(var) genes (green) and other genes in the genome (gray; right). (C) An example gene (Ago2) whose expression rank negatively correlates with the magnitude of TE-mediated H3K9me2 enrichment between species within species complex, but differs significantly between species complex (left). The distributions of regression coefficients for the effect of a gene’s expression on the magnitude of TE-mediated H3K9me2 enrichment significantly differs between Su(var) genes (green) and other genes in the genome (gray, right). (D) The reduced dosage of Su(var) genes influences the epigenetic silencing effect of 1360. For each candidate Su(var), we performed three replicates (three independent crosses), and one dot represents one cross. (E) The mean magnitude of TE-mediated enrichment of H3K9me2 associates with the abundance of H3K9me2-enriched Kmers. **p<0.01 and *p<0.05 for Mann-Whitney U test or Kolmogorov-Smirnov test (see text).

Species-specific repressive chromatin landscape associates with between-species differences in epigenetic effects of TEs

If varying strength of TE-mediated epigenetic effects indeed contributes to between-species differences in euchromatic TE copy number, as suggested by our observations (Figure 5A), species-specific differences in genetic factors that modulate the epigenetic effects of TEs could be the ultimate drivers for varying genomic TE abundance. Consistent with this possibility, we observed that the strength of TE-mediated epigenetic effects strongly depends on the host genetic background (Figure 2). The similarities between TE-mediated enrichment of repressive marks and the well-studied ‘position effect variegation’' (PEV, Gowen and Gay, 1934) suggest a feasible path to narrow down the causes for variation in the epigenetic effects of TEs. Repressive epigenetic marks at constitutive heterochromatin are long known to spread to juxtaposed euchromatic genes, a phenomenon known as PEV (reviewed in Elgin and Reuter, 2013). Several genetic factors have been identified to modulate the strength and extent of PEV (reviewed in Girton and Johansen, 2008; Elgin and Reuter, 2013). These include Su(var) genes, whose protein products act as structural or enzymatic components of the heterochromatin (Girton and Johansen, 2008; Elgin and Reuter, 2013; Swenson et al., 2016), and heterochromatic repeats, which are targets of heterochromatin formation (Girton and Johansen, 2008; Elgin and Reuter, 2013). The ratio between the dosage of these two players (targets and regulators) was postulated to influence the nucleation and formation of constitutive heterochromatin and, accordingly, the spread of repressive marks from constitutive heterochromatin to the euchromatic genome (Locke et al., 1988).

In order to test the hypothesized role of Su(var)s and heterochromatic repeats in shaping between-species differences in TE-mediated epigenetic effects, multi-species comparisons that span a reasonable phylogenetic resolution would be needed. Also, it would be important to demonstrate that Su(var)s that regulate the spreading of heterochromatic marks from constitutive heterochromatin play a similar role in the epigenetic effects of TEs in the euchromatic genome. This is because, though both enriched with H3K9me2/3, heterochromatin is comprised of large blocks of repetitive sequences enriched for repressive epigenetic marks (Riddle et al., 2011), whereas euchromatic TEs, though epigenetically silenced, are relatively short (Lee, 2015; Lee and Karpen, 2017) and usually surrounded by sequences enriched for active epigenetic marks (Kharchenko et al., 2011).

To test our hypothesis that species-specific differences in genetic factors contribute to varying epigenetic effects of TEs, we investigated whether there are positive associations between TE-mediated epigenetic effects and the dosage of Su(var) genes, which promotes the spreading of repressive marks from constitutive heterochromatin (Girton and Johansen, 2008; Elgin and Reuter, 2013; Swenson et al., 2016). We used the expression of Su(var) genes in 16–18 hr embryos, a developmental stage matching our epigenomic data, as a proxy for the dosage of Su(var) protein products. Because we only found that the magnitude, but not the extent, of TE-mediated H3K9me2 enrichment negatively associates with TE copy number (Figure 5A), our following analysis focused on the magnitude of TEs’ epigenetic effects, with the aim of identifying the ultimate cause for varying genomic TE abundance. We estimated the Spearman rank correlation coefficients (ρ) between the expression rank of a gene (rank from the highest expressed genes) and the average magnitude of TE-mediated H3K9me2 enrichment across genomes (see Figure 5B for an example, also see Materials and methods). A negative correlation indicates that a Su(var) gene has higher expression (and thus lower expression rank) in a genome with stronger epigenetic effects of TEs. We found that Su(var) genes, as a group, have significantly lower ρ than other genes in the genome (Mann-Whitney U test, p=0.0073) and have a shifted distribution of ρ toward smaller values (Kolmogorov-Smirnov test, p=0.033, Figure 5B). These observations suggest that the expression levels of Su(var)s correlate more positively with the magnitude of TE-mediated H3K9me2 enrichment than other genes in the genome. Among analyzed Su(var)s, the ρ for bin3, promoter of small-RNA-mediated silencing (Singh et al., 2011), is among the top 5% of all genes. Lhr and HP4, both of whose protein products are structural components of heterochromatin (Greil et al., 2007), and Hsc70-4, whose protein product is an interactor of core heterochromatin protein HP1a (Swenson et al., 2016), are among top 10% genome-wide (Supplementary file 1).

Intriguingly, the expression rank of Ago2, a key gene in initiating epigenetic silencing and a Su(var) (Deshpande et al., 2005), is among the top 10% genome-wide that positively correlates with TE-mediated enrichment of H3K9me2 (Supplementary file 1), an association that is opposite to prediction. Upon further examination, we found that the expression of Ago2 significantly differs between species in the two species complexes, which drives the positive correlation (Figure 5C). Accordingly, we performed regression analyses that associate the magnitude of TE-mediated H3K9me2 enrichment and gene expression rank while accounting for the effects of species complex (magnitude of TE-mediated epigenetic effects ~ gene expression + species complex). Similar to the analysis based on ρ (see above, Figure 5B), the regression coefficients are significantly smaller for Su(var)s than those for other genes (Mann-Whitney U test, p=0.063; Kolmogorov-Smirnov test, p=0.032; Figure 5C), further corroborating our findings. With the regression analysis, we found that the association between Ago2 expression rank and the magnitude of TE-mediated epigenetic effects is negative and is among 10% genome-wide. We also found other Su(var) genes whose regression coefficients are among the top 10% of all genes, which includes RpLP0, whose protein product is a core interactor of HP1a (Frolov and Birchler, 1998), and Su(var)3–3, which codes for an eraser of active histone modification (Rudolph et al., 2007, Supplementary file 1).

To confirm that these identified candidate Su(var)s also modulate TE-mediated local enrichment of H3K9me2 in the euchromatic genomes, we leveraged a previously published reporter system that allows the quantification of TE-mediated epigenetic effects (Sentmanat and Elgin, 2012). In this system, a DNA-based TE, 1360, results in the enrichment of H3K9me2/3 at the immediate adjacent mini-white reporter gene in the euchromatic genome. The same study also found associations between the presence of 1360, the enrichment level of H3K9me2 at the mini-white reporter gene, and the reduced amounts of red-eye pigmentation (Sentmanat and Elgin, 2012). Accordingly, the eye pigmentation level in this system serves as a convenient readout for the magnitude of TE-mediated H3K9me2 enrichment. By using existing loss-of-function Su(var) mutants, we found that reduced dosage (hemizygous) of candidate Su(var)s leads to significantly elevated levels of eye pigmentation, or reduced TE-mediated silencing effects (Figure 5D), supporting the roles of candidate Su(var)s in modulating the epigenetic effects of euchromatic TEs. With the assumption that the functional roles of Su(var) genes are conserved across Drosophila species studied, our findings suggest that the observed species-specific expression of Su(var)s could drive between-species differences in TE-mediated H3K9me2 enrichment in the euchromatic genome.

We also investigated whether the epigenetic effects of TEs negatively associated with the abundance of heterochromatic repeats, which was found to weaken the spreading of H3K9me2/3 from constitutive heterochromatin (reviewed in Girton and Johansen, 2008; Elgin and Reuter, 2013). We identified Kmers that are enriched with H3K9me2 in our ChIP-seq data and quantified their abundance using Illumina sequencing with PCR-free library preparation, in an effort to avoid biases in quantifying simple repeats (Wei et al., 2018) (see Materials and methods). Consistent with the prediction that repressive chromatin landscape weakens with increased abundance of heterochromatic repeats, we found that the magnitude of TE-mediated H3K9me2 enrichment is negatively associated with the abundance of H3K9me2-enriched repeats between species, though the comparison is not significant (Figure 5E, Spearman rank correlation coefficient ρ=–0.595, p=0.13). Overall, our observations support the prediction that species-specific repressive chromatin landscape, which is shaped by the expression level of Su(var) genes and abundance of heterochromatic repeats, associates with differences in TE-mediated epigenetic effects between species.

Discussion

The replicative nature of TEs has made them successful at occupying nearly all eukaryotic genomes. Yet, their ‘success’, or genomic abundance, drastically varies across the phylogenetic tree (Huang et al., 2012; Elliott and Gregory, 2015; Wells and Feschotte, 2020) and between closely related species (Hu et al., 2011; Rius et al., 2016; Legrand et al., 2019), raising important questions regarding the evolutionary causes and functional consequences of varying genomic TE abundance. It is worth noting that most studies about the evolution of TE abundance, including ours, mainly concern insertions in the euchromatic genome. Although TEs are highly abundant in the heterochromatic genome, these TEs are typically fragmented (Hoskins et al., 2007; Hoskins et al., 2015) and would no longer be part of the ‘life cycle’ of TEs. The role of euchromatic TEs in determining the abundance and evolutionary dynamics of TEs in the heterochromatic genome is still a largely unaddressed question, mainly due to the challenges associated with assembling TEs in repeat-rich heterochromatic sequences (Salzberg and Yorke, 2005; Treangen and Salzberg, 2011, but see Khost et al., 2017; Chang and Larracuente, 2019; Chang et al., 2019).

Several hypotheses have been proposed to explain the large variation in euchromatic TE abundance. Phylogenetic signals may explain some of the differences between distantly related taxa (Wells and Feschotte, 2020) and, sometimes, between species within a taxa (Szitenberg et al., 2016). Variation in TE activities was also postulated to contribute to the wide variability of TE abundance (Chen et al., 2019; Wong et al., 2019, but see Ho et al., 2021). Still another plausible cause is systematic differences in the strength of natural selection removing TEs between species. Investigations of this hypothesis have been largely focused on population genetic parameters that influence the efficacy of natural selection, such as mating systems (Wright and Schoen, 1999; Dolgin et al., 2008; Boutin et al., 2012; Arunkumar et al., 2015; Agren et al., 2014) and effective population size (Lynch and Conery, 2003; Mérel et al., 2021). Here, our multi-species study tested step-by-step yet another evolutionary mechanism by which the strength of selection removing TEs could differ—through differences in the epigenetic effects of TEs that are driven by species-specific host chromatin landscape.

In this study, we investigated the prevalence, variability, and evolutionary importance of ‘the epigenetic effects of TEs’—TE-mediated local enrichment of repressive epigenetic marks—in the euchromatic genomes of six Drosophila species. These species come from two important and well-studied complexes within the melanogaster species subgroup (melanogaster and yakuba complexes), providing good phylogenetic resolution to decipher the evolutionary role of epigenetic mechanisms in shaping species-specific TE abundance. We observed that TE-mediated local enrichment of repressive marks is prevalent in all species studied and widely varies both within and between genomes. Interestingly, in addition to the intrinsic biological properties of TEs, our analyses revealed that host genetic background plays a critical role in determining the strength of TE-mediated epigenetic effects. These TE-mediated enrichments of repressive marks alter the epigenetic states of neighboring euchromatic sequences, including actively transcribing genes. Importantly, TEs exerting such epigenetic effects are selected against across multiple species, and the strength of this TE-mediated epigenetic effect negatively correlates with genomic TE abundance within species complex. Our findings extend previous studies based on few species and provide one of the first support for the importance of the inadvertent harmful effects of TE epigenetic silencing in shaping divergent genomic TE landscapes.

Curiously, we only found a negative association between genomic TE abundance and the magnitude of TE-mediated H3K9me2 enrichment across species, but not the extent of the spreading, even though these two indexes of a TE significantly correlate within genomes (Figure 1—figure supplement 2). Similarly, the evidence for selection against TE-mediated epigenetic effects is stronger for the magnitude of TE-induced H3K9me2 enrichment than for the extent (Figure 4). Previous empirical (Belyaeva and Zhimulev, 1991; Talbert and Henikoff, 2000) and theoretical (Erdel and Greene, 2016) studies have suggested that the spreading of repressive heterochromatic marks is a nonlinear process. Accordingly, the extent of the spreading of repressive marks from TEs could be sensitive to the genomic context and subject to substantial stochasticity. This echoes our observations that the transcription of genes influences the extent of TE-mediated H3K9me2 spreading, but not the magnitude of H3K9me2 enrichment (Figure 3D and Figure 3—figure supplement 10; also see below). The magnitude of TE-mediated H3K9me2 enrichment, which is estimated in windows right next to TE boundary, may thus be a more direct measurement of the strength of TE-mediated epigenetic effects and more indicative of the associated harmful impacts.

One of our most surprising findings is perhaps the limited evidence in supporting that TE-mediated epigenetic effects reduce neighboring gene expression on a genome-wide scale. While TE-mediated epigenetic effects increase the enrichment of repressive epigenetic marks at adjacent genes, we found limited associations between such effects and gene expression within genomes (Figure 3A and B). Comparing the expression of homologous alleles with and without TEs showing epigenetic effects also reached the same conclusion (Figure 3C). These observations echo previously reported lack of genome-wide associations between TE-induced enrichment of repressive epigenetic marks and the expression of neighboring genes in several model organisms (Quadrana et al., 2016; Stuart et al., 2016; Lee and Karpen, 2017; Choi and Purugganan, 2018) and are consistent with the wider observation that TEs upstream or downstream to genes are not predominantly associated with reduced gene expression (Goubert et al., 2020; Ullastres et al., 2021; Rech et al., 2022, reviewed in Kelleher et al., 2020; Choi and Lee, 2020). The limited impacts of TE-mediated H3K9me2 enrichment on gene expression could have resulted from the complex relationship between repressive epigenetic modification and gene expression (de Wit et al., 2007; Yasuhara and Wakimoto, 2008; Riddle et al., 2011; Meng et al., 2016; Caizzi et al., 2016), the varying sensitivity of genes to the enrichment of repressive marks (Rudolph et al., 2007; Vogel et al., 2009; Riddle et al., 2011), and the presence of other types of variants that also modulate gene expression (Stranger et al., 2007). Interestingly, our findings suggested another possibility—the transcription of genes, in return, influences TE-mediated epigenetic effects across genomes. Specifically, we found a more restricted spreading of repressive marks from TEs on the side facing a gene than on the intergenic side, and this difference is mainly observed for TEs near the promoters of highly expressed genes (Figure 3—figure supplement 12). It is worth noting that this difference in the extent of H3K9me2 spreading between the two sides of a TE is unlikely driven by the presence of insulator sequences (Gaszner and Felsenfeld, 2006) because of the limited associations between the presence of insulator sequences and the reduction in TE-mediated H3K9me2 spreading on the genic side (Figure 3—figure supplement 13). On the other hand, active histone modifications enriched at transcriptionally active genes could antagonize the assembly of heterochromatin (Allshire and Madhani, 2018), potentially restraining the spreading of repressive marks from silenced TEs. Consistently, in mice, active histone modification was reported to spread from an actively transcribing candidate gene into an adjacent TE (Rebollo et al., 2012). Similarly, active gene transcription at Drosophila miranda neo-Y chromosome was found to impede the formation of heterochromatin (Wei et al., 2020). Future analysis of the distributions of active histone modifications may help reveal the mechanistic cause for our observed genome-wide dependencies of the extent of TE-mediated epigenetic effects on the transcriptional activities of adjacent euchromatic genes.

If not due to reducing the expression of adjacent genes, what functional consequences of TE-mediated epigenetic effects could have impaired host fitness and led to the observed selection against these TEs? Euchromatic TEs enriched with repressive epigenetic marks were reported to spatially interact with constitutive heterochromatin through phase separation mechanisms, a process observed to alter 3D structures of genomes and inferred to lower host fitness (Lee et al., 2020). In addition, the epigenetic effects of TEs could shift the usual DNA repair process in the gene-rich euchromatic genome, perturbing the maintenance of local genome integrity. This is because double-stranded breaks happening in constitutive heterochromatin are repaired through a distinct cellular process from those in the euchromatic genome, owing to the enrichment of repressive epigenetic modifications (Chiolo et al., 2011; Janssen et al., 2016; Janssen et al., 2019). Still, the variance of gene expression was shown to be shaped by natural selection (Metzger et al., 2015; Duveau et al., 2018). Given the variegating properties of the spreading of repressive marks (reviewed in Elgin and Reuter, 2013), TE-mediated epigenetic effects could have shifted the variance, instead of the mean, of neighboring gene expression, impacting host fitness.

Building upon our findings of the prevalence and variability of TE-mediated epigenetic effects and selection against such effects, our multi-species analysis provides strong support for the previously proposed role of TE-mediated epigenetic effects in determining genomic TE abundance (Figure 6). This interpretation is based on the assumption that TE copy number is at equilibrium and transposition rates of TEs are similar across species, leaving the strength of selection removing TEs being the major determinant of genomic TE abundance (Charlesworth and Charlesworth, 1983). By combining transcriptomic analysis and Drosophila genetics experiment, we further revealed that between-species differences in TE-mediated epigenetic effects could be driven by species-specific expression levels of host genetic factors that modulate heterochromatin. These findings connect the evolution of host genome (genomic TE abundance) with chromatin landscape through the inadvertent harmful effects of the epigenetic silencing of TEs (Figure 6). Curiously, the negative association between TE-mediated H3K9me2 enrichment and genomic TE abundance was only observed within species complex, pointing that the importance TE-mediated epigenetic effects in determining genomic TE abundance may be on a short evolutionary time scale. Specifically, even though we found stronger epigenetic effects of TEs in species of the yakuba complex than those in the melanogaster complex, TEs are more abundant in the former. This accumulation of TEs in the genomes of yakuba complex is in line with previously observed greater net gains of intron sequences (Presgraves, 2006) and duplicated genes (Hahn et al., 2007) in D. yakuba than those in the lineage leading to melanogaster complex, which could suggest a higher tolerance of D. yakuba genome with gained sequences. Such lineage-specific mechanism that shapes the gain and loss of sequences over the entire genome may also play an important role in determining genomic TE abundance and potentially override the influence of selection against varying TE-mediated epigenetic effects on a longer evolutionary time scale.

Proposed role of the chromatin landscape in determining genomic transposable element (TE) abundance.

Our observations suggest that higher expression of Su(var)s, which promote a repressive chromatin environment, would result in stronger TE-mediated epigenetic effects (e.g., species on the left). With mechanisms that are yet to be revealed (see Discussion), the stronger epigenetic effects of TEs would reduce individual fitness, resulting in stronger selection against TEs and thus their lower population frequencies. Under the assumptions that the rate of TE increase through transposition is similar across species and the changes in TE copy number is at equilibrium, genomic TE abundance is determined by the strength of selection against TEs. Accordingly, the stronger epigenetic effects of TEs and the associated stronger selection removing them could drive an overall lower genomic TE abundance (e.g., species on the left). Observations made in this study are denoted with *.

It is worth noting that teasing apart the causality of the observed associations between the epigenetic effects and genomic abundance of TEs is challenging. A reduction in effective population size and accordingly weakened effectiveness of natural selection at removing TEs could drive the accumulation of TEs in some of the species (e.g., Mérel et al., 2021). Under this scenario, selection may favor epigenetic regulation that limits excessive inadvertent spreading of repressive marks from euchromatic TEs, which could also result in our observed weaker TE-mediated epigenetic effects in species with abundant TEs. Moreover, even if the difference in genomic TE abundance is indeed the ‘consequence’, varying efficacy of selection in purging deleterious TE insertions could also contribute to the observed between-species difference in TE abundance. Specifically, within species complex, estimated nucleotide polymorphism, an indicator for effective population size (Charlesworth, 2009), largely follows a similar rank order to our observed magnitude of TE-mediated H3K9me2 enrichment (nucleotide polymorphism, melanogaster complex: D. simulans > D. mauritiana > D. melanogaster (Langley et al., 2012; Meiklejohn et al., 2018) and yakuba complex: D. tessieri > D. yakuba > D. santomea (Bachtrog et al., 2006); the only difference in the rank order being the epigenetic effect of TEs is stronger in D. mauritiana than in D. simulans, Figure 5A). More effective selection at removing TEs in species with larger effective population size (but also stronger epigenetic effects of TEs) is another plausible driver for varying TE abundance that our analysis could not rule out.

Our study focuses on late-stage embryos and would capture the averaged epigenetic effects of TEs across heterogeneous cell types at this developmental stage. If TE-mediated epigenetic effects vary minimally across tissues or cell types, our estimate serves as a good proxy to study the associated fitness consequence. Yet, if TE-mediated enrichment of repressive marks substantially varies between tissues and cell types that have different importance in shaping individual functions, our estimate could fail to capture the most relevant variation in TE-mediated epigenetic effects that determine individual fitness. Even though we did find evidence supporting selection against TE-mediated enrichment of repressive marks that was estimated from whole embryos, future investigation on the variability of such effects across tissues and cell types will be important for the finer dissection of the role of TE-mediated epigenetic effects in individual fitness and, thus, genomic TE abundance.

Proper maintenance of heterochromatin ensures genome function and integrity (reviewed in Janssen et al., 2018). Intriguingly, the expression level (see above), copy number (Levine et al., 2012; Ross et al., 2013; Helleu and Levine, 2018), and amino acid sequences (Levine et al., 2012; Sasaki et al., 2019) of heterochromatin structural and enzymatic components have been observed to vary substantially within and between species. The abundance and composition of heterochromatic repeats also rapidly turnover between genomes (Larracuente, 2014; Wei et al., 2014; Wei et al., 2018). While variation in these genetic factors was postulated to contribute to the evolution of heterochromatin functions (e.g., Ferree and Barbash, 2009; Helleu and Levine, 2018; Mills et al., 2019), our discoveries further extend their evolutionary impacts from the gene-poor genomic dark matter into the gene-rich euchromatin. Heterochromatin modulators may not only determine the local epigenetic impacts of euchromatic TEs on flanking sequences, but also shape TE abundance and, accordingly, the structure and function of the euchromatic genomes. The evolution of euchromatin and heterochromatin, two genomic compartments usually presumed to function independently, may be more interconnected than previously thought by the inadvertent deleterious epigenetic effects of a widespread genetic parasite.

Materials and methods

Drosophila strains

Request a detailed protocol

Drosophila strains used in this study include D. melanogaster ORw1118 (Sexton et al., 2012), D. simulans w501 (Drosophila Species Stock Center, strain 1) and Mod6 (strain 2), D. mauritiana w12, D. yakuba NY73PB (strain 1) and Tai 18E2 (strain 2), D. santomea AG01482, and D. teissieri GT53W gen14. Flies were cultured on standard medium at 25°C and 12 hr/12 hr light/dark cycles.

ChIP-seq and RNA-seq experiment

Request a detailed protocol

We collected 16–18 hr embryos from each strain to perform ChIP-seq and RNA-seq experiments (with two biological replicates). Before embryo collection, young adults (3–7 days of age) were allowed to lay eggs on fresh apple juice plates for 1 hr at 25°C. We then collected embryos for 2 hr on another fresh apple juice plates and stored those collection plates at 25°C to enrich for 16–18 hr embryos. Chromatin isolation and immunoprecipitation were performed following the modEncode protocol (http://www.modencode.org/) with the following modification. Before splitting the input and IP samples, we added 4 µl of SNAP-ChIP K-MetStat Panel (EpiCypher) to every 10 µg of chromatin. This concentration allowed ~200–500 barcoded H3K9me2 nucleosome reads with ~15 million Illumina reads, which is the targeted sequencing depth of our samples. The SNAP-ChIP K-MetStat Panel serves as a spike-in and allows normalization between samples. We used H3K9me2 antibody (abcam 1220), which was validated by the modEncode and shows H3K9me2-specific binding (Egelhofer et al., 2011), to perform the ChIP experiment. Quantifying ChIP-seq reads corresponding to different histone modifications of the SNAP-ChIP K-MetStat panel (see below) also found this antibody has high specificity to H3K9me2 (Figure 1—figure supplement 5). ChIP-seq libraries were prepared using NEBNext Ultra DNA Library Prep Kit for Illumina (NEB) following the manufacturer’s protocol. RNAs were extracted using the RNeasy Plus kit (Qiagen) and library prep with Illumina TrueSeq mRNA stranded kit (Illumina) following the manufacturer’s protocols. Both ChIP-seq and RNA-seq samples were sequenced on Illumina Hi-Seq4000 with 100 bp, paired-end reads.

PacBio assemblies

Request a detailed protocol

We generated PacBio assemblies for D. melanogaster ORw1118 strain with following procedures. High molecular weight DNA was extracted from 450 adults (mixed sexes) using Genomictip Kit (Qiagen) following the manufacturer’s protocol. Extracted DNA then underwent SMRTbell library preparation with target insert size as 20 kb and sequenced on one SMRTcell using P6-C4 chemistry on a Pacific Biosciences RSII platform. We used Canu (version 2.0, Koren et al., 2017) to assemble a draft genome from the PacBio reads. We then used pbmm2 SMRT Analysis (version 7.0.0) to align the PacBio raw reads to the genome and used Arrow, also from SMRT Analysis, to polish the genome. The final genome has N50=12.1 MB and 178 MB for the total genome size. PacBio assemblies for D. simulans and D. mauritiana are from Chakraborty et al., 2021, and for D. yakuba (NCBI: GCA_016746365.2 and GCA_016746335.2), D. santomea (NCBI: GCA_016746245.2), and D. teissieri (GCA_016746235.2) are from the Andolfatto lab (Columbia University).

Annotation of TEs in PacBio assemblies

Request a detailed protocol

We ran Repeatmodeler2 (version 2.0.1; Flynn et al., 2020) on PacBio genome assemblies. NCBI BLASTDB and the analyzed genomes were used as inputs for the BuildDatabase function in Repeatmodeler2. We ran Repeatmodeler2 with options ‘-engine ncbi -LTRStruct’ and selected the output repeats with class as LTR, LINE, DNA, or unknown as potential TE sequences. We required a TE to be at least 500 bp to be included in our analyses.

To assign family identity to identified TEs, we used an iterated blast approach. We first blasted the TE sequences (using blastn; Camacho et al., 2009) to TEs annotated in D. melanogaster reference genome (version 6.32) and ‘canonical’ TE sequences of Drosophila (retrieved from Flybase September 2019) with following parameters: -evalue 1e-5; -perc_identity 80 for the melanogaster complex and -perc_identity 60 for the yakuba complex. When an identified TE has blast hits to multiple TE families, we only assigned the TE to a family when at least 80% of the covered query belongs to one and only one family. We then added the annotated TEs to the ‘blast database’ and repeat the process three more times in order to allow the identification of TEs that are diverged from those in the reference D. melanogaster genome or canonical TEs annotated in various Drosophila species. TE insertions that are within 500 bp were merged if they are from the same TE family or excluded from the analysis if they belonged to different families. We excluded DINE-1, which are mostly fixed in the melanogaster complex species (Kapitonov and Jurka, 2003) but underwent a recent burst of activities in D. yakuba and likely other closely related species (Yang and Barbash, 2008). We also excluded telomeric TEs (HeT-A, TART, and TAHRE), which predominantly locate at the end of chromosomes that are largely heterochromatic. Because most of the TEs included in our blast database are from D. melanogaster, it is plausible that our family assignment process is biased against assigning family identity to TEs in the species of the yakuba complex. Accordingly, except for analyses that require TE family identity (e.g., Figure 2B, Figure 2C, Figure 2—figure supplements 1 and 2 and Figure 2—figure supplement 4, and the estimation of TE population frequencies, see below), we included all TEs that are at least 500 bp in the analysis, irrespective whether we could assign their family identities or not. To identify whether TEs were full-length or truncated, TE sequences were aligned to the annotated D. melanogaster canonical sequences for the respective families using MAFFT (v7.505, Katoh and Standley 2013). TEs whose length is at least 70% of the annotated D. melanogaster canonical sequences were considered intact while others were deemed truncated. This D. melanogaster-centric analysis may bias against assigning TEs in other species as full length.

Identification of euchromatin/heterochromatin boundaries

Request a detailed protocol

Because we are interested in the epigenetic effects of euchromatic TEs, our analysis excluded those that are in or close to heterochromatin. To determine the euchromatin-heterochromatin boundaries in each genome, we generated H3K9me2 fold-enrichment tracks using MACS v2.7.15 (Zhang et al., 2008). We then visualized the H3K9me2 enrichment genome-wide using IGV (version 2.10.2, Thorvaldsdóttir et al., 2013) to identify the sharp transition in H3K9me2 enrichment and used positions that are 0.5 Mb ‘inward’ (toward the euchromatin) from the transition as the euchromatin-heterochromatin boundaries. These conservative euchromatin-heterochromatin boundaries are expected to minimize the influence of constitutive heterochromatin in influencing our analysis.

ChIP-seq analysis

Request a detailed protocol

In order to align reads originated from the spike-in control (SNAP-ChIP K-MetStat Panel, Epicypher), we combined PacBio genomes with DNA sequence barcode of the SNAP-ChIP K-MetStat Panel, which were serve as the ‘reference genome’ for Illumina sequence alignment. Raw reads were trimmed using Trimmomatic v0.35 (Bolger et al., 2014) before aligning to the reference genome with bwa mem v 0.7.16a (Li and Durbin, 2009). Read with low mapping quality (q<30) and reads that mapped to multiple locations were removed using Samtools (Li, 2011). We used the abundance of reads corresponding to the spike-in H3K9me2 nucleosome to normalize the background level of H3K9me2, following Lam et al., 2019. Briefly, we calculated the enrichment of spike-in H3K9me2 nucleosome reads as: Esi = (barcode fragments in ChIP)/(barcode fragments in input). We then used bedtools v2.25.0 (Quinlan and Hall, 2010) to obtain read coverage across the genome for ChIP and input samples. For 25 bp nonoverlapping windows across the genome, we calculated the per-locus enrichment as Elocus = (fragment coverage in ChIP)/(fragment coverage in input). The HMD for H3K9me2 in a particular 25 bp window was then estimated as Elocus/Esi. Windows with fewer than five fragment coverage in input samples were treated as missing data.

Estimation of TE-mediated H3K9me2 enrichment

Request a detailed protocol

For each TE, we normalized the local H3K9me2 HMD level in its flanking sequence by the median HMD for regions 20–40 kb upstream and downstream of each TE, following Lee and Karpen, 2017. The reasons behind this approach come from the observations that TE-mediated spreading of repressive epigenetic marks is usually within 10 kb in Drosophila (Lee, 2015; Lee and Karpen, 2017). We then divided the 20 kb upstream and downstream from a TE into 1 kb nonoverlapping windows and, for each window, calculated the median of normalized H3K9me2 HMD among its 40 25bp-HMD units (see above, m-HMD); at least 10 HMD estimates are required to calculate m-HMD for a 1 kb window. To estimate the magnitude of TE-mediated local enrichment of H3K9me2, we calculated the m-HMD for the 1 kb left and right flanking regions separately and took the average of the two sides. To estimate the extent of H3K9me2 spreading, we examined whether m-HMD is above one, which indicates that the H3K9me2 HMD level for the window is higher than that of the local background. We scanned across windows, starting from those right next to TEs and identified the farthest window in which the m-HMD was consecutively above one. We then used the average for the two sides as the extent of H3K9me2 spreading. The estimates for HMD magnitude or extent from two replicates positively correlate (Spearman rank correlation coefficient ρ=0.19–0.84, p<10–8, Figure 1—figure supplement 6). Curiously, the strength of correlation for the estimated magnitude or extent of TE-mediated epigenetic effects among replicates is low for some samples (e.g., magnitude of the effect for D. melanogaster). We performed IDR (irreproducible rate) analysis between replicates (Li et al., 2011) and found limited associations between IDR and the strength of between-replicate correlation (Figure 1—figure supplement 7). For instance, D. melanogaster shows the lowest correlation for the magnitude of TE-mediated H3K9me2 enrichment between replicates, but has decent consistency between replicates by IDR analysis. A plausible cause for this discrepancy is that TE-mediated enrichment of repressive marks does not have the typical characteristic of ‘enrichment peaks’ (Lee and Karpen, 2017; Lee et al., 2020) and thus could be oftentimes missed by custom pipelines that were designed to identify ‘sharp peaks’ (reviewed in Park, 2009; Nakato and Shirahige, 2017). For instance, only 12% of TEs showing local enrichment of H3K9me2 was detected by peak calling in a previous study (Lee et al., 2020). On the other hand, IDR analysis estimates the reproducibility of called peaks between replicates and might not include the majorities of regions with TE-mediated enrichment of repressive marks. More importantly, the spreading of repressive marks from constitutive heterochromatin is a variegated phenotype and was observed to vary between cells and individuals (reviewed in Elgin and Reuter, 2013), which could explain the variability between replicates. While the rank order of significance for called peaks are mostly consistent between replicates and falls along the diagonal line for significant called peaks for most samples (Figure 1—figure supplement 7), we noticed that this trend in D. simulans strain 2 is weak. We thus tried excluding this sample in analysis that tested the associations between genomic TE abundance and TE-mediated epigenetic effects (Figure 5A) and reached similar conclusion (see text).

To proceed with following analysis, we averaged the estimates from two replicates to generate the magnitude and extent of H3K9me2 enrichment. It is worth noting that the estimated extent of H3K9me2 spreading with different thresholds of m-HMD strongly correlate (Spearman rank correlation coefficient ρ=0.69–0.87, p<10–16, Figure 1—figure supplement 8), suggesting the robustness of such estimate. In the analysis, a TE was considered showing ‘epigenetic effects’ when the magnitude of H3K9me2 enrichment is above one or has at least 1 kb spreading of H3K9me2 (see text). Because the magnitude and extent of H3K9me2 enrichment do not follow normal distribution even after log transformation (Anderson-Darling normality test, for all test p<2.2e–16), we chose to perform nonparametric Spearman rank correlation tests when investigating the associations between TE-mediated epigenetic effects and factors of interests. The nonparametric Spearman rank correlation tests are also less sensitive to the effects of outliers, which is critical for our goal of identifying genome-wide patterns.

For two species (D. simulans and D. yakuba), we also investigated the epigenetic effects of TEs by comparing the enrichment of H3K9me2 at homologous sequences with and without TE insertions and contrast that to single species-based method (see above). We first used minimap v2.17 (Li, 2018) to generate assembly-to-assembly alignments of the PacBio genomes and identify homologous sequences for the two strains of a species (paftools.js from minimap). We ran MACS v2.7.15 (Zhang et al., 2008) to identify shared peaks of H3K9me2 enrichment using liberal significant threshold (with broad-cutoff p=0.5). TEs in these shared peaks were then excluded from the analysis because we could not determine if the enrichment of H3K9me2 for these regions were induced by TEs. We also excluded focal TEs whose homologous sequences were within 1 kb of another TE in the alternative strain. The magnitude of H3K9me2 enrichment was estimated as the m-HMD of the focal TE in the focal strain standardized by the m-HMD in the 1 kb TE-flanking regions in the homologous sequence in the alternative strain. The extent of H3K9me2 spreading was measured as the distance for the farthest windows from TEs that the m-HMD was consecutively higher in focal strain than that in the alternative strain. The estimates based on one or two genomes significantly correlate (Figure 1—figure supplement 1; see text).

Inference of genic H3K9me2 enrichment

Request a detailed protocol

Reference genome sequences and annotations for D. melanogaster, D. simulans, D. mauritiana, and D. yakuba were downloaded from NCBI Datasets (D. melanogaster, version Release 6 plus ISO1 MT; D. simulans, version ASM75419v2; D. mauritiana, version ASM438214v1; D. yakuba, version dyak_caf1). The annotations were lifted to PacBio assemblies by aligning the PacBio assemblies to NCBI references using minimap v2.17 (Li, 2018). Because the annotations for D. santomea and D. teissieri were not available when we performed the analyses, we used MAKER v2.31.8 (Holt and Yandell, 2011, p. 2) to annotate the PacBio assemblies. Specifically, we used Trinity-v2.85 (Grabherr et al., 2011) to de novo assemble the mRNA-seq for the genome as EST evidence for MAKER. We also supplied the protein sequences from D. melanogaster, D. sechellia, D. simulans, D. mauritiana, D. erecta, and D. yakuba as protein homology evidence for MAKER. We then ran MAKER with the default settings to obtain annotations for D. santomea and D. teissieri. To study the relationships between the epigenetic effects of TEs and the epigenetic states of neighboring genes, we assigned each TE to its closest gene and distinguished whether it inserted at the 5’ or 3’ side of the gene. Genic H3K9me2 enrichment is estimated as the average of HMD of the gene body, excluding genes with TEs inserted. When comparing the H3K9me2 enrichment level of homologous genic alleles with and without adjacent TEs, we calculated z-score as: (mean HMD of allele with nearby TE – mean HMD of allele without nearby TE in the alternative strain)/(standard deviation of both strains). Because genic H3K9me2 enrichment does not follow normal distribution, even after log transformation (Anderson-Darling normality test, for all test p<2.2e–16), we chose to perform nonparametric Spearman rank correlation tests when investigating the associations between genic HMD and factors of interests.

Association between TE-mediated epigenetic effects and TE abundance across species

Request a detailed protocol

To examine whether the epigenetic effects of TEs associate with TE abundance across species while accounting for species complex effects, we performed linear regression analyses using the following model: TE abundance ~ species complex (melanogaster or yakuba complex) + epigenetic effect. ANOVA F-test was used to examine whether the epigenetic effect was significant. We also performed PGLS (Grafen, 1989; Martins and Hansen, 1997) analysis using a tree adapted from Turissini and Matute, 2017; Chakraborty et al., 2019, with arbitrary branch lengths: ((((Dsim_strain1:0.1,Dsim_strain2:0.1):0.15, Dmau:0.25):3,Dmel:3.25):7.25, (((Dyak_strain1:0.1,Dyak_strain2:0.1):0.9,Dsan:1):1.75, Dtei:2.75):7.75). Although the significance level was sensitive to the within-species branch lengths, the sign of the coefficients remained unchanged. To fit the above model, we used gls function from nlme package with a Brownian correlation structure based on the tree (imported using ape package) in R.

Gene expression analysis

Request a detailed protocol

To estimate gene expression abundance, we mapped the raw RNAseq reads to the annotated genomes with STAR v2.6.0 (Dobin et al., 2013) with options --quantMode TranscriptomeSAM GeneCounts --chimFilter None. In order to compare the expression levels between different species, we used reciprocal best blasts between D. melanogaster and one other species to identify one-to-one orthologs using blastn (version 2.8.1). We then obtained a set of shared orthologs among six species to compare expression levels. For this set of shared orthologs, we estimated the RPKM (reads per kilobase per million reads) as the averaged RPKMs from two replicates (RPKMs from two replicates strongly correlate; Spearman rank correlation coefficient ρ>0.98 for all genomes, p<10–22). Genes with 0 RPKM in both replicates were excluded from the analysis. We then ranked genes from the highest to lowest RPKMs in each strain to get expression rank.

Estimation of TE population frequencies

Request a detailed protocol

Raw Illumina reads for D. melanogaster (Lack et al., 2015), D. simulans (Rogers et al., 2014), D. mauritiana (Garrigan et al., 2012), and D. yakuba (Rogers et al., 2014) were downloaded from SRA (SRP006733, SRP040290, SRP012053, and SRP029453, respectively). We used FastQC v0.11.7 (https://qubeshub.org/resources/fastqc) to check the read quality and TrimGalore v0.6.0 (Babraham Bioinformatics, 2022) to remove adapter and low-quality sequences. Illumina reads were then mapped to the corresponding PacBio assembly using bwa mem v0.7.16a (pair-end mode). We further used samtools to filter out reads with mapping quality smaller than 50 (MAPQ < 50).

To call the presence/absence of TEs, we followed the basic ideas developed in Cridland et al., 2013; Lee and Karpen, 2017, with following modifications. Briefly, we parsed out reads that uniquely mapped to the ±500 bp around TEs annotated in the PacBio assembly using seqtk v1.3-r107-dirty (https://github.com/lh3/seqtk; Li, 2022). Parsed reads were assembled into contigs using phrap v1.090518 (Ewing and Green, 1998) with parameters from Cridland et al., 2013. We mapped the resultant contig to PacBio assembly using bwa mem (default options). If a single contig mapped across 20 bp upstream and downstream of an annotated TE, a TE is called absent. To determine whether a TE is present, we evaluated whether two separate contigs spanned across the start and end of a TE’s boundaries, respectively, with at least 30 bp inside the TE and 20 bp outside the TE and a minimum total alignment length of 50 bp. We also evaluated contigs that mapped within annotated TEs. A TE is called present if there are two contigs spanning both the start/end of the TE insertion respectively or there is one contig spanning either the start or end of the TE and one contig aligned within TE insertion. A TE is considered missing data if none of the above criteria were met.

In the initial runs, we noticed that the failure to identify some TEs in the tested strain is due to the imprecise TE boundaries annotated by RepeatModeler2. In order to reduce the rates of missing data, we refined the called TE boundaries with following procedures. TE absence calls could be viewed as deletion structural variants (SVs) with respect to the PacBio genome. We thus used Lumpy-sv v0.2.13 (Layer et al., 2014) to call SVs in each strain. Deletions that overlap with the annotated TEs were extracted using sytyper v0.0.4 (Chiang et al., 2015) and merged using svtools v0.5.1 (Larson et al., 2019) to refine TE boundaries, which were then used in the above pipeline. For most TEs, the differences between the updated and RepeatModeler2 boundaries are short and only 5.2–7.9% of TEs boundaries in each genome were significantly updated (≥20 bp from original boundaries). Code for the TE calling pipeline can be found at https://github.com/harsh-shukla/TE_freq_analysis, (Huang, 2022 copy archived at swh:1:rev:24218fab83996f657e489402c5ff1c4cc06bfe9c).

Quantification of the heterochromatic repeats

Request a detailed protocol

To identify the repeat sequences enriched in the heterochromatic regions of the genome, we first used KMC (version 3.1.1, Kokot et al., 2017) to quantify the 12-mers in our H3K9me2 IP and matching input samples. In order to compare 12-mers abundance between IP and input libraries, we normalized the 12-mers counts by the number of reads mapped uniquely to the PacBio genome with at least 30 mapping quality score. 12-mers that have at least a threefold enrichment in an IP sample when compared to its matching input sample were considered as heterochromatic repeats.

Because PCR amplification of sequencing libraries were shown to influence the quantification of simple repeats (Wei et al., 2018), we sequenced the genomes of focused strains with Illumina PCR-free library preparation. We extracted DNA with 40 females using DNeasy Blood & Tissue Kit (Qiagen), following the manufacturer’s protocol. Extracted DNA was then prepared into Illumina sequencing libraries with PCR-free protocol and sequenced with 150pb paired-end reads by Novogene (Sacramento, CA). We then ran KMC3 on these libraries to quantify the abundance of heterochromatic repeats identified above. In order to compare across strains, the number of heterochromatic repeats were further normalized with the total number of reads from either the orthologous region.

Drosophila mutant crosses and eye pigmentation assay

Request a detailed protocol

To investigate the mutant effects of identified candidate Su(var)s, we followed approaches outlined in Sentmanat and Elgin, 2012. Specifically, we used a strain developed in Sentmanat and Elgin, 2012, which has a TE 1360 placed next to mini-white in y; w background. It was found that 1360-mediated epigenetic effects influence the expression level of mini-white and thus the intensity of eye pigmentation. Our analysis only tested the effects of Su(var)s whose existing mutants do not have any eye markers to avoid confounding effects on the quantification of eye pigmentation. We crossed 3- to 5-day-old virgin females of this strain to males of the Su(var) mutant strains or a control strain and then quantified the eye pigmentation level in 25–30 three- to five-day-old F1 males following methods described in Sentmanat and Elgin, 2012. Three independent crosses were performed for each mutant. Strains used in the analysis include BDSC 6398 (Hsc70-4 mutant), BDSC 11537 (RpLP0 mutant), BDSC 30640 (Bin 3 mutant), BDSC 36511 (Ago 2 mutant), and BDSC 6559 (y; w, control).

Data availability

Short-read data for ChIP-seq, RNA-seq, and Illumina whole-genome sequencing have been deposited to NCBI BioProject under accession number PRJNA855483. PacBio long-read data and genome assembly of D. melanogaster ORw1118 strain have been deposited to NCBI BioProject under accession number PRJNA855235.

The following data sets were generated
    1. Lee YCG
    (2022) NCBI BioProject
    ID PRJNA855235. Evolutionary causes and consequences of transposable elements' epigenetic effects.
    1. Lee YCG
    (2022) NCBI BioProject
    ID PRJNA855483. Evolutionary causes and consequences of transponsable elements' epigenetic effects.
The following previously published data sets were used
    1. Lack JB
    (2015) NCBI
    ID SRP006733. "1000" Drosophila Genomes - DPGP3.
    1. Rogers RL
    (2014) NCBI
    ID SRP040290. Drosophila simulans strain:multiple Genome sequencing.
    1. Rogers RL
    (2014) NCBI
    ID SRP029453. Drosophila yakuba strain:multiple Genome sequencing.

References

    1. Clark AG
    2. Eisen MB
    3. Smith DR
    4. Bergman CM
    5. Oliver B
    6. Markow TA
    7. Kaufman TC
    8. Kellis M
    9. Gelbart W
    10. Iyer VN
    11. Pollard DA
    12. Sackton TB
    13. Larracuente AM
    14. Singh ND
    15. Abad JP
    16. Abt DN
    17. Adryan B
    18. Aguade M
    19. Akashi H
    20. Anderson WW
    21. Aquadro CF
    22. Ardell DH
    23. Arguello R
    24. Artieri CG
    25. Barbash DA
    26. Barker D
    27. Barsanti P
    28. Batterham P
    29. Batzoglou S
    30. Begun D
    31. Bhutkar A
    32. Blanco E
    33. Bosak SA
    34. Bradley RK
    35. Brand AD
    36. Brent MR
    37. Brooks AN
    38. Brown RH
    39. Butlin RK
    40. Caggese C
    41. Calvi BR
    42. Bernardo de Carvalho A
    43. Caspi A
    44. Castrezana S
    45. Celniker SE
    46. Chang JL
    47. Chapple C
    48. Chatterji S
    49. Chinwalla A
    50. Civetta A
    51. Clifton SW
    52. Comeron JM
    53. Costello JC
    54. Coyne JA
    55. Daub J
    56. David RG
    57. Delcher AL
    58. Delehaunty K
    59. Do CB
    60. Ebling H
    61. Edwards K
    62. Eickbush T
    63. Evans JD
    64. Filipski A
    65. Findeiss S
    66. Freyhult E
    67. Fulton L
    68. Fulton R
    69. Garcia ACL
    70. Gardiner A
    71. Garfield DA
    72. Garvin BE
    73. Gibson G
    74. Gilbert D
    75. Gnerre S
    76. Godfrey J
    77. Good R
    78. Gotea V
    79. Gravely B
    80. Greenberg AJ
    81. Griffiths-Jones S
    82. Gross S
    83. Guigo R
    84. Gustafson EA
    85. Haerty W
    86. Hahn MW
    87. Halligan DL
    88. Halpern AL
    89. Halter GM
    90. Han MV
    91. Heger A
    92. Hillier L
    93. Hinrichs AS
    94. Holmes I
    95. Hoskins RA
    96. Hubisz MJ
    97. Hultmark D
    98. Huntley MA
    99. Jaffe DB
    100. Jagadeeshan S
    101. Jeck WR
    102. Johnson J
    103. Jones CD
    104. Jordan WC
    105. Karpen GH
    106. Kataoka E
    107. Keightley PD
    108. Kheradpour P
    109. Kirkness EF
    110. Koerich LB
    111. Kristiansen K
    112. Kudrna D
    113. Kulathinal RJ
    114. Kumar S
    115. Kwok R
    116. Lander E
    117. Langley CH
    118. Lapoint R
    119. Lazzaro BP
    120. Lee S-J
    121. Levesque L
    122. Li R
    123. Lin C-F
    124. Lin MF
    125. Lindblad-Toh K
    126. Llopart A
    127. Long M
    128. Low L
    129. Lozovsky E
    130. Lu J
    131. Luo M
    132. Machado CA
    133. Makalowski W
    134. Marzo M
    135. Matsuda M
    136. Matzkin L
    137. McAllister B
    138. McBride CS
    139. McKernan B
    140. McKernan K
    141. Mendez-Lago M
    142. Minx P
    143. Mollenhauer MU
    144. Montooth K
    145. Mount SM
    146. Mu X
    147. Myers E
    148. Negre B
    149. Newfeld S
    150. Nielsen R
    151. Noor MAF
    152. O’Grady P
    153. Pachter L
    154. Papaceit M
    155. Parisi MJ
    156. Parisi M
    157. Parts L
    158. Pedersen JS
    159. Pesole G
    160. Phillippy AM
    161. Ponting CP
    162. Pop M
    163. Porcelli D
    164. Powell JR
    165. Prohaska S
    166. Pruitt K
    167. Puig M
    168. Quesneville H
    169. Ram KR
    170. Rand D
    171. Rasmussen MD
    172. Reed LK
    173. Reenan R
    174. Reily A
    175. Remington KA
    176. Rieger TT
    177. Ritchie MG
    178. Robin C
    179. Rogers Y-H
    180. Rohde C
    181. Rozas J
    182. Rubenfield MJ
    183. Ruiz A
    184. Russo S
    185. Salzberg SL
    186. Sanchez-Gracia A
    187. Saranga DJ
    188. Sato H
    189. Schaeffer SW
    190. Schatz MC
    191. Schlenke T
    192. Schwartz R
    193. Segarra C
    194. Singh RS
    195. Sirot L
    196. Sirota M
    197. Sisneros NB
    198. Smith CD
    199. Smith TF
    200. Spieth J
    201. Stage DE
    202. Stark A
    203. Stephan W
    204. Strausberg RL
    205. Strempel S
    206. Sturgill D
    207. Sutton G
    208. Sutton GG
    209. Tao W
    210. Teichmann S
    211. Tobari YN
    212. Tomimura Y
    213. Tsolas JM
    214. Valente VLS
    215. Venter E
    216. Venter JC
    217. Vicario S
    218. Vieira FG
    219. Vilella AJ
    220. Villasante A
    221. Walenz B
    222. Wang J
    223. Wasserman M
    224. Watts T
    225. Wilson D
    226. Wilson RK
    227. Wing RA
    228. Wolfner MF
    229. Wong A
    230. Wong GK-S
    231. Wu C-I
    232. Wu G
    233. Yamamoto D
    234. Yang H-P
    235. Yang S-P
    236. Yorke JA
    237. Yoshida K
    238. Zdobnov E
    239. Zhang P
    240. Zhang Y
    241. Zimin AV
    242. Baldwin J
    243. Abdouelleil A
    244. Abdulkadir J
    245. Abebe A
    246. Abera B
    247. Abreu J
    248. Acer SC
    249. Aftuck L
    250. Alexander A
    251. An P
    252. Anderson E
    253. Anderson S
    254. Arachi H
    255. Azer M
    256. Bachantsang P
    257. Barry A
    258. Bayul T
    259. Berlin A
    260. Bessette D
    261. Bloom T
    262. Blye J
    263. Boguslavskiy L
    264. Bonnet C
    265. Boukhgalter B
    266. Bourzgui I
    267. Brown A
    268. Cahill P
    269. Channer S
    270. Cheshatsang Y
    271. Chuda L
    272. Citroen M
    273. Collymore A
    274. Cooke P
    275. Costello M
    276. D’Aco K
    277. Daza R
    278. De Haan G
    279. DeGray S
    280. DeMaso C
    281. Dhargay N
    282. Dooley K
    283. Dooley E
    284. Doricent M
    285. Dorje P
    286. Dorjee K
    287. Dupes A
    288. Elong R
    289. Falk J
    290. Farina A
    291. Faro S
    292. Ferguson D
    293. Fisher S
    294. Foley CD
    295. Franke A
    296. Friedrich D
    297. Gadbois L
    298. Gearin G
    299. Gearin CR
    300. Giannoukos G
    301. Goode T
    302. Graham J
    303. Grandbois E
    304. Grewal S
    305. Gyaltsen K
    306. Hafez N
    307. Hagos B
    308. Hall J
    309. Henson C
    310. Hollinger A
    311. Honan T
    312. Huard MD
    313. Hughes L
    314. Hurhula B
    315. Husby ME
    316. Kamat A
    317. Kanga B
    318. Kashin S
    319. Khazanovich D
    320. Kisner P
    321. Lance K
    322. Lara M
    323. Lee W
    324. Lennon N
    325. Letendre F
    326. LeVine R
    327. Lipovsky A
    328. Liu X
    329. Liu J
    330. Liu S
    331. Lokyitsang T
    332. Lokyitsang Y
    333. Lubonja R
    334. Lui A
    335. MacDonald P
    336. Magnisalis V
    337. Maru K
    338. Matthews C
    339. McCusker W
    340. McDonough S
    341. Mehta T
    342. Meldrim J
    343. Meneus L
    344. Mihai O
    345. Mihalev A
    346. Mihova T
    347. Mittelman R
    348. Mlenga V
    349. Montmayeur A
    350. Mulrain L
    351. Navidi A
    352. Naylor J
    353. Negash T
    354. Nguyen T
    355. Nguyen N
    356. Nicol R
    357. Norbu C
    358. Norbu N
    359. Novod N
    360. O’Neill B
    361. Osman S
    362. Markiewicz E
    363. Oyono OL
    364. Patti C
    365. Phunkhang P
    366. Pierre F
    367. Priest M
    368. Raghuraman S
    369. Rege F
    370. Reyes R
    371. Rise C
    372. Rogov P
    373. Ross K
    374. Ryan E
    375. Settipalli S
    376. Shea T
    377. Sherpa N
    378. Shi L
    379. Shih D
    380. Sparrow T
    381. Spaulding J
    382. Stalker J
    383. Stange-Thomann N
    384. Stavropoulos S
    385. Stone C
    386. Strader C
    387. Tesfaye S
    388. Thomson T
    389. Thoulutsang Y
    390. Thoulutsang D
    391. Topham K
    392. Topping I
    393. Tsamla T
    394. Vassiliev H
    395. Vo A
    396. Wangchuk T
    397. Wangdi T
    398. Weiand M
    399. Wilkinson J
    400. Wilson A
    401. Yadav S
    402. Young G
    403. Yu Q
    404. Zembek L
    405. Zhong D
    406. Zimmer A
    407. Zwirko Z
    408. Jaffe DB
    409. Alvarez P
    410. Brockman W
    411. Butler J
    412. Chin C
    413. Gnerre S
    414. Grabherr M
    415. Kleber M
    416. Mauceli E
    417. MacCallum I
    418. Drosophila 12 Genomes Consortium
    (2007) Evolution of genes and genomes on the Drosophila phylogeny
    Nature 450:203–218.
    https://doi.org/10.1038/nature06341
    1. Elliott TA
    2. Gregory TR
    (2015) What’s in a genome? the C-value enigma and the evolution of eukaryotic genome content
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 370:20140331.
    https://doi.org/10.1098/rstb.2014.0331
    1. Ewing B
    2. Green P
    (1998)
    Base-calling of automated sequencer traces using phred. II. Error probabilities
    Genome Research 8:186–194.
    1. Finnegan DJ
    (1992) Transposable elements
    Current Opinion in Genetics & Development 2:861–867.
    https://doi.org/10.1016/s0959-437x(05)80108-x
    1. Grafen A
    (1989) The phylogenetic regression
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 326:119–157.
    https://doi.org/10.1098/rstb.1989.0106
    1. Hill WG
    2. Robertson A
    (1966)
    The effect of linkage on limits to artificial selection
    Genetical Research 8:269–294.
    1. Kent TV
    2. Uzunović J
    3. Wright SI
    (2017) Coevolution between transposable elements and recombination
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 372:20160458.
    https://doi.org/10.1098/rstb.2016.0458
    1. Lee YCG
    2. Langley CH
    (2010) Transposable elements in natural populations of Drosophila melanogaster
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 365:1219–1228.
    https://doi.org/10.1098/rstb.2009.0318
    1. Wright SI
    2. Schoen DJ
    (1999)
    Transposon dynamics and the breeding system
    Genetica 107:139–148.

Decision letter

  1. Magnus Nordborg
    Reviewing Editor; Gregor Mendel Institute, Austria
  2. Molly Przeworski
    Senior Editor; Columbia University, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

[Editors' note: this paper was reviewed by Review Commons.]

https://doi.org/10.7554/eLife.81567.sa1

Author response

1. General Statements

This study builds upon a previous study published in eLife (Lee and Karpen 2017), which established the importance of transposable element (TE)-mediated local enrichment of repressive marks in TE evolution in the model species, D. melanogaster. In the current study, we investigated this TE-mediated phenomenon in six closely related Drosophila species to have a good phylogenetic resolution, aiming to address step-by-step how such TE-mediated effects may shape the evolution of genomic TE content, a long-standing question in evolutionary genomics. Our analyses successfully connect species-specific host chromatin regulation, TE-mediated epigenetic effects, the strength of natural selection against TEs, and genomic TE abundance unique to individual species, while uncovering the evolutionary causes for the wide variety of TE profiles between species.

We appreciate the reviewers for their overall positive response and helpful comments on the previous versions of the manuscript. We appreciate the opportunity to correct our previous oversights, revise analysis to corroborate findings, and improve interpretations of our observations. Thanks for these suggestions, we believe the revised manuscript is much improved and hope it is now acceptable for publication in eLife. Please see below for our pointby-point responses to reviewers’ comments.

In addition to analyses suggested by the reviewers, we also updated one of our analyses to avoid non-independence of samples. Specifically, in Figures 3A and B and Supplementary Figures S1-6, we investigated the associations between TE-mediated epigenetic effects and genic epigenetic states/expression. In our previous analysis, a TE could be associated with more than one gene, if it is the nearest TE to more than one gene. Yet, this could lead to nonindependence of observations. In this revised manuscript, we restricted the analysis to TE and its nearest gene, excluding this potential issue. The results are qualitatively the same and quantitatively similar, and our conclusions from these results remain unchanged.

2. Point-by-point description of the revisions

Reviewer #1 (Evidence, reproducibility and clarity (Required)):

Summary:

In this work titled "Species-specific chromatin landscape determines how transposable elements shape genome evolution", Huang, Shukla and Lee examined the epigenetic effect of euchromatic TE insertions in 6 Drosophila species (3 in melanogaster species complex and 3 in yakuba species complex). They profiled the H3K9me2 mark and transcriptome from 16-18hr embryos and, using genomics, population genetics and statistical tools, they extended claims made previously in D. mel to 5 more species: (1) H3K9me2 spreads to nearby regions, exerting an epigenetic effect; (2) TEs with epigenetic effects are selected against. More importantly, this work aims to identify rules that explain the differences in TE-mediated epigenetic effects across the 6 species examined and how their distinct chromatin landscape might affect the way in which TE-mediated epigenetic effects shape genome evolution. The analysis seems to be well done, rather rigorous, and considers important counterarguments in several places throughout the manuscript. The figure is clear, and the method is well detailed. I have some comments that should be addressed before publication, involving additional bioinformatic analysis and text edits to either tone down or re-frame the conclusions.

Major comments:

1. Exclusive usage of Spearman rank correlation for key conclusions requires further support. I worry that using Spearman rank correlation loses too much of the quantitative information in e.g., gene expression. Have authors ever looked at Pearson's correlation (perhaps upon log-transformed levels when necessary) and does it lead to similar or different results? While this concern holds for all conclusions based on Spearman rank correlation, it is particularly a concern for analysis involving gene expression, especially for the surprising claim that TE-mediated K9 spreading does not affect adjacent gene expression.

We agree with the reviewer that, compared to parametric tests (e.g., Pearson tests), nonparametric tests have more restricted statistical power. An insignificant Spearman correlation test could either be due to true biological reasons or simply a lack of statistical power. Yet, we found that the magnitude and extent of TE-mediated H3K9me2 enrichment and the genic HMD do not follow normal distribution even after log transformation (Anderson-Darling normality test, all p-value < 10-16). Accordingly, we are unable to perform the suggested Pearson correlation tests. Another reason why we would like to stick with the Spearman tests is because it is less sensitive to outliers than the Pearson correlation tests. This is especially important because our goal is to identify genome-wide average patterns.

It is worth noting that, in addition to Spearman rank correlation analysis, we also compared the expression level between homologous alleles with and without TEs showing H3K9me2 spreading in D. simulans and D. yakuba (Figure 3C and Figure 3 – Supplementary Figure S7-9) and found consistent results. We included additional discussion to explain why we chose to use nonparametric tests (p50, line 12-18; p51, line 27-31) and emphasize that, in addition to the correlation analysis, other tests also support limited impacts of TE-mediated H3K9me2 enrichment on gene expression (p39, line 13-15).

We agree with the reviewer that the observation of the limited associations between TE-mediated H3K9me2 spreading and gene expression is surprising. We hope our addressing the reviewer’s comment 3 and additional discussions and citations of other relevant/suggestive studies that reached similar conclusions (p39, line 15-22) help address this issue.

2. Choice of binning procedure over correlation analysis seems unjustified. I understand the binary binning in e.g., Figure 1C, but for quantitative variables e.g., distance to TE/gene, allele frequency in Figures3 and 4, I worry that the binning procedure based on somewhat arbitrary cutoffs could mask the real relationship. Are the conclusions robust against the cutoff chosen?

Is there actual correlation between the two variables examined?

For previous Figure 3A and 3B (now Figure 3 – Supplementary Figures 1-4), following the reviewer’s suggestion, we performed regression analysis that uses the distance between genes and TEs as continuous variables to test whether the associations between TE-mediated epigenetic effects and genic H3K9me2 enrichment/expression depends on TE-gene distance. Such analysis reached the same conclusion and is now included in p15, line31 – p16, line 4; p16, line 16-19. We chose to present the analyses that bin genes based on their distance to TEs because it is challenging to interpret them as three-dimensional graphs; current Figure 3 – Supplementary Figures 1-4 have TE-mediated epigenetic effects as X-axis, genic attribute as Y-axis, and bins of gene-TE distance as a binary color. It is worth noting that, in this revised version of the manuscript, we chose to present the correlation coefficients between TE-mediated H3K9me2 enrichment and genic epigenetic states/expression (new Figures 3A and 3B). This is because these figures provide a good visual summary for all the comparisons (magnitude/extent, genic H3K9me2 enrichment/expression, and TEs 5’ and 3’ to genes) and are easier to be interpreted than original Figures 3A and 3B (now Figure 3 – Supplementary Figures 1-4).

In Figure 3D, we estimated the z-score to compare the extent of TE-mediated spreading of H3K9me2 on the genic and intergenic sides between TEs close and distant from genes (the binning is based on 50% percentile of the distance between TEs and genes). Following the reviewer’s suggestion, we performed correlation tests between the z-score and TE-gene distance and found consistent results – we identified significant negative associations between the z-score of a TE and its distance from the nearest gene. These results are now presented in p19, line18-23.

For the population frequencies analysis (Figure 4), we were unclear about how we categorized low-frequency and high-frequency TEs. A TE is “low-frequency” if it is only identified in the focused strain, but not in the population (in other words, a private TE insertion). This information is now included in p23, line 19-21. Surveys on TE population frequencies in Drosophila found highly skewed frequency spectra and that the majority of the TEs appear as singletons (i.e., present in only one individual in the population sample); analysis typically categorizes TEs according to whether they are singleton (low-frequency in our analysis) or not (high-frequency in our analysis). Accordingly, we stick with this binning in our analysis.

3.1. TE position relative to genes seems oversimplified. The manuscript compares TEs that are close to or far from genes, but what happens to TEs in the genes (either in UTRs, CDSs or introns)? Do these different TE positions inside and outside genes exert different effects on K9 magnitude, K9 spreading extent, gene expression and/or allele frequency in the population?

(We split the reviewer’s comment #3 into two and answered them individually).

Following the reviewer’s suggestion, we compared the magnitude and extent of TE-mediated H3K9me2 enrichment between intergenic, intronic, and exonic TEs. While we don’t find any consistent significant differences for the magnitude of TE-mediated epigenetic effects, intergenic TEs exert a larger extent of H3K9me2 enrichment than other TEs in multiple genomes. These results are now included in p13, line 4-15 and Figure 2 – Supplementary Figure S4.

Our goal is to identify how TE-mediated epigenetic effects associate with genic epigenetic status, gene expression, and TE population frequencies. The difference in TE-mediated epigenetic effects between intergenic and genic TEs would have been captured by the estimated magnitude and extent of H3K9me2 enrichment. Accordingly, for the downstream analysis, we did not further distinguish between genic and intergenic TEs with one exception. Because TEs inside genes could influence gene expression through mechanisms other than altering the epigenetic states, a major confounding factor of previous studies (reviewed in (Choi and Lee 2020 PLoS Genetics, Kelleher et al. 2020 Trends in Genetics)), we excluded TEs inside genes when testing the associations between TE-mediated epigenetic effects and gene expression. We further emphasized this point in p15, line 15-18.

3.2. Are there any differences to TEs being 5' or 3' to the genes, since TEs 5' to the genes are more likely to affect the promoter? I worry that simple comparison of close vs. far might skew the actual relationships.

We appreciate the reviewer for pointing this out. Our additional analysis following the reviewer’s suggestions led to some new observations that further corroborate our conclusions, mainly whether neighboring gene expression could influence the extent of TE-mediated epigenetic effects.

We agree with the reviewer that TE-mediated epigenetic effects on the epigenetic states and expression of genes could differ when TEs are 5’ or 3’ to genes, especially given that TEs 5’ to genes should be more likely to influence the promoter. We performed correlation analyses between TE’s epigenetic effects and genic H3K9me2 enrichment or expression separately for TEs 5’ and 3’ to genes and reached similar conclusions to those of analyzing all TEs together (Figure 3B). Specifically, we found significant associations between TE’s epigenetic effects and genic H3K9me2 enrichment for both TEs 5’ and 3’ to genes, such associations depend on TE-gene distance, and we did not notice any major difference between TEs 5’ and 3’ to genes. On the other hand, we observed limited associations between TE-mediated epigenetic effects and gene expression, except for one incidence. These results are included in p16, line 23-p15, line 16 and Figure 3B.

We also tested whether the impact of adjacent gene expression on the extent of TE-mediated enrichment of H3K9me2 differs between TEs 5’ or 3’ to genes. We found that the difference in the extent of TE-mediated H3K9me2 enrichment is greater for TEs 5’ than 3’ to genes, likely due to the effects mentioned by the reviewer and providing stronger support for the hypothesized effects of gene expression on the extent of TE-mediated local H3K9me2 enrichment. We included these results in p20, line 5-10, Figure 3 – Supplementary Figure S12.

4. Transcriptional activity of TEs was not considered. Many TEs are (5'-) truncated and thus transcriptionally inactive. There is literature that suggests K9 installation and spreading only happen on transcriptionally active TEs. Have authors considered separating active/inactive or intact/truncated TEs? Combining TEs with different transcriptional activity (or K9 biology) despite the same length, family or position could confound the analysis and lead to different conclusions.

We agree with the reviewer that this is an important axis of the analysis that we did not consider. We categorized TEs whose length is shorter than 70% of the annotated, canonical sequence in D. melanogaster as truncated, and otherwise as full length. Consistent with what is suggested by the reviewer, full-length TEs tend to show a larger magnitude and extent of H3K9me2 enrichment. These results are now included in p11, line 18-26 and Figure 1 – Supplementary Figure S4. Because our following analysis focuses on how TE-mediated epigenetic effects, instead of individual TE attributes, influence genic H3K9me2

enrichment/expression and selection, we did not further consider full-length/truncated TEs for the downstream analysis.

5. Usage of 16-18hr whole-embryo averaged transcriptome and H3K9me2 epigenome require justification. The 16-18hr embryo is presumably consisted of a very heterogenous population of cells, how does such a population average reflect bona fide biology? Authors should either justify this or comment explicitly in text on potential problems such datasets could have on major conclusions.

We agree with the reviewer that 16-18hr embryos consist of heterogeneous cell types, which are precursors of different tissues and organs of a whole animal. By examining the epigenetic effects of TEs in whole embryos, our analysis would capture the averaged effects across cell types. If the relative importance of tissues/organs in determining individual fitness is proportional to the abundance of their specific cell precursors, our estimate could provide a reasonable estimate of the epigenetic effects of TEs that may be relevant to individual performance and thus fitness. Yet, if the fitness of an individual is predominantly determined by a tissue/organ that originated from a rare cell type and TE-mediated epigenetic effects vary substantially between cell types, our analysis of whole embryos may not be a good indication of the associated fitness consequence. Even though we did find evidence supporting selection against TE-mediated enrichment of repressive marks estimated from whole embryos, future investigation on the variability of such effects across tissues and cell types, which we aim to address in the future, will be important to further dissect the role of TE-mediated epigenetic effects in individual fitness and, thus, genomic TE abundance.

Following the reviewer’s suggestion, we emphasized that our dataset should be viewed as an average of all cell types and stated the potential caveats in p42, line 22-p43, lin2.

6. Euchromatic TE and genome-wide TE abundance should not be used interchangeably. This work is done solely on euchromatic TEs, but throughout the text, authors comment on how a factor examined might shape genome-wide TE abundance or even genome evolution. I wonder how much euchromatic TE copy number reflects genome-wide TE abundance? I'd expect heterochromatin harbors most TE copies and to be the major place that determines the overall TE abundance. Authors should either discuss how euchromatic TEs do reflect genome-wide TE abundance (and genome evolution) or change all related claims to better describe analysis done on euchromatic TEs.

We agree with the reviewer that this is a very important point that we did not specify previously. The genome-wide TE abundance we discussed in the text and also in the field mainly concerns euchromatic TEs. This is because of two major reasons. One is the technical difficulties associated with assembling and annotating TEs in the highly repetitive heterochromatic regions of the genome. The other reason is that, in the few species where heterochromatin has been assembled, most TEs in the heterochromatic regions are fragmented and have lost of transposition ability, making them no longer important players in the evolutionary dynamics of TEs (i.e., TE number increase through transposition and decrease through selection and excision). Following the reviewer's suggestions, we clearly defined that our study focuses on euchromatic TEs in the opening paragraph (p2, line 9-30), in the first paragraph of result (p7, line 21-22), and in the opening paragraph of discussion (p37, line 7-16).

We also emphasized that our study focused on euchromatic TEs and revised our claim at places throughout the manuscript (e.g., p6, line 8; p29, line 8, p37, line18).

7. Potential contradiction between key conclusion and the "species complex effect". In Figures5A and 5C, the species complex effect is obvious and, when we compare the two species complexes, one would conclude the opposite of what the authors attempted to argue. For instance, K9 and TE copy number correlate between mel and yak species complexes. Does this actually mean that there is different biology (or forces) at play at different evolutionary timescales? Authors' major claims seem valid within a shorter period of time (<3 million years) within each of the two species complexes, but they do not appear to hold when one looks at longer timescales. If so, the authors should modify key conclusions.

We appreciate the reviewer for pointing this out and agree with the reviewer’s interpretation. Indeed, according to our data, the importance of TE-mediated epigenetic effects in shaping genomic TE abundance mainly acts on a short evolutionary time scale.

Previous studies have found greater net gains of intron sequence and duplicated genes in D. yakuba than in D. melanogaster and/or D. simulans. These observations suggest that the D. yakuba genome may be more tolerant of gained sequence. If similar mechanisms apply to the other two species in the yakuba species complex, lineage-specific mechanisms that shape the evolution of the whole-genome may also play a vital role in determining genomic TE abundance, and might potentially override the influence of selection against TE-mediated epigenetic effects. We revised our conclusion (by emphasizing it supports association within species complex, p30, line 8-12 (in Results)) and included these discussions in and p41, line 16-29 (in Discussions).

8. Biological difference between H3K9me2 magnitude and spreading extent seems unclear. Roughly speaking, the beginning of the manuscript made the two seem similar, the middle part found an effect on the extent only, and the last part found an effect on the magnitude only.

While I appreciate the authors examining two aspects of H3K9me2 spreading, I am unsure what either aspect truly reflects. I worry that some readers would find the authors cherry-picking either aspect to support their claims – perhaps authors should've commented a priori the difference between the two aspects such that conclusions made on one, but not the other, aspect of TE-mediated epigenetic effects are more convincing.

We agree with the reviewer that we did not make clear distinctions between the magnitude and extent of H3K9me2 enrichment, which could confuse readers about why we would need both indexes.

Previous studies on the spreading of repressive marks from constitutive heterochromatin suggested that the magnitude and extent of H3K9me2 enrichment may not be perfectly associated, mainly because the local genomic context could influence the extent of the spread. If similar molecular mechanisms are also applicable to epigenetically silenced TEs in euchromatin, the magnitude and extent should provide different information about TE-mediated local enrichment of repressive marks. Before conducting this study, we had no a priori expectation for which index would be more relevant to the questions we aim to address. Accordingly, we feel it is most appropriate to present both of them. Especially, in several analyses where the results differ between the two indexes, they provide important information about the underlying mechanisms (see below).

According to our observations, the major difference between these two indexes is that the magnitude of TE-mediated H3K9me2 enrichment is more significantly associated with lower TE population frequencies (i.e., strong evidence supporting being selected against) and is correlated with genomic TE abundance. On the other hand, the impact of neighboring gene expression only influences the extent, but not the magnitude, of TE-mediated H3K9me2 enrichment. In the discussion, we interpreted these two findings together as the result of the extent being more sensitive to the local genomic context (which is similar to previous studies) than the magnitude, making the latter a more direct estimate for the functional and thus fitness impact (which we still don’t know, see p40, line 18-p41, line 2 for discussion) of TE-mediated H3K9me2 enrichment.

To make better distinctions between these two indexes, we (1) included definitions for what these two indexes are at where they were first mentioned in Results (instead of just in Materials and methods; p10, line 3-9), (2) discussed previous observations to explain why we used both indexes and provided our a priori guess that the extent may be more sensitive to local genomic context (p10, line 19-p11, line 8), (3) moved discussions regarding the difference and relative importance of these two indexes from the end of Discussion to the front of Discussion (p38, line 23-p39, line 7), and (4) emphasize their difference throughout the text when necessary. We hope these changes addressed the reviewer’s concerns.

9.1. Chromatin landscape was examined by Su(var) gene expression and satellite repeat abundance, which seems insufficient and incomplete. Chromatin landscape is too broad of a term when authors only examined Su(var) gene expression and simple repeat abundance. Have authors look at other heterochromatin genes (e.g., ZAD-ZNF, piRNA pathway) besides Su(var)?

(We split the reviewer’s comment #9 into two and answered them individually). Most of the piRNA pathway genes are expressed in the female germline and at the early embryonic stage. Accordingly, nearly none of them are expressed in our dataset, which was generated using 16-18hr late-stage embryos. In Drosophila, ZAD-ZNF proteins have a wide range of functions. For a ZAD-ZNF gene that has a demonstrated role in heterochromatin function (Oddjob, a Su(var)), we included it in our analysis.

It is worth noting that, instead of looking at the initiation of the epigenetic silencing of TEs, we mainly focus on the inadvertent side effect associated with the silencing – the enrichment of repressive marks at flanking non-TE sequences. This phenomenon is similar to the well-studied Position Effect Variegation (PEV). Accordingly, our analysis focuses on Su(var), whose role in determining the strength of PEV has been well characterized. To explain this point, we included more motivation for looking at Su(var) at p31, line 8-15. We agree with the reviewer and tone down several places that our study mainly concerns repressive chromatin landscape (e.g., p38, line 1) in addition to other existing places.

9.2 What is the relationship between simple repeat abundance, TE abundance, and genomewide repeat content? This latter point is related to the point #6 above about euchromatic TE vs. genome-wide TE. I kept wondering, in the context of this work, whether there is biological difference between satDNA and TE – are they similar because they cause similar epigenetic effects or do they differ?

As described in our response to reviewer’s comment #6, our study mainly concerns the evolutionary dynamics of TEs in the euchromatic genome. We hope changes made while addressing that comment clarify this point.

We are unaware of studies that report H3K9me2/3 enrichment mediated by satellites in the Drosophila euchromatic genomes, except for those extreme cases through experimental manipulation (e.g., a large tandem array of constructs). The replication/expansion mechanisms and potential fitness impacts also differ between TEs and satellite repeats. Accordingly, we would argue the evolutionary dynamics would differ between these two types of repetitive sequences.

Minor comments:

– Figure 2 beautifully shows the epigenetic effect depends on the TE family and other factors. I wonder if the authors would want to expand Figures3-5 to see if any TE family exerts a particularly strong (or weak) epigenetic effect in different species?

Because there are only a few TE families that have at least five copies included in the analysis for all genomes, our power to perform the suggested analysis is limited. Also, in Figures 3-5, our main goal is to identify the associations between the strength of TE-mediated epigenetic effects and various factors (e.g., gene expression, genomic TE abundance), irrespective of TE identity.

Due to these two reasons, we did not perform the analysis suggested by the reviewer.

Nevertheless, we hope the data presented in Figure 2C provides a glimpse into the reviewer’s question. We added some more discussions in the text regarding this p12, line 21-26, pointing to TE families that are of particularly strong/weak epigenetic effects across species.

– Figure 6 cartoon probably should have K9 spread less towards the gene side.

We agree with the reviewer and revised the figure accordingly.

Reviewer #1 (Significance (Required)):

The senior author has previously found that euchromatic TE exerts an epigenetic effect via H3K9me2 spreading, which confers a fitness cost. Now, this work builds on the previous finding to attempt a conceptual advance that such a TE-mediated epigenetic effect differs across six species examined and depends on species-specific chromatin landscape. The examination of six species and the connection between molecular biology, population genetics and potentially genome evolution are rare among existing literature and thus present novelty. It should interest people working on TE, heterochromatin, epigenetics and population genetics. My expertise is in chromatin, epigenetics and RNAi. I do not have sufficient expertise to evaluate the regression analysis and TE population frequency estimation done in this work.

Reviewer #2 (Evidence, reproducibility and clarity (Required)):

Summary: In this article, Huang and colleagues investigate the degree to which TEs can induce heterochromatin formation within the euchromatin across species. They hypothesize that different TE abundances across species might be explained by variation in the degree to which TEs impact a selective effects. In this article, they show that there does seem to be a big difference in the magnitude of H3K9me2 induction between the mel clade and the yak clade. This is extremely interesting and very worthy of publication. The authors do a great job exploring a variety of alternate hypotheses.

Major Comments:

I have no major comments that may significantly alter the conclusions. I found their approach robust. However, I will admit there is one thing I find difficult to understand regarding their model and conclusion. In the introduction they state:

"These observations spurred our previous hypothesis that varying epigenetic effects of TEs could result in between-species differences in the strength of selection against TEs, eventually contributing to divergent genomic TE abundance (Lee and Karpen 2017)."

And in the final figure (figure 6) they provide a schematic of this hypothesis (note typo in the worded "mediated" in the figure). In particular, they show (indicating the * where they state the paper found) that (1) The magnitude of Heterochromatin induction can lead to (2) lower TE population frequencies and this can lead to (3) lower genomic TE abundance. But doesn't 5A show that the yak clade has HIGHER magnitude of heterochromatin induction (as shown earlier) but also higher estimated TE copy number? So, 5A shows the opposite conclusion with respect to genomic copy number. I feel the disconnect might arise from mixing up the notions of low frequency of TE insertions from the notion of genomic copy number. It is possible to have a VERY large number of TEs that are segregating at very low numbers. So, selection against TEs (as evident by low allele frequency) doesn't easily translate to total copy number. I hope this makes sense.

We appreciate the reviewer for pointing this out and this is indeed something we should address. Part of the issue raised by the reviewer is hopefully addressed in our response to reviewer #1’s comment #7 (please see above). In short, the strong species complex effect suggests that the role of TE-mediated epigenetic effects in determining genomic TE abundance may be acting on a short evolutionary time scale and could be overridden by other general processes that determine gain/loss of sequences (irrespective of types). This may explain why TEs in yakuba complex species show stronger epigenetic effects but are also more abundant.

We agree with the reviewer that our model (presented in Figure 6) has unstated underlying assumptions. According to the population genetic theory of TE copy number, the equilibrium copy number is determined by the strength of selection removing TEs and TE increase through transposition rate. Our model assumes that the transposition rates are similar, and TE copy numbers are at equilibrium between species. Here, the population frequencies of TEs serve as an indicator for the strength of selection removing TEs. It is possible to have a large number of TEs segregating at low population frequencies, but the average copy number across individuals of a species (i.e., the number of TEs in each genome) would still be small. We hope this interpretation makes sense to the reviewer and are open to suggestions to further revise the model. We included additional discussions regarding the assumptions of our model, which can be found at p41, line 7-10 and in Figure 6 legend.

Minor Comments: I have some minor comments for clarity.

1) Since the paper is motivated to explain differences between species in genomic copy number, this information should be provided in the first figure in 1A (alongside the species name, with the two sim and yak samples also both provided numerically)

We revised Figure 1A according to the reviewer’s suggestion.

2) Something should be made of the fact that some of these arguments can be circular and the arrow of causality is hard to tease apart. Genomic TE abundance may be explained by species level differences in the epigenetic impact of TEs. But, it may also be the case that as TE abundance varies, there may be evolution of the epigenetic impact. Causality is tricky here and worth being up front about.

We agree with the reviewer, and we included this possibility in Discussion. In this revised manuscript, we rewrote the discussion to further emphasize this point (p42, line 1-8).

3) Some clarity in the text (not just in the methods) should be provided for the reader about how the "extent of spreading" was calculated.

We agree with the reviewer and added these details in the Results section when this index was first mentioned (p10, line 3-9).

4) In figure S1, what does it mean for points to have zero magnitude of enrichment but a great extent of enrichment? Shouldn't all points for zero magnitude be at the origin? I am looking at the scattered points that run along the Y axis in S1 of just the mel group (not seen in yak group)

We greatly appreciate the reviewer for pointing this out.

We identified an error in the coding for generating this plot. In this script, a few “NA”s in the data file were treated as 0 when calculating the average magnitude of the two replicates. This error affected the results of Figure 1- Supplementary Figure S1, Figure 2B and 2C, Figure 2 – Supplementary Figures S1 & S2. We updated these figures and revised the correlation coefficients in Figure 1 – Supplementary Figure S1. After removing these TEs with erroneous estimates, the correlation between the magnitude and extent of TE-mediated epigenetic effects becomes stronger. Other updated results reached the same conclusions as previous results.

We thoroughly checked our coding, and all the other results throughout the study are unaffected. For other analyses, we used a different script in R, which removes the “NA”s.

While addressing this comment, we noticed that we forgot to provide details about why a TE may have NA estimates for TE-mediated epigenetic effects (low local sequencing coverage) and how we treated them in the analysis (remove them from the analysis). We included these details in p48, line 20-21; p48, line 31.

5. In S2, how they compare the correlations across two genomes is confusing. Is that for shared TEs? Also, do estimates from two strains (Y axis) include values from X axis? It seems like this is saying X is correlated with the average of X and Y, which will always be true if Y is a random variable. Overall, I was a fairly confused about what S2 was showing.

We agreed with the reviewer that neither the main text nor the legend for Figure1 –

Supplementary Figure 1 (previous Figure1 – Supplementary Figure 2) contains sufficient details for how these two estimates were generated.

The concepts for estimating the magnitude and extent of TE-mediated epigenetic effects are similar for using one genome and two genomes. The major difference is the baseline used. For one genome, we used H3K9me2 enrichment > 1 as a criterion to estimate the extent. In contrast, for two genomes, the comparison is made between homologous sequences of the focal genome with the focused TE insertion to the alternative genome without the TE insertion (i.e., whether H3K9me2 enrichment in strain 1 > H3K9me2 enrichment in strain 2). For the magnitude of TE-mediated H3K9me2 enrichment, the two-genome estimates were further normalized by the enrichment level in the homologous sequence without the TE insertion. We agreed with the reviewer that these two estimates are related due to shared variables, but the relationship between them is not a simple arithmetic calculation. Accordingly, we feel that presenting this correlation result is still important.

We included some more details in the main text (p10, line 9-13) and Figure 1 – Supplementary Figure S1 legend that hopefully could help readers interpret the result. We also noted that detailed methods can be found in Materials and methods.

When revising according to this comment, we noticed that we did not include method details about how we identify homologous regions in the two PacBio genome assemblies. We added these details (using minimap) in p50, line 22-25.

6) In figure 2A. Do not includes lines connecting values between species. Just used colored dots This is not a series. I think this issue is elsewhere as well.

We appreciate the reviewer for pointing out our oversight and revised Figures 2A and 2C accordingly.

Reviewer #2 (Significance (Required)):

This article provides some extremely important insight into the manner in which heterochromatin induction can shape TE dynamics across species. What explains variation in genome size and TE content is one of the big questions in evolutionary genomics. The results are somewhat complex, but that is because biological systems are complex.

My expertise is in evolutionary genomics of TEs. I am not an expert in chromatin biology and regulation of gene expression.

Reviewer #3 (Evidence, reproducibility and clarity (Required)):

In this manuscript, Huang et al. expand upon results presented in two previous influential publications led by the corresponding author of this manuscript: Lee (2015) and Lee and Karpen (2017). Lee (2015) reported that H3K9me3 spreads from euchromatic TE insertions, these insertions lead to reduced expression of nearby genes, and this phenomenon depends on the piRNA pathway. Lee and Karpen (2017) reported a significant negative correlation between both H3K9me2 magnitude and extent and TE population frequencies in D. melanogaster and demonstrated strong spreading of H3K9me2 from TE insertions, an average of 4.5 kb into flanking regions. They also report stronger epigenetic effects of TEs in D. simulans compared to D. melanogaster which they propose could be due to higher expression of Su(var) 3-9.

Here, Huang et al. examine H3K9me2 magnitude and extent (i.e. spreading) across six species of Drosophila, including two different strains of D. simulans and D. yakuba. Across all species and strains, the authors observe a strong enrichment of H3K9me2 at euchromatic TE insertions. They find a large amount of variation in H3K9me2 magnitude and extent across TE classes and families, and among species and strains, but this variation does not appear to show any meaningful and/or consistent pattern. While they find that genic H3K9me2 is positively associated with the magnitude and extent of TE-derived H3K9me2, only the magnitude effect depends on gene distance. They find no relationship between TE epigenetic effects and expression of adjacent genes and some evidence of asymmetric H3K9me2 spreading due to neighboring highly expressed genes. Across all species, the magnitude of H3K9me2, but not the extent, was negatively correlated with TE population frequency. Similarly, the magnitude of H3K9me2, but not the extent, was correlated with TE abundance across species. Finally, they find that the magnitude of H3K9me2 (extent not tested) is positively correlated with expression of Su(var) genes across species.

The analyses performed here are rigorous and the appropriate replicates have been included. The authors also clearly go to great lengths to try to address confounding factors that could explain the correlations they are seeing.

1. My main concern related to evidence and reproducibility is that the results shown for H3K9me2 magnitude and spreading are highly variable between species and strains. It is not clear to me how much of this variation is due to biology versus technical issues related to ChIP, however I am surprised at how low the correlations among replicates are for these measures (Figure 1, Supplementary Figure S4). I would encourage the authors to perform additional QC measures for their ChIP experiments and replace any low quality replicates so as to minimize the contribution of technical variation. Some QC suggestions: the K-MetStat Panel spike-in should allow for an assessment of antibody specificity for the K9me2 modification and the Irreproducibility Discovery Rate (IDR) framework can quantify reproducibility among replicates.

We appreciate the reviewer’s helpful comment and performed analyses following the reviewer’s suggestion. We used ChIP-seq reads mapping to SNAP-ChIP K-MetStat Panel to assay the specificity of the used antibody (Abcam 1220). Consistent with previous investigations by modEncode (Egelhofer et al. 2011) and others (Sentmanat and Elgin 2012), our analysis found that this antibody shows high specificity against H3K9me2. We included this result in p45, line 23-25 and Figure 1 – Supplementary Figure S5.

We performed the suggested IDR analysis and found no direct associations between IDR and the correlations between replicates for the index of TE-mediated epigenetic effects (the diagnostic IDR plots are now in Figure 1 – Supplementary Figure S7). For instance, D. melanogaster shows the lowest correlation for the magnitude of TE-mediated H3K9me2 enrichment between replicates, but has decent consistency between replicates by IDR analysis.

A plausible technical reason for this discrepancy is that IDR analysis relies on called peaks. Custom pipelines, which were designed to identify “sharp peaks,” are likely to miss the enrichment of repressive marks, which typically have broad and lower level of enrichment (reviewed in Park 2009 Nature Review Genetics). Accordingly, TE-mediated H3K9me2 enrichment in euchromatic genome could be missed in the IDR analysis. Indeed, a previous study found that only 12% of TE-induced H3K9me2 enrichment was detected by peak calling (Lee et al. 2020 PLoS Genetics). Another more important biological reason is that the spreading of repressive marks from constitutive heterochromatin has been observed to be a variegating phenotype and could differ between cells and individuals. Accordingly, the variability for the magnitude and extent of TE-mediated local enrichment of heterochromatic marks could have contributed to the low correlation between replicates for some samples. We included these discussions in p49, line 10-29.

Thanks to the reviewer’s comment, we noticed that D. simulans strain 2 has much fewer peaks showing low IDR between replicates. To be conservative, we thus tried excluding this sample for the key analysis for our study and still found negative associations between the magnitude of TE-mediated H3K9me2 and genomic TE abundance. We included this analysis in p29, line 2224; p29, line 28-29; p49, line 29-p50, line 3.

2.1. Related to the comment above: does the "extent" measurement used here differ from the "spreading" measurement used in Lee (2015) and Lee and Karpen (2017). If so, why was a new metric used?

(We split the reviewer’s comment #2 into two and answered them individually).

We are a bit unsure about the “spreading” measurement the reviewer is referring to. In Lee (2015), we did not estimate TE-mediated local enrichment of H3K9me3 and instead investigated the associations between genic H3K9me3 enrichment and the presence of adjacent TEs. In Lee and Karpen (2017), we used the same measurement as the current study – the extent and magnitude of TE-mediated local enrichment of H3K9me2.

We agreed with the reviewer that we did not include enough descriptions of the two indexes for the epigenetic effects of individual TE and this can cause significant confusion. To address this, we added text in p10, line 3-9.

2.2. Could that explain why spreading no longer seems to be an important contributor to TE epigenetic effects?

We agree with the reviewer that we did not clearly define “epigenetic effects of TEs” and, without that, the term may implicitly imply functional consequences (e.g., gene expression). We used the term “epigenetic effects” to describe TE-mediated local enrichment of H3K9me2. We added definition for this term when first mentioning this term in Introduction (p5, line 28) and Results (p10, line 3-5).

The spreading of H3K9me2 beyond TE boundaries could lead to local enrichment of H3K9me2. The strength of such an event could be described by the enrichment level (magnitude) and extent (extent) of the local enrichment of H3K9me2 outside TE sequence. We see that the word “spreading” might also be interpreted as the “extent” of the H3K9me2 enrichment. Accordingly, we replaced this word throughout the manuscript with a more accurate description (“local enrichment of H3K9me2”) to avoid confusion and misinterpretation (e.g., p1, line 23; p4, line 15; p5, line 1; p38, line 10).

According to Figure 1B and C, we observed that the presence of TEs not only lead to H3K9me2 enrichment at immediate TE-flanking regions, but also could extend for several kbs, suggesting that the “spreading” (used in reviewer’s context) of repressive marks from TEs is still an important “epigenetic effects.” We also observed that the extent of TE-mediated enrichment of H3K9me2 is positively correlated with genic H3K9me2 enrichment (Figure 3A). It is when we investigated the role of TE-mediated local enrichment of H3K9me2 that we did not find the predicted associations between the extent of such effect and genomic TE abundance, but did find significant associations between the magnitude of such effect and TE abundance. It is worth noting that Lee and Karpen 2017 made similar observations with two species (i.e., TE copy number is associated with the magnitude of epigenetic effects, but not with the extent), which we noted in p30, line 1-3. In addressing the reviewer #1’s comment 8, we also made more clear distinctions between the two estimates (magnitude and extent of TE-mediated H3K9me2 enrichment; see above).

3. In terms of clarity: if the authors demonstrate that their ChIP replicates are highly reproducible and the "extent" metric is not meaningfully different from the previously used "spreading" metric, then it is hard for me to think of technical or analytical reasons why the findings reported here differ substantially from those reported previously. The previous model was intuitive: an underappreciated aspect of the deleterious nature of TE insertions is that euchromatic TEs can initiate heterochromatin formation which spreads into nearby genes leading to harmful reductions in their expression. This study suggests that magnitude of H3K9me2 is the major contributor to the "epigenetic effect" of TEs and that there is no relationship between TE-induced H3K9me2 and neighboring gene expression. I would encourage the authors to use the discussion to acknowledge these contradictions and try to reconcile these results with those from the previous two papers. Is there reason to believe the previous studies were incorrect? Otherwise, it is very hard to know what to take away from this study.

We agreed with the reviewer that without clearly defining “epigenetic effects,” readers could get the impression that the current study contradicts previous findings. Hopefully, our revision according to the reviewer’s comment #2 above would clarify the meaning of “epigenetic effects of TEs.” Also, as noted above, similar to Lee and Karpen (2017), we used two indexes to describe the epigenetic effects of TEs, or TE-mediated local enrichment of H3K9me2 in the euchromatic genome. Both the previous study with a limited number of species and our study found that (1) TEs lead to a significant increase in the magnitude and extent of local enrichment of H3K9me2, and (2) TE copy number is negatively associated with the magnitude, but not extent, of the effect. Our findings do not contradict previous studies.

The generally limited associations between TE-mediated epigenetic effects and lowered gene expression, as mentioned by the reviewer, was also reported in other model species, including in D. melanogaster (Lee and Karpen 2017 eLife), Arabidopsis thaliana (Quadrana et al. 2016 eLife; Stuart et al. 2016 eLife), rice (Choi and Purugganan 2018 MBE), and mouse (Pezic et al. 2014 Genes & Development). This has been a perplexing finding in the field of TE evolution, which was extensively discussed in two recent reviews (Kelleher et al. 2020 Trends in Genetics; Choi & Lee 2020, PLoS Genetics; discussed in our Introduction p5, line 9-16). Our current study investigated this question in multiple species and used both within-genome and across-genome (i.e. between homologous alleles) comparisons, further revealing the complex relationship between TE-mediated epigenetic effects and neighboring gene expression. We hope that our new analyses addressing reviewers' comments (reviewer #1’s comment 3 and reviewer #3’s comment 4) further help corroborate our findings. We also included additional discussions regarding the limited associations between TE-mediated epigenetic effects and gene expression (p39, line 9-22) and provided possible reasons why we did not find the expected negative associations between them (p39, line 22 – p40, line 16).

It is worth noting that our analysis focuses on identifying the genome-wide average pattern; it is still plausible that the local enrichment of H3K9me2 induced by individual TE occasionally lowers nearby gene expression (e.g., genes with positive z-scores in Figure 3C). In fact, many early examples of TE-mediated enrichment of repressive marks were discovered by the phenotypic consequences of reduced neighboring gene expression (reviewed in (Choi and Lee 2020 PLoS Genetics)). We included this discussion in p18, line 14-21.

In this study, we tested the model proposed by Lee and Karpen (2017) (Figure 6). We were able to find support for most of the models except for the generally expected, but unsupported prediction that TEs’ epigenetic effect would lower neighboring gene expression. We proposed several possibilities for this missing link, which can be found in p40, line 18 – p41, line 2.

Other major comments:

4. Figure 3C: For the gene expression analysis, it would be helpful to see a scatterplot where each point is a gene and the axes are the log2 expression level of the gene in the strain with an adjacent TE versus the log2 expression of the same gene in the strain without an adjacent TE. The slope of the regression line would show whether there is any tendency for genes with adjacent TEs to have lower expression.

We performed the suggested analysis following the reviewer’s suggestion (p18, line 9-11 and Figure 3 – Supplementary Figure 9). In these plots, we have the log2 expression level of genes in the strain with a TE on the X-axis and in the strain without TE on the Y-axis. We plotted genes whose nearest TEs with and without epigenetic effects with different colors. The slope is close to 1 (i.e., no difference in gene expression with and without TE presence), and there is no difference in slope for genes whose nearest TEs with and without epigenetic effects. These observations are consistent with our z-score analysis, finding a limited influence of TE-mediated epigenetic effects on gene expression.

5. It is interesting that there is clear evidence supporting enrichment of H3K9me2 at TEs leading to genic H3K9me2 enrichment but not expression downregulation. It would be very interesting if the authors could extend this analysis to include all genes from the euchromatic chromosome arms showing H3K9me2 enrichment in at least one strain or species, irrespective of TEs. Overall, are differences in genic H3K9me2 enrichment between strains/species associated with differences in gene expression?

This is a very interesting question raised by the reviewer. Yet, addressing this question is beyond the scope of this study. The suggested analysis would address the general role of the enrichment of H3K9me2 in gene expression, which is a broad topic.

It may be worth noting that, according to a previous analysis (Lee et al. 2020 PLoS Genetics), most, if not all, of the polymorphic H3K9me2 enrichment in the euchromatic genome is due to the presence/absence of TEs within species. Few H3K9me2 enriched regions in the euchromatic genome will thus be TE-free and could be suitable for performing the suggested analysis. On the other hand, the between-species comparison of gene expression is likely to be confounded by factors other than the enrichment of H3K9me2 (e.g., the evolution of regulatory sequences) and thus be challenging to interpret. We thus did not proceed with the suggested analysis.

6. I'm confused by the difference between Figure 3D and Figure 3 – Supplementary Figure S7. If I understand correctly, for D. melanogaster, there is no difference in the extent of spreading between the two sides of a TE when the highly expressed gene is near the TE, nor is there a difference when the highly expressed gene is far from the TE (Figure 3D). However, when expression measurements are used from a different strain (Figure 3 – S7), there is a clear difference in the extent of spreading on the gene side, specifically when the highly expressed gene is close to the TE. What explains these contradictory results? Are a large proportion of the genes identified as highly expressed when using the other strain's data no longer highly expressed in the strain of interest? (presumably due to the effect of the nearby TE?). If so, that seems notable as it lends support to an association between spreading and gene downregulation.

We appreciate the reviewer for pointing out this incongruence. Following the reviewer’s suggestion, we counted the number of genes that have a TE nearby and are categorized as highly expressed-

In both our data and modEncode: 179

Only in the modEncode data: 93

Only in our data: 83

There is no significant difference in the number of genes that are categorized as highly expressed in only modEncode data or only in our data. Accordingly, the suggested possibility (the difference between Figure 3D and previous Figure 3 – Supplementary Figure S7 is due to TE’s epigenetic effects lowering gene expression) may not be applicable to our observations.

The main reason why the result in Figure 3D differs from those of previous Figure 3 –

Supplementary Figure S7 is probably due to the fact that we used different statistics. In Figure 3D, we estimated the z-score for the difference in the extent of TE-mediated H3K9me2 enrichment between the genic and intergenic sides, normalized by the standard deviation of the estimates. In previous Figure 3 – Supplementary Figure S7, we simply compared the extent of TE-mediated H3Kme2 enrichment on the two sides of a TE. The z-score should be a more appropriate statistic given that it is normalized, so comparable across TE insertions. By using the z-score to perform the analysis in previous Figure 3 – Supplementary Figure S7, we found the same results like those in Figure 3D, and there is no statistically significant difference in the extent of TE-mediated H3K9me2 enrichment between the two sides of a TE in D. melanogaster.

Author response image 1

Given that there is no major difference between Figure 3D and updated Figure 3 – Supplementary Figure S7, and the latter does not add additional information, we decided to exclude analysis based on categorizing genes according to modEncode data (i.e., previous Figure 3 – Supplementary Figure S7).

7. The authors report that H3K9me2 spreading from TEs is counteracted by the presence of a highly expressed gene nearby (either upstream or downstream), leading to asymmetric spreading. However, the spreading extent index is calculated by averaging the extent of spreading between the two flanking regions of the TE. Is this the best choice given the presence of asymmetric spreading? Would using the longer of the two distances be better?

We think the reviewer is suggesting we use the “intergenic side,” which is less likely influenced by the expression of adjacent genes. Our response is based on this interpretation.

According to our data presented in Figure 3D, the extent of TE-mediated H3K9me2 enrichment is not symmetrical even for TEs far from genes. This is consistent with findings from constitutive heterochromatin that the spreading of repressive marks is a stochastic process. Also, our goal is to understand the potential functional and thus fitness consequences of TE-mediated H3K9me2 enrichment; such effect on the genic side may thus be of higher importance. Accordingly, it is challenging to predict a priori which side of the extent of H3K9m2 enrichment is more important, and we still chose to use the average of the two sides for our analysis.

8. An important prediction related to the proposed relationship between differences among species in expression of Su(Var)s and differences in the magnitude of H3K9me2 at TEs is that the same pattern should be present for other heterochromatic loci. Do other (non-TE) heterochromatic loci show stronger H3K9me2 enrichment in species with higher expression of Su(Var)s?

We agree with the reviewer that this is an interesting question. Yet, the suggested analysis and interpretation of the results are outside the scope of the current study. In addition, the effects of reduced dosage effects of Su(var)s were mainly observed at the boundaries of heterochromatin (e.g., at translocated genes located at the heterochromatin/euchromatin boundary and associated with PEV phenotype, reviewed in Elgin and Reuter 2013). In order to perform such analysis, one needs to first define the heterochromatin (HC)/euchromatin (EU) boundaries epigenetically using H3K9me2 enrichment. A potential issue is that the enrichment level of H3K9me2 at genes at HC/EU boundaries would highly depend on how the boundaries were first defined by H3K9me2, making the interpretation of the analysis challenging. Because of these reasons, we did not proceed with the suggested analysis.

9. It is not clear to me how the 1360 reporter experiment is applicable to the Su(Var)s expression results. The authors only consider H3K9me2 magnitude in their comparison to Su(Var)s expression rank yet, if I understand correctly, the 1360 reporter experiment measures spreading/extent of H3K9me2 which leads to downregulation of the adjacent mini-white gene. The reporter experiment described here suggests that reducing Su(Var) gene dosage leads to a reduction of H3K9me2 spreading along with a concomitant reduction in reporter gene silencing, however the other results reported here suggest that spreading of H3K9me2 is not an important contributor to the epigenetic effects of TE nor does spreading lead to gene downregulation.

We agree with the reviewer that our motivation for performing the 1360 reporter experiment was unclear and needs to be further clarified. As discussed in response to the reviewer’s comments 2 & 3, we did find that euchromatic TEs lead to local enrichment of H3K9me2, which is characterized by higher H3K9me2 enrichment level (magnitude) and an extent beyond the TE sequence. Even though on a genome-wide average, we did not find predominantly negative impacts of TE’s epigenetic effects on neighboring gene expression, our analysis does not preclude individual cases where a gene’s expression is reduced by TE-mediated enrichment of H3K9me2 (see discussions in comment 3 above and p18, line 14-21).

The 1360 reporter consists of a TE right next to the reporter gene. Accordingly, for this reporter system, the magnitude, instead of the extent, of TE-mediated H3K9me2 enrichment may be more important. In this specific reporter system, there are associations between the presence of 1360, the enrichment level of H3K9me2 at the mini-white reporter gene, and the reduced amounts of red-eye pigmentation (likely due to reduced gene expression). Accordingly, the eye pigmentation level in this system serves as a convenient readout for the magnitude of TE-mediated H3K9me2 enrichment. We clarified these points in p33, line 23-29.

Minor Comments

1. Lee (2015) examined H3K9me3 while Lee and Karpen (2017) and this manuscript study H3K9me2 levels. I think it is important to explain why one modification was used versus the other. Do you expect them to be redundant or might they play different roles in spreading versus magnitude of heterochromatin?

We agree with the reviewer that we did not have enough context for why we chose to study H3K9me2. Both H3K9me2 and H3K9me3 are found to be highly enriched in the heterochromatic regions of the genome in D. melanogaster (Riddle et al. 2011 Genome Research; Kharchenko et al. 2011 Nature). According to Lee and Karpen (2017), both H3K9me2 and H3K9me3 captured TE-mediated local enrichment of repressive epigenetic marks. The major reason that our study focused on H3K9me2 is due to the availability of suitable antibodies. Antibody against H3K9me2 is monoclonal and has been validated by modEncode to show high consistency (Egelhofer et al. 2011), while commonly used H3K9me3 ChIP antibody available at the time when this study started (Abcam ab8898) is polyclonal and likely cross-react with other histone modifications, in particular H3K27me3 (personal communications with modEncode team). Because H3K9me2 and H3K9me3 would both provide the needed information about TE-mediated local enrichment of heterochromatic marks, we proceeded with only H3K9me2. To provide more context, we added text at p7, line 6-12.

2. For the TE population frequency analysis, I was not able to find the cutoffs for what were considered high frequency and low frequency TEs

This is an error at our end; we forgot to include details about the categorization of TEs based on population frequencies. High-frequency TEs are those that are found in both the studied strain and the population sample, while low-frequency TEs are those that are unique to the studied strain and absent in the population sample. We included this information p23, line 19-21.

3. The Su(Var) section is a bit confusing to interpret because gene rank should correlate negatively with H3K9me2 but expression level should correlate positively and thus, both positive and negative correlations are referenced in this section when describing the same pattern.

Could you switch the gene rank order (i.e. lowest to highest) so that the expected correlation is in the same direction regardless of whether you are referring to expression level or expression rank?

In order to be consistent, we used expression rank for all expression analyses presented in Figure 3 (TE’s impact on gene expression) and Figure 5 (Su(var) expression). To help readers interpret the results, we added additional texts to either define expression rank (“rank from the highest expressed genes”) or interpret the result (“lower expression rank (i.e., higher expression)”, e.g., p16, line 7; p17, line 26). We hope these could help readers and avoid confusion.

4. Is there a relationship between local recombination rate and H3K9me2 extent and/or magnitude? If so, that could also be a potential confounder with respect to TE population frequency. I believe this was addressed in Lee and Karpen (2017) but it would be useful to reiterate here.

We repeated a similar analysis as Lee and Karpen (2017) in D. melanogaster, the only species in which a high-resolution recombination map is publicly available. Consistently, we found no associations between TE-mediated epigenetic effects and local recombination rate, further corroborating the previous analysis. We included these results in p24, line 17-29.

Reviewer #3 (Significance (Required)):

This study builds upon two previous studies that have made a large impact on the field. This work is significant because it suggests that the epigenetic effects of euchromatic TE insertions are more complicated than previously suggested by work that was only conducted in D.

melanogaster. It also suggests that the mechanisms by which TEs exert these effects and the mechanisms that make these effects harmful need to be revised and studied further.

This work will be of general interest to those whose work involved evolutionary biology, population genetics, epigenetics and chromatin, and transposable elements.

My expertise is in the evolutionary and functional genomics of Drosophila.lex and 3 in yakuba species comple

https://doi.org/10.7554/eLife.81567.sa2

Article and author information

Author details

  1. Yuheng Huang

    Department of Ecology and Evolutionary Biology, University of California, Irvine, Irvine, United States
    Contribution
    Data curation, Formal analysis, Validation, Investigation, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1586-953X
  2. Harsh Shukla

    Department of Ecology and Evolutionary Biology, University of California, Irvine, Irvine, United States
    Contribution
    Software, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Yuh Chwen G Lee

    Department of Ecology and Evolutionary Biology, University of California, Irvine, Irvine, United States
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    grylee@uci.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0081-7892

Funding

National Institutes of Health (R00GM121868)

  • Yuh Chwen G Lee

National Institutes of Health (R35GM14292)

  • Yuh Chwen G Lee

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We greatly appreciate JJ Emerson and Peter Andolfatto for generously providing the PacBio genomes and corresponding Drosophila strains, Sally Elgin for providing 1360-mini-white strains, Giacomo Cavalli for sharing ORw1118 line, Bloomington Drosophila Stock Center for providing Su(var) mutants, and David Acevedo and Jasmine Osei-Enin for technical assistance. Patrick Reilly and Mahul Chakraborty provided helpful guidance on the PacBio genome assembly, and Kevin Brick for discussion about ChIP-seq normalization. We thank University of California High-Throughput Genomics Facility and High Performance Cluster at UC Irvine for sequencing and computational resources. We also appreciate Aneil Agrawal, Ching-Ho Chang, Jae Choi, Brandon Gaut, Mia Levine, Aline Muyle, and members of the Lee lab for providing helpful discussions and comments on the manuscript and Leila Lin for assistance in drawing the model figure. We also thank three reviewers for their helpful comments. This work was supported by NIH R00GM121868 and R35GM14292 to YCGL.

Senior Editor

  1. Molly Przeworski, Columbia University, United States

Reviewing Editor

  1. Magnus Nordborg, Gregor Mendel Institute, Austria

Publication history

  1. Preprint posted: March 12, 2022 (view preprint)
  2. Received: July 4, 2022
  3. Accepted: July 15, 2022
  4. Version of Record published: August 23, 2022 (version 1)

Copyright

© 2022, Huang et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,619
    Page views
  • 483
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Yuheng Huang
  2. Harsh Shukla
  3. Yuh Chwen G Lee
(2022)
Species-specific chromatin landscape determines how transposable elements shape genome evolution
eLife 11:e81567.
https://doi.org/10.7554/eLife.81567

Further reading

    1. Epidemiology and Global Health
    2. Evolutionary Biology
    Erin Brintnell, Art Poon
    Insight

    Combining clinical and genetic data can improve the effectiveness of virus tracking with the aim of reducing the number of HIV cases by 2030.

    1. Evolutionary Biology
    Julian M Wagner, C Jaco Klok ... Jon F Harrison
    Research Article Updated

    The scaling of respiratory structures has been hypothesized to be a major driving factor in the evolution of many aspects of animal physiology. Here, we provide the first assessment of the scaling of the spiracles in insects using 10 scarab beetle species differing 180× in mass, including some of the most massive extant insect species. Using X-ray microtomography, we measured the cross-sectional area and depth of all eight spiracles, enabling the calculation of their diffusive and advective capacities. Each of these metrics scaled with geometric isometry. Because diffusive capacities scale with lower slopes than metabolic rates, the largest beetles measured require 10-fold higher PO2 gradients across the spiracles to sustain metabolism by diffusion compared to the smallest species. Large beetles can exchange sufficient oxygen for resting metabolism by diffusion across the spiracles, but not during flight. In contrast, spiracular advective capacities scale similarly or more steeply than metabolic rates, so spiracular advective capacities should match or exceed respiratory demands in the largest beetles. These data illustrate a general principle of gas exchange: scaling of respiratory transport structures with geometric isometry diminishes the potential for diffusive gas exchange but enhances advective capacities; combining such structural scaling with muscle-driven ventilation allows larger animals to achieve high metabolic rates when active.