
==== Front
Mol Biol Evol
Mol Biol Evol
molbev
Molecular Biology and Evolution
0737-4038
1537-1719
Oxford University Press UK

39189646
10.1093/molbev/msae181
msae181
Discoveries
AcademicSubjects/SCI01130
AcademicSubjects/SCI01180
Prevalent Fast Evolution of Genes Involved in Heterochromatin Functions
Lin Leila Department of Ecology and Evolutionary Biology, University of California, Irvine, CA, USA

https://orcid.org/0000-0002-1586-953X
Huang Yuheng Department of Ecology and Evolutionary Biology, University of California, Irvine, CA, USA

McIntyre Jennifer Department of Ecology and Evolutionary Biology, University of California, Irvine, CA, USA

https://orcid.org/0000-0001-9361-1190
Chang Ching-Ho Division of Basic Sciences, Fred Hutchinson Cancer Center, Seattle, WA, USA

Colmenares Serafin Department of Cell and Molecular Biology, University of California, Berkeley, CA, USA

https://orcid.org/0000-0002-0081-7892
Lee Yuh Chwen G Department of Ecology and Evolutionary Biology, University of California, Irvine, CA, USA

Mulligan Connie Associate Editor
*Corresponding author: E-mail: grylee@uci.edu.
9 2024
27 8 2024
27 8 2024
41 9 msae18115 3 2024
14 8 2024
20 8 2024
18 9 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com.

Abstract

Heterochromatin is a gene-poor and repeat-rich genomic compartment universally found in eukaryotes. Despite its low transcriptional activity, heterochromatin plays important roles in maintaining genome stability, organizing chromosomes, and suppressing transposable elements. Given the importance of these functions, it is expected that genes involved in heterochromatin regulation would be highly conserved. Yet, a handful of these genes were found to evolve rapidly. To investigate whether these previous findings are anecdotal or general to genes modulating heterochromatin, we compile an exhaustive list of 106 candidate genes involved in heterochromatin functions and investigate their evolution over short and long evolutionary time scales in Drosophila. Our analyses find that these genes exhibit significantly more frequent evolutionary changes, both in the forms of amino acid substitutions and gene copy number change, when compared to genes involved in Polycomb-based repressive chromatin. While positive selection drives amino acid changes within both structured domains with diverse functions and intrinsically disordered regions, purifying selection may have maintained the proportions of intrinsically disordered regions of these proteins. Together with the observed negative associations between the evolutionary rate of these genes and the genomic abundance of transposable elements, we propose an evolutionary model where the fast evolution of genes involved in heterochromatin functions is an inevitable outcome of the unique functional roles of heterochromatin, while the rapid evolution of transposable elements may be an effect rather than cause. Our study provides an important global view of the evolution of genes involved in this critical cellular domain and provides insights into the factors driving the distinctive evolution of heterochromatin.

heterochromatin
protein evolution
transposable elements
repetitive sequences
arms race
positive selection
Cancer Research Foundation 10.13039/100002002 NIH 10.13039/501100012264 NIH R35GM139653 NIH R35-GM14292
==== Body
pmcIntroduction

Heterochromatin, first discovered as the darkly stained chromosomal regions that remain condensed throughout the cell cycle (Heitz 1928), is a distinct cytological domain that is conserved across eukaryotic cells (Liu et al. 2020). While a small fraction of heterochromatin was found to be cell-type specific, or “facultative,” the majority of these chromosomal blocks remain condensed across different cell types and are known as “constitutive” heterochromatin (Heitz 1928). Constitutive heterochromatin (referred to as “heterochromatin” for simplicity hereafter) is usually located around centromeres and telomeres (Janssen et al. 2018), and its underlying DNA is depleted of functional genes. Instead, heterochromatin is mainly composed of repetitive sequences, including satellite repeats (Peacock et al. 1978) and transposable elements (TEs) (Hoskins et al. 2007, 2015). Accordingly, heterochromatin is oftentimes assumed to be functionally inert and nicknamed the “dark matter” of the genome. Yet, studies have found heterochromatin playing critical roles in many chromosomal functions, such as maintaining genome stability, mediating chromosome segregation, and ensuring proper DNA repair of repeat-rich sequences, across eukaryotes (reviewed in Feng and Michaels 2015; Allshire and Madhani 2018; Janssen et al. 2018; Kendek et al. 2021). Not surprisingly, disruption of heterochromatin functions has been linked to various diseases (Hahn et al. 2010), aging progression (Villeponteau 1997; Lee, Kim, et al. 2020), and cancer (Janssen et al. 2018).

The critical functions and interspecies conservation of heterochromatin domains would naturally lead to the expectation that genes modulating the function of heterochromatin should be highly conserved. Surprisingly, however, a handful of these genes show nonneutral evolution of their amino acid sequence and rapid turnover of gene copy number, both of which are consistent with a history of positive selection. The evolution of some of these genes was of interest due to them being core structural components of heterochromatin, such as the heterochromatin protein 1 (HP1) family (Levine et al. 2012; Helleu and Levine 2018) and its paralogs (Vermaak et al. 2005; Ross et al. 2013). The rest of the genes were observed to have roles in the evolution of small RNA pathways (Obbard et al. 2006, 2009), transcription factors (Kasinathan et al. 2020), or germline stem cells (Flores et al. 2015), as well as speciation (Barbash et al. 2003). Because satellite repeats and TEs are enriched in heterochromatic sequence and can defy Mendel's law to increase their transmission to the next generation (Cosby et al. 2019; Henikoff et al. 2001), an arms race between heterochromatic sequence and host genes suppressing such selfish behaviors has been a common theme explaining the rapid evolution of these genes (e.g. Vermaak et al. 2005; Satyaki et al. 2014; Helleu and Levine 2018; Brand and Levine 2021). Further supporting such conjecture, the composition, abundance, and location of repetitive sequences in heterochromatin are found to diverge rapidly even between closely related species (Wei et al. 2014, 2018; Kim et al. 2021; de Lima and Ruiz-Ruano 2022).

However, this handful of genes represents only a small fraction of genes involved in modulating heterochromatin functions. The unique and essential functions of heterochromatin depend on its specific enrichment of repressive di- and tri-methylation at histone H3 Lysine 9 (H3K9me2/3; Allshire and Madhani 2018) and the “reader” of these marks, HP1a, which serves as the foundational structural component of heterochromatin (Eissenberg and Elgin 2014) and leads to DNA compaction (Verschure et al. 2005) as well as transcriptional suppression (Li et al. 2003). Interestingly, the binding of HP1a to H3K9me2/3 further recruits “writers” of H3K9me (histone methyltransferase; Schotta 2002), and this positive feedback mechanism between readers and writers leads to the propagation of such repressive epigenetic marks independent of the underlying DNA sequence (Allshire and Madhani 2018). This unique ability of heterochromatin to “spread” is best exemplified by the position effect variegation (PEV) first discovered in Drosophila (Muller 1930)—the mosaic expression of euchromatic genes translocated to heterochromatin-proximal regions due to the stochastic spreading of repressive H3K9me2/3 from heterochromatin (Girton and Johansen 2008; Elgin and Reuter 2013). By studying mutants that either enhance or weaken PEV, many genes involved in heterochromatin functions were identified. These include not only the structural (e.g. Shaffer et al. 2006) and enzymatic (e.g. Czermin et al. 2001) components for heterochromatin but also those that antagonize the initiation and/or maintenance of heterochromatin (e.g. writers of antagonizing histone modifications; Bao et al. 2007). Still, other genes were discovered through their colocalization with heterochromatin domain either cytologically (e.g. Swenson et al. 2016) or epigenomically (e.g. Alekseyenko et al. 2014) and were subsequently identified to be involved in heterochromatin function.

Given the large numbers and vastly diverse functional roles of genes involved in heterochromatin function, we conducted an expansive survey to determine whether the previously reported rapid evolution of a small subset of these genes is anecdotal or a common feature of proteins involved in the functions of this unique chromatin environment. To do so, we compiled an exhaustive list of candidate genes that have been shown, or are likely, to be involved in heterochromatin functions, including PEV modifiers, histone-modifying enzymes influencing H3K9me2/3 enrichment, and genes whose protein products localize to heterochromatin. Our investigation finds that these genes evolve exceptionally fast, a pattern that is not general to genes interacting with other repressive chromatin marks and suggests a selective pressure unique to constitutive heterochromatin. We further dissect the domains and protein properties targeted by positive selection and specifically test the premise that the rapid evolution of genes involved in heterochromatin functions could be driven by fast-changing repetitive sequences. Based on our findings, we propose an evolutionary model where the rapid changes in these genes are the unavoidable consequence of the unique functional roles of heterochromatin, and that the evolution of repetitive sequences may be the consequence, instead of the cause. Our study provides an important global view of the evolution of genes involved in this critical cellular domain and sheds light on the drivers behind the unique evolution of heterochromatin.

Results

Identification of Candidate Genes Involved in Heterochromatin Functions

With the goal of identifying the global evolutionary patterns for genes involved in the functions of heterochromatin, we selected candidate genes based on three criteria to maximize inclusivity (supplementary table S1, Supplementary Material online). The first category contains genes demonstrated to modulate PEV (Girton and Johansen 2008; Elgin and Reuter 2013; Swenson et al. 2016). Genes whose mutations weaken PEV and, consequently, heterochromatin functions are known as suppressors of variegation (Su(var)), while those that augment PEV are Enhancers of Variegation (E(var)). Here, we included 59 and seven previously identified Su(var)s and E(var)s, respectively.

The second category focuses on histone-modifying enzymes that can influence the enrichment levels of H3K9me2/3. Enrichment of H3K9me2/3 requires “eraser” proteins that first remove antagonizing active marks, such as H3K9ac (deacetylases; Czermin et al. 2001) and H3K4me3 (demethylase; Rudolph et al. 2007). This is followed by writer proteins depositing H3K9me2/3 (H3K9 methyltransferases; Schotta 2002). We consider these histone-modifying enzymes to enhance H3K9me2/3 enrichment and, thus, heterochromatin functions. On the contrary, erasers for H3K9me2/3 (H3K9 demethylase; Herz et al. 2014), as well as writers for acetylation (acetyltransferase; Kuo and Andrews 2013) and S10 kinase (Deng et al. 2008), are known to antagonize the deposition and maintenance of H3K9me2/3. H4K20me3 is another conserved hallmark of heterochromatin (Schotta et al. 2004), and we included the corresponding methyltransferases. In total, our list contains 16 and seven histone-modifying enzymes that enhance or weaken heterochromatin functions, respectively.

The last category includes genes whose protein products colocalize with heterochromatin, suggesting they play roles in this subnuclear compartment. We surveyed the literature for cytology-based evidence of colocalization with HP1a through either immunofluorescence (e.g. Greil et al. 2007) or live imaging (e.g. Swenson et al. 2016) as well as epigenomics through chromatin immunoprecipitation (e.g. Alekseyenko et al. 2014; Kasinathan et al. 2020). Because a previous comprehensive survey found that some HP1a-binding proteins identified through immunoprecipitation-mass spectrometry are not necessarily enriched in heterochromatin cytologically (Swenson et al. 2016), genes implicated in HP1a-binding but lacking localization evidence are not included in our analysis. Hereafter, we referred to this largest category with 63 genes as “heterochromatin proteins.”

In total, we studied the evolution of 106 genes involved in heterochromatin functions, which we termed “heterochromatin-related genes” hereafter (see Fig. 1a for the number of genes in each category and supplementary table S1, Supplementary Material online for the full list of genes and references). It is worth noting that 43 candidate genes belong to more than one category because these three categorizations are not mutually exclusive (supplementary table S1, Supplementary Material online). We consider genes in categories of Su(var)s, histone-modifying enzymes enhancing H3K9me2/3 or H4K20me3 enrichment, or heterochromatin proteins as “mediating” heterochromatin functions, while those defined as E(var)s or histone-modifying enzymes weakening H3K9me2/3 enrichment as “antagonizing” heterochromatin functions. To determine whether the evolutionary histories of heterochromatin-related genes are exceptional, we compared them with Polycomb group genes (Kassis et al. 2017), whose protein products are enriched at facultative heterochromatin. Even though also generally associated with transcriptional suppression, facultative heterochromatin is restricted to developmentally regulated genes in euchromatin and is enriched for another type of repressive histone modification (H3K27me3) (Bannister and Kouzarides 2011; Bell et al. 2023). Polycomb genes thus allowed us to determine whether the exceptional evolutionary histories, if any, of heterochromatin-related genes are due to their involvement in a repressive chromatin environment or specifically for heterochromatin structure.

Fig. 1. Positive selection on protein sequences of heterochromatin-related genes found by evolutionary tests at both short and long evolutionary time scales. a) Barplot showing the number of genes in each category of heterochromatin-related genes and the Polycomb control genes. Those expected to mediate (Su(var)s, +HC enzymes, and HC protein) or antagonize (E(var)s and -HC enzymes) heterochromatin functions are in hues of green and orange, respectively. −HC enzyme: histone-modifying enzymes weakening H3K9me2/3 enrichment; +HC enzyme: histone-modifying enzymes enhancing H3K9me2/3 enrichment. b) Phylogenetic tree for the species included in the study. Note that the branch lengths are not to scale. c, e) Violin plots for α, the proportion of amino acid substitutions fixed by positive selection inferred by the MK test c) and dN/dS ratio, the relative rates of nonsynonymous substitutions e), for different categories of genes. Numbers of significant genes are in parentheses. d, f) Barplots showing the proportion of genes with significant MK tests d) and accelerated rates of dN/dS driven by positive selection f). The number of genes in each category is either on top of the barplot a) or in parenthesis d, f). Mann–Whitney test c, e) and binomial test d, f): *P < 0.05, **P < 0.01, and ***P < 0.001.

Heterochromatin-related Genes Experience Pervasive Positive Selection on Protein Sequences across Short- and Long-Time Scales

We performed evolutionary genetic tests to identify heterochromatin-related genes with signatures of positive selection at two time scales—contrasting the polymorphism within Drosophila melanogaster and divergence between D. melanogaster and Drosophila simulans (∼4 million years divergence between lineages) under the framework of unpolarized McDonald–Kreitman (MK) tests (McDonald and Kreitman 1991) and phylogenetic analysis by maximum likelihood (PAML) tests (Yang 2007) using 16 Drosophila species with an estimated divergence time of ∼25 million years (Suvorov et al. 2022) (see Fig. 1b for studied species). In addition to comparing the evolution of heterochromatin-related genes to that of Polycomb genes, we made the contrast to the evolution of a thousand randomly sampled genes to represent genome-wide estimates (referred to as “random genes” hereafter).

On a short evolutionary time scale, we first estimated the proportion of amino acid substitutions that might have been driven by positive selection (α; Smith and Eyre-Walker 2002). Heterochromatin-related genes, as a group, have significantly larger α than those of Polycomb (Mann–Whitney U tests, P = 0.0244, median = 0.0847 [heterochromatin-related genes] and −0.424 [Polycomb]) and random genes (Mann–Whitney U tests, P < 10−8, median = −0.516 [random genes]). Breaking down our candidate genes according to functions revealed that most categories have larger α than Polycomb control or random genes, even though the comparisons are statistically significant only for Su(var) and heterochromatin protein (Fig. 1c). These observations suggest the possibility that heterochromatin-related genes experienced frequent positive selection. Consistently, we identified that 21 out of 104 (20.2%) heterochromatin-related genes with sufficient polymorphism data to perform MK tests have evidence of adaptive evolution (rejection of the null hypothesis and an excess of nonsynonymous fixed differences between species). This is a significantly larger proportion than that of the Polycomb genes (5.56%; Fisher's exact test, P = 0.191; binomial test, P < 10−6, odds ratio = 4.30) and random genes (6.94%, Fisher's exact test, P < 10−3; binomial test, P < 10−5, odds ratio = 2.90), while we found no difference in this proportion between Polycomb genes and random genes (Fisher's exact test and binomial test, P = 1). Specifically, Su(var) (22.0%) and heterochromatin proteins (26.2%), both of which are expected to mediate heterochromatin functions, have much larger proportions of genes with evidence of positive selection than either Polycomb or random genes (Fig. 1d). Indeed, candidate genes mediating heterochromatin functions, as a group, show significantly more evidence of adaptive evolution than the Polycomb (20.2%, Fisher's exact test, P = 0.190; binomial test, P < 10−6, odds ratio = 4.30) as well as random genes (Fisher's exact test, P < 10−3; binomial test, P < 10−4, odds ratio = 2.91). It is worth noting that the number of codons involved in the MK tests is similar across categories of genes (Mann–Whitney U tests, P > 0.05 for all comparisons), suggesting that these observations are unlikely driven by differences in the statistical power of the MK tests.

On a longer evolutionary time scale, we estimated the relative rates of amino acid substitution (nonsynonymous divergence/synonymous divergence [dN/dS]) across 16 studied Drosophila species. There is an overall trend that heterochromatin-related genes have larger dN/dS ratios than the Polycomb control genes (medians = 0.0974 [heterochromatin-related genes] vs. 0.0852 [Polycomb controls]) and random genes (medians = 0.0802), and this difference is statistically significant for comparisons to random genes (Mann–Whitney U test, P = 0.167 [compared to Polycomb control] and 0.0032 [compared to random genes]). Nevertheless, only the category of “heterochromatin protein,” the largest category, has a statistically larger dN/dS ratio than that of the Polycomb genes (Mann–Whitney U test, P = 0.0204) and random genes (Mann–Whitney U test, P < 10−5; Fig. 1e). To test whether the changes in dN/dS ratios of heterochromatin-related genes over the phylogenetic tree might have been driven by positive selection, we compared the log-likelihood for two alternative models estimating dN/dS at each site: the null model that assumes dN/dS ratio is beta-distributed and not greater than 1 (M8a) and the alternative model in which dN/dS ratio is allowed to be greater than 1 (M8 models; Yang 2007). A significantly larger likelihood of the alternative model is consistent with a history of positive selection acting on the gene (Swanson et al. 2003). We found that the substitution patterns of 42 heterochromatin-related genes (39.6%) better fit the alternative model, suggesting frequent positive selection acting on them. This proportion is significantly larger than that of Polycomb genes (26.3%; Fisher's exact test, P = 0.314; binomial test, P = 0.00276, odds ratio = 1.84) and the random genes (24.8%; Fisher's exact test, P = 0.00159; binomial test, P = 0.000975, odds ratio = 1.99), while there was no difference in this proportion between the Polycomb and random genes (Fisher's exact test, P = 0.795; binomial test, P = 0.796). Similar to observations made on a short evolutionary time scale, Su(var)s and heterochromatin proteins have stronger evidence of positive selection compared to Polycomb and random genes (Fig. 1f), as do all genes mediating heterochromatin functions (41.4%, Fisher's exact test, P = 0.305; binomial test, P = 0.00125, odds ratio = 1.97 [compared to Polycomb control] and Fisher's exact test, P = 0.0115; binomial test, P = 0.000272, odds ratio = 1.67 [compared to random genes]). Interestingly, while histone-modifying enzymes that enhance heterochromatin, as a group, do not exhibit exceptional rates of protein evolution, likely due to the low statistical power associated with its small number of genes, all methyltransferases for H3K9me2/3 (Su(var)3-9, egg, and G9a) show accelerated rates of protein evolution (supplementary table S1, Supplementary Material online). Overall, we found that heterochromatin-related genes, especially those mediating heterochromatin functions, experienced frequent positive selection on their protein sequence at both short- and long-time scales, even when compared to Polycomb genes. This observation suggests that the widespread adaptive evolution of heterochromatin-related genes is likely driven by selection acting on their functions specific to the heterochromatin environment, rather than by general selective pressure acting on genes interacting with repressive epigenetic marks.

Rapid Evolution of Heterochromatin-related Genes Also Manifests as Changes in Gene Copy Number

In addition to changes in amino acid sequences, the rapid evolution of genes can be in the form of changes in gene copy number (Hastings et al. 2009), and dramatic turnover of gene copy number for a small set of heterochromatin-related genes was previously reported (Levine et al. 2012; Lewis et al. 2016; Lee et al. 2017; Helleu and Levine 2018). Changes in gene copy number may hold particular significance for heterochromatin functions—several enzymatic and structural proteins of heterochromatin have been shown to exhibit dosage-dependent effects (Elgin and Reuter 2013). Accordingly, changes in gene copy number of heterochromatin-related genes could have immediate functional consequences.

We performed reciprocal BLAST searches using D. melanogaster amino acid sequences as queries to identify homologs and paralogs in other species, followed by manual curations (see Materials and Methods). With these, we found 18 heterochromatin-related genes having differences in gene copy number among 16 Drosophila species studied (17.0% of the heterochromatin-related genes), with 17 genes having gains of copies and one gene loss when compared to D. melanogaster. Similar to our analyses for the evolution of amino acid sequences (Fig. 1), the proportion of heterochromatin-related genes with copy number variation (CNV) is significantly more than that of Polycomb genes (5.26%, Fisher's exact test, P = 0.302, binomial test, P < 10−4, odds ratio = 3.68). Su(var) genes and heterochromatin proteins again exhibit exceptionally large proportions of genes with CNV (Fig. 2a). CNV of several of these genes has been studied before (e.g. HP1 [Levine et al. 2012], AGO2 [Lewis et al. 2016], Cav and Nap-1 [Lee et al. 2017], and mh [Brand et al. 2024]) and our following discussions mainly focused on those first identified by our study.

Fig. 2. Rapid evolution of heterochromatin-related genes manifested in multiple forms. a) Barplot showing the proportion of heterochromatin-related genes and Polycomb control genes with CNV. The numbers of genes with CNV are in parentheses. b) Summary of the duplication and loss events as well as sex-biased expression for heterochromatin-related genes identified to have CNV. Each column represents one orthologous copy. The number of duplication events is noted at the bottom. The colors of each square represent sex-biased expression, except for the 19 copies of mod duplicates in D. eugracilis, which is shown at the top left of b). c) Venn diagram depicting the number of genes identified to have exceptional evolution with one of the three evolutionary tests as well as their overlaps, with those identified by more than one test listed. a.a., amino acid. Binomial tests: **P < 0.01; ***P < 0.001.

By using synteny to infer the orthologous relationships between duplicates in different species (see Materials and Methods), we found that most observed CNVs for heterochromatin-related genes involve a single or two duplication events (Fig. 2b; supplementary fig. S1 and table S2, Supplementary Material online). The most notable exception is the nucleolin homolog, mod, which has 18 duplications in Drosophila eugracilis alone. mod plays important roles in morphogenesis (Graba et al. 1994) and spermatogenesis (Park et al. 2023) and has been shown as a dosage-dependent Su(var) (Garzino et al. 1992), suggesting these identified CNVs could readily alter heterochromatin functions. In addition, we observed complex duplication/loss events leading to the CNV observed in mof (supplementary fig. S2, Supplementary Material online), an acetyltransferase that influences H3K9me2/3 enrichment (Feller et al. 2015) and is known to be involved in dosage compensation (Hilfiker et al. 1997). In addition to the duplication happening on the lineage leading to Drosophila kikkawai, there is likely a duplication event of mof in the common ancestor of the oriental lineage (see Fig. 1b), followed by a subsequent loss in the subsets of lineages that leads to the paraphyletic presence of a mof duplicate (supplementary fig. S2, Supplementary Material online).

Duplicated genes are commonly found to have male-biased expression, which is suggested as a resolution to sexual genetic conflicts (Gallach and Betrán 2011). We are interested in examining whether duplicates of heterochromatin-related genes exhibit similar trends of sex-biased expression. By using publicly available transcriptome data for the species studied, we categorized the parental ortholog and duplicated paralog (with respect to D. melanogaster) as male-biased, female-biased, or unbiased in expression (see Materials and Methods). With the exception of MTA-like and Su(var)2-10, most identified duplicates indeed show male-biased expression, irrespective of whether the parental copy is male or female biased in expression (Fig. 2b). Interestingly, for mod in D. eugracilis, the parental copy has a female-biased expression, while the 18 duplicates exhibit a mixture of male-biased, female-biased, and unbiased expression. Another notable case is del, which is part of a germline complex enabling the transcriptions of TE-targeting small RNAs from loci in pericentromeric heterochromatin (Mohn et al. 2014) and shows female-biased expression. While the del duplicate in D. simulans shows male-biased expression, the orthologous duplicated copy in the sister species Drosophila sechellia has no biases in expression. In summary, the rapid evolution of heterochromatin-related genes, when compared to Polycomb genes, is also reflected in their changes in gene copy number.

Heterochromatin-related Genes with Diverse Functions are Recurrent Targets of Positive Selection

The three evolutionary tests we conducted detect signals of positive selection at different time scales and in different forms (amino acid sequences vs. gene copy number). We were interested in investigating whether some heterochromatin-related genes may exhibit rapid evolution with multiple tests. Consistent with analyses conducted separately for different evolutionary tests, Su(var)s and heterochromatin proteins have an especially large proportion of genes with evidence of rapid evolution in any one of the tests (binomial test, P < 0.001 for both comparisons; supplementary fig. S3, Supplementary Material online), and this difference becomes even more prominent when considering genes with exceptional evolution for more than one test (binomial test, P < 0.001 for both comparisons; supplementary fig. S3, Supplementary Material online). Three genes show evidence of rapid evolution for all three tests (Fig. 2c): AGO2, the core effector in endogenous siRNA silencing pathway that silences TEs in somatic cells (Ghildiyal et al. 2008), Su(var)2-10, which critically links piRNA targeting and the transcriptional silencing of TEs (Ninova et al. 2020), and mh, whose evolution is the focus of a recent study (Brand et al. 2024). Other heterochromatin-related genes with more than one evidence of rapid evolution (Fig. 2c) include those related to silencing of TEs (e.g. del, moon, sov, Trf2, and wde), dosage compensation (e.g. Su(var)3-7), maintenance of chromatin structure (e.g. Nipped-B), and structural components of heterochromatin (e.g. HP5). While most of these genes are either Su(var) or heterochromatin proteins, several histone-modifying enzymes were also found to undergo positive selection with multiple evolutionary tests, such as the H3K4me3 demethylase NO66 and acetyltransferase mof. These observations suggest that positive selection recurrently acts on heterochromatin-related genes with diverse functions.

Positive Selection Targets Both Ordered Domains and Intrinsically Disordered Regions of Proteins Encoded by Heterochromatin-Related Genes

To identify the potential evolutionary drivers for the observed rapid evolution of heterochromatin-related genes, we investigated the locations of positive selection within windows of each rapidly evolving gene. We identified windows with evidence of adaptive evolution by performing sliding MK tests and located sites with high Bayes Empirical Bayes (BEB) posterior probability of positive selection under the PAML framework (Yang et al. 2005). Locations of these windows or sites with evidence of positive selection were then contrasted with the annotated and/or predicted domains of heterochromatin-related genes (see Materials and Methods). We found that signatures of positive selection are present in domains with diverse functions (Fig. 3), and several of them are especially pertinent to heterochromatin functions. These include chromo domains that directly interact with histone methylation (in HP6, rhi, and Kis) and Jmjc demethylase domain in NO66, the H3K4me3 demethylase. We also found signatures of positive selection located within various nucleic acid-binding domains (e.g. C2H2 type and UBZ4-type Zinc-Finger DNA-binding domain) and domains mediating protein–protein interactions (e.g. BESS domains).

Fig. 3. Signatures of positive selection fall within domains with diverse functions and IDRs of heterochromatin-related genes. Genes with evidence of positive selection on protein sequences are categorized according to which types of domains/regions their signatures of positive selection fall within. Annotated and/or predicted structured domains are shown as horizontal lines with the following functional categories: nucleic acid binding (orange), protein–protein interaction (yellow), histone modification binding (pink), and other domains (purple). IDRs are marked with gray. Highlighted vertical windows represent those under positive selection as identified by sliding MK tests (green) and PAML inferences (blue). Note that genes with no positively selected windows/sites overlapping with either annotated domains or IDRs were not shown.

In addition to structured domains with well-characterized functions, proteins also contain regions that lack fixed 3D structures, known as intrinsically disordered regions (IDRs). IDRs are increasingly appreciated for playing critical roles in protein functions (Forman-Kay and Mittag 2013) and are frequently found in proteins enriched in phase-separated cellular compartments (Nott et al. 2015; Pak et al. 2016; Lin et al. 2017). Indeed, IDR-mediated phase properties have been suggested to be critical for heterochromatin functions (Larson et al. 2017; Strom et al. 2017). Accordingly, we investigated whether IDRs in heterochromatin-related genes may also have signatures of positive evolution by using flDPnn (Hu et al. 2021) to predict IDRs in D. melanogaster protein sequences. Interestingly, we found that signatures of rapid evolution also frequently fall within IDRs (Fig. 3), and, for multiple rapidly evolving heterochromatin-related genes, fast-evolving sites/windows only fall within IDRs, but not other structured domains (Fig. 3, right column).

We also investigated whether varying IDR properties of heterochromatin-related genes contribute to the evolutionary differences between gene categories by estimating the percentage of amino acids falling within predicted IDR domains (percentage of IDRs). A higher percentage of IDRs has been found in proteins involved in the formation of phase-mediated cellular domains (e.g. membrane-less organelles [Sawyer et al. 2019] and the heterochromatin domain [Guthmann et al. 2019]). We found that the percentage of IDR for heterochromatin-related genes in D. melanogaster is significantly greater than that of random genes (Mann–Whitney U test, P = 0.00214), which is especially true for Su(var) and heterochromatin proteins (Mann–Whitney U test, P = 0.0448 [Su(var)s] and 0.0141 [heterochromatin proteins]; supplementary fig. S4a, Supplementary Material online). There is a lack of difference in percentage of IDR between our candidate genes and Polycomb control (Mann–Whitney U test, P > 0.05 for all tests, supplementary fig. S4a, Supplementary Material online), which may be due to the fact that Polycomb proteins were shown to also undergo phase separation (Tatavosian et al. 2019). Surprisingly, despite previous reports on the reduced evolutionary constraints in IDRs (Brown et al. 2002; Khan et al. 2015), the percentage of IDRs does not differ between heterochromatin-related genes with or without evidence of positive selection (Mann–Whitney U test, P > 0.05; supplementary fig. S4b, Supplementary Material online).

To further investigate the possible associations between the percentage of IDRs and the evolution of candidate genes, we estimated the phylogenetic signals (Blomberg and Garland 2002) of the percentage of IDRs across studied species. A small Blomberg's K (Blomberg et al. 2003), our chosen index for phylogenetic signal, suggests the percentage of IDR evolves faster than that of random expectation (Brownian motion of trait evolution; Kamilar and Cooper 2013). Intriguingly, we found no difference in K values between heterochromatin-related genes and control genes, among categories of heterochromatin-related genes that showed different rates of adaptive evolution, and between heterochromatin-related genes with and without evidence of positive selection over long evolutionary time scales (supplementary fig. S4c and d, Supplementary Material online; Mann–Whitney U test, P > 0.05 for all comparisons). These observations indicate that, despite the observed pervasive adaptive evolution of heterochromatin-related genes, their percentage of IDRs evolves similarly to that of other genes in the genome. In summary, positive selection acts on both structured domains with diverse functions and IDRs of heterochromatin-related genes, while stabilizing selection might have maintained the IDR content of these genes in the face of rapidly changing amino acid sequences.

Rates of Protein Evolution of Heterochromatin-Related Genes Significantly Associate with Genomic Abundance of Transposable Elements

The antagonistic interaction with repetitive sequences enriched in heterochromatin was suggested as the driver for the rapid evolution of a handful of heterochromatin-related genes (Vermaak et al. 2005; Satyaki et al. 2014; Helleu and Levine 2018). To test whether this conjecture may be broadly applicable to heterochromatin-related genes, we estimated the correlation between the rates of amino acid substitution (dN/dS) of heterochromatin-related genes and the changes in abundance of heterochromatic repetitive sequences using methods developed in Lartillot and Poujol (2011). This approach corrects for the phylogenetic nonindependence of the quantitative traits of interest (here, repeat abundance) and models their evolution following the Brownian process. To estimate repeat abundance, we performed Illumina sequencing with PCR-free library preparation to avoid the sequencing bias against AT-rich sequences, which is commonly found in repetitive sequences (Wei et al. 2018). Using these sequencing data, we quantified the genomic abundance of satellite repeats (both simple and complex satellites) and TEs, both of which should be dominated by the DNA content of heterochromatin (see Materials and Methods).

Among the 42 heterochromatin-related genes with evidence of positive selection over the long evolutionary time scale, nearly 40% (38.1%) of them have rates of protein evolution tracking changes in TE abundance across species, with significant (P < 0.05) phylogenetically controlled correlations between dN/dS and TE abundance (supplementary table S3, Supplementary Material online). This proportion is much higher than that of the 248 random genes with accelerated rates of protein evolution (20.5%, Fisher's exact test, P = 0.0173; binomial test, P = 0.0112; Fig. 4a). Intriguingly, 93.8% of these correlations of the heterochromatin-related genes are in the negative direction, which indicates that these genes evolved faster in species with lower TE abundance (see Fig. 4c for an example). We also examined this proportion separately for each category of heterochromatin-related genes (Fig. 4b). Quite strikingly, even though only a small number of histone-modifying enzymes were found to have evidence of positive selection (26.1%), all of their rates of protein evolution were significantly associated with TE abundance (Fig. 4b). These histone-modifying enzymes include all three H3K9me2/3 methyltransferase (egg, G9a, and Su(var)3-9), NO66 (H3K4me3 demethylase), KDM4B (H3K9me2/3 demethylase), and HDAC6 (zinc-dependent deacetylase) (supplementary table S3, Supplementary Material online). In stark contrast to analysis focusing on TE abundance, only 7.14% of heterochromatin-related genes with evidence of positive selection show a significant correlation between dN/dS and the abundance of total satellite repeats, a proportion that is not significantly different from that of randomly sampled genes (Fig. 4a; Fisher's exact test, P = 0.322 and binomial test, P = 0.360). Analysis focusing on the abundance of only simple satellites gave consistent results (supplementary fig. S5, Supplementary Material online, P > 0.05 for both Fisher's exact test and binomial test). Overall, we found that the rates of protein evolution of heterochromatin-related genes negatively correlated with the abundance of TEs, but not total repeats, across species, suggesting a possible source of selective force shaping the evolution of these genes (see Discussions).

Fig. 4. Significant negative associations between dN/dS of heterochromatin-related genes and genomic TE abundance. a) Stacked bar plots showing the proportion of positively selected genes whose dN/dS correlates with repeat abundance (left) or TE abundance (right) for heterochromatin-related genes and randomly sampled genes. The color of the genes in the side box represents the gene's category: E(var)s (yellow), Su(var)s (green), histone-modifying enzyme enhancing heterochromatin (light green), and heterochromatin protein (dark green). b) Bar plots showing the proportion of positively selected heterochromatin-related genes whose dN/dS significantly coevolve with TE abundance for different gene categories. c) X-Y plots showing the associations between dN/dS (X-axis) and normalized TE abundance for an example gene (sov). Blue points are estimates for the extant species, and gray points are for internal nodes. Lines across each dot denote 95% confidence intervals for dN/dS (extant species and internal nodes) and normalized TE content (internal nodes). Binomial test: n.s. P > 0.05; *P < 0.05.

Discussion

Heterochromatin is a highly conserved cellular compartment with essential functions across complex eukaryotes (Allshire and Madhani 2018; Janssen et al. 2018). Nevertheless, our evolutionary analyses revealed that genes involved in heterochromatin function are highly labile, exhibiting pervasive evidence of rapid evolution both in the forms of amino acid substitutions and gene copy number changes at both short and long evolutionary time scales. Importantly, the rapid evolution of these genes is likely driven by functions specific to constitutive heterochromatin, instead of mechanisms general to proteins interacting with repressive chromatin. Evidence of positive selection on protein evolution is especially prominent for heterochromatin-related genes that should enhance heterochromatin function, and, strikingly, all three methyltransferases responsible for the enrichment of H3K9me2/3, the characteristic histone modification of heterochromatin, are under positive selection. Our further characterization of the various aspects of the evolution of heterochromatin-related genes provided an important avenue to identify the possible source of selective forces acting on heterochromatin, which we discussed below.

Close examinations of the signatures of positive selection of heterochromatin-related genes found them not only located within structured domains with well-known functions but also inside unstructured IDRs (Fig. 3). IDRs were recently found to be critical for the phase properties and thus functions of heterochromatin (Larson et al. 2017; Strom et al. 2017), and our observation suggests that positive selection could also act on heterochromatin functions mediated by such properties. Interestingly, despite the previous suggestions that IDRs are evolutionarily less constrained (Brown et al. 2002; Khan et al. 2015), the proportion of IDR sequences in positively selected heterochromatin-related genes is similar to other candidate genes (Fig. 4b). Even more, we identified strong phylogenetic signals that are consistent with the presence of stabilizing selection preserving the percentage of IDRs even for heterochromatin-related genes under strong positive selection. Such findings may indicate that the percentage of IDRs of a protein, but not the underlying amino acid sequence, is evolutionarily constrained and plays an essential role in the function of studied heterochromatin-related genes. Indeed, a few IDRs with rapidly evolving sequences were found to possess conserved molecular features (Moesa et al. 2012; Zarin et al. 2019), with some of them experimentally demonstrated to be functionally equivalent (Zarin et al. 2017). This raises another question of why the underlying amino acid sequences rapidly evolve, considering the need to maintain the IDR content of heterochromatin-related genes (see below).

Several of our findings suggest that TEs, but not satellite repeats, must be involved in the pervasive rapid evolution of heterochromatin-related genes. First, the associations between rates of protein evolution of heterochromatin-related genes and repeat abundance were mainly observed for TEs. In addition, many heterochromatin-related genes under positive selection have functions related to TE suppression. In particular, the interactions between the protein products of several of these genes (e.g. del, moon, rhi, Trf2, and sov; Klattenhoff et al. 2009; Mohn et al. 2014; Andersen et al. 2017; Andreev et al. 2022) and repressive H3K9me2/3 are responsible for licensing piRNA clusters, which are TE-enriched loci generating the majorities of piRNAs targeting TEs and located within pericentromeric heterochromatin. Similar to the rapid evolution of the DNA and protein components of centromeres (i.e. the centromeric drive hypothesis; Henikoff et al. 2001; Malik et al. 2002), changes in heterochromatic TE sequences may alter their interactions with proteins encoded by heterochromatin-related genes. Consequently, selection may favor evolutionary changes in proteins that revert the strength of DNA–protein interaction to that prior to the changes of TE sequences. Subsequent changes in TE sequences could initiate another cycle of this interaction and ultimately drive the fast evolution of the genes involved. However, only some of the identified signatures of positive selection fall within nucleic acid-binding domains (Fig. 3), and the majority of the positively selected heterochromatin-related genes lack sequence specificity. Moreover, unlike viruses, the “success” of TEs is tightly intertwined with host fitness, owing to the fact that they propagate in the host germline and are inherited vertically (Haig 2016). Accordingly, an arms race between TEs and host proteins suppressing them was suggested to be unlikely to drive the pervasive adaptive evolution of host genes (Blumenstiel et al. 2016; Cosby et al. 2019).

What might have driven the fast evolution of heterochromatin-related genes then? The genomic autoimmunity hypothesis (Blumenstiel et al. 2016) provides a plausible explanation—molecular mechanisms suppressing TE activities are expected to be constantly juggling between maximizing TE silencing while minimizing inadvertent off-target suppression of functional genes, and this alternating selective pressure could drive the rapid evolution of genes involved. Under this model, for example, variants of heterochromatin-related genes that enhance the generation of piRNAs from piRNA clusters, and thus strengthen TE silencing, may lead to the inadvertent production of piRNAs that target other functional elements, a possibility with empirical support (Andersen et al. 2017; Kelleher 2021). Selection will then instead favor variants increasing the stringency of this licensing to avoid fitness costs of the off-target effect. The repeated alternation of selective targets could accordingly drive the rapid evolution of heterochromatin-related genes.

The genomic autoimmunity hypothesis is even more suitable to explain the rapid evolution of heterochromatin-related genes involved in modulating the intrinsic and unique molecular properties of heterochromatin—a tendency to “spread” to nearby loci in a sequence-independent manner. The positive feedback between writers and readers of H3K9me2/3 promotes the propagation of repressive chromatin marks (Allshire and Madhani 2018; Bell et al. 2023), and the extent of this spreading depends on the concentrations of readers and writers at the suppressed loci (Locke et al. 1988). Heterochromatin-mediated silencing thus needs to be carefully balanced to prevent inadvertent silencing of functional elements while maintaining sufficient suppression at the euchromatin-heterochromatin boundaries. Interestingly, such a balance needs to be maintained not only at the ends of chromosomes (around pericentromeric and subtelomeric heterochromatin) but also around epigenetically silenced, H3K9me2/3-enriched TEs in the euchromatic genome. The spreading of repressive marks from silenced euchromatic TEs into functional genes is widely documented (reviewed in Choi and Lee 2020) and has been inferred to impair individual fitness (Lee 2015; Lee and Karpen 2017; Huang et al. 2022). Consistently, several genes known to directly mediate TE epigenetic silencing are found to be under positive selection in our analysis (e.g. wde and Su(var)2-10; Ninova et al. 2020). In addition, histone-modifying enzymes, which are directly involved in the reader–writer positive feedback loops, should frequently be caught in cycles of alternating selection for enhanced or weakened TE epigenetic silencing. Indeed, we found that all three writers for H3K9me2/3 show evidence of rapid evolution (Su(var)3-9, egg, and G9a) and all histone-modifying enzymes with accelerated rates of protein evolution have evidence of coevolution with genomic TE abundance (Fig. 4). Furthermore, the deleterious off-target effects of TE silencing could result from the long-range spatial interactions between H3K9me2/3-enriched euchromatic TEs and pericentromeric heterochromatin, which is mediated by phase separation mechanisms and is selected against (Lee, Ogiyama, et al. 2020). Repeated alteration in selective pressure to ensure these proteins confer sufficient phase properties for proper heterochromatin functions while avoiding such off-target effects may similarly drive rapid changes in protein sequences while preserving the percentage of IDR for proteins encoded by heterochromatin-related genes, as we have observed.

Intriguingly, we found predominantly negative associations between TE abundance and rates of protein evolution of heterochromatin-related genes, which may initially seem counterintuitive. However, this pattern can be explained by our proposed model that augments the genomic autoimmunity hypothesis by incorporating how the alternating selective pressure on heterochromatin-mediated silencing may concurrently influence genomic TE abundance (Fig. 5a). Strong heterochromatin-mediated silencing should lead to rampant off-target effects, which not only select for variants that weaken the strength of silencing (Huang and Lee 2024) but also select against individual TE copies due to the deleterious spreading of repressive marks to TE-adjacent functional sequences. TE abundance should accordingly decrease (Fig. 5a1; Huang et al. 2022). When variants that reduce heterochromatin-mediated silencing become fixed, the maintenance of heterochromatin at epigenetically silenced TEs weakens, resulting in increased TE replication and more new TE insertions (Fig. 5a2). The consequentially increased TE abundance should then drive selection for enhanced heterochromatin-mediated silencing, returning to the initial state of strong silencing (Fig. 5b). Accordingly, as heterochromatin-related genes experience cycles of alternating selection targets and gaining amino acid substitutions, genomic TE abundance also fluctuates. Yet, whether TE abundance eventually increases or decreases depends on the relative strength of selection against TEs with off-target effects and the changes in TE replication rates (Fig. 5b). Selection coefficients for the harmful off-target effects of TE epigenetic silencing are yet to be estimated, but they could be strong if such effects perturb the expression of nearby vital genes (e.g. Coronado-Zamora and González 2023) or disrupt global 3D genome organization (e.g. Lee, Ogiyama, et al. 2020). On the other hand, replication rates of Drosophila TEs are generally low (10−5 to 10−4; Charlesworth and Langley 1989; Adrion et al. 2017; Wang et al. 2023). Changes in these rates are likely weaker than selection removing TEs through their off-target effects, leading to decreased TE abundance over cycles of alternating selective pressure on the strength of heterochromatin-mediated silencing and thus giving rise to negative associations between genomic TE abundance and rates of protein evolution between species. It is worth noting that many other processes, such as recent demographic changes (e.g. Mérel et al. 2021), also contribute to between-species differences in TE abundance. If these forces systematically influence genomic TE abundance in the species studied (e.g. correlate with the evolution of heterochromatin-related genes), similar associations could arise in the absence of the proposed model.

Fig. 5. Model for the correlated evolution of heterochromatin-related genes and genomic TE abundance. a) The predicted impacts of different strengths of heterochromatin-mediated silencing on TE epigenetic silencing and genomic TE abundance. a, 1) Strong heterochromatin-mediated silencing results in strong off-target effects due to the spreading of repressive marks from silenced euchromatic TEs, leading to positive selection for variants of heterochromatin-related genes that weaken the strength of silencing. In addition, because of the low rates of TE replication from strong silencing and negative selection against TEs with off-target effects, TE abundance should reduce. a, 2) On the other hand, when heterochromatin-mediated silencing is weak, selection against TE-mediated off-target effects is weak and rates of TE replication should be high, leading to increased TE abundance. At the same time, selection may favor variants that enhance TE silencing to reduce rates of TE replication. b) Proposed sequences of events that lead to correlated evolution of heterochromatin-related genes and genomic TE abundance. The model should be applicable to any of the states and the following description starts with strong heterochromatin-mediated silencing (i). TE abundance is expected to reduce due to the suppressed rates of TE replication and selection against TEs with off-target effects (Phase 1). Concurrently, there should be selection for variants of heterochromatin-related genes that weaken heterochromatin-mediated silencing to reduce deleterious off-target effects. When variants that weaken silencing are fixed, the maintenance of TE epigenetic silencing decreases (ii), leading to both reduced occurrence of off-target effects and increased rates of TE replication. TE abundance should thus increase (Phase 2). The high TE abundance and replicative activity would select for variants of heterochromatin-related genes that enhance TE silencing, going back to the initial state of strong heterochromatin-mediated silencing (e.g. iii here). The relative strength of a decrease in TE abundance through selection against off-target effects and increase in TE abundance through TE replication determines whether there are negative (Scenario 1) or positive (Scenario 2) associations between TE abundance and the evolution of heterochromatin-related genes across species.

Our proposed model, which suggests that the necessity of how the need of maintaining a balanced TE silencing drives the fast evolution of heterochromatin-related genes, may help explain some intriguing patterns we observed. The majority of our identified duplicates of heterochromatin-related genes (48.78%) have male-biased expression (Fig. 2), indicating their potential functional significance in male germlines. While the resolution of sexual genetic conflict (Gallach and Betrán 2011; Wyman et al. 2012), a common explanation for male-biased expression of newly duplicated genes, may underlie such an observation, it may also stem from the need to properly silence the male-specific TE landscape (Chen et al. 2021), resulting in sexually dimorphic chromatin environment. Another intriguing observation is the limited evidence of positive selection for heterochromatin-antagonizing genes (Fig. 1), which is in stark contrast to the pervasive signatures of adaptive evolution of genes mediating heterochromatin functions. Protein products of heterochromatin-antagonizing genes are not only found at heterochromatin–euchromatin boundaries but broadly distributed throughout the euchromatic genome (Dimova et al. 2003; Kellner et al. 2012), playing diverse roles in transcriptional regulation (Meignin and Davis 2008; Regnard et al. 2011; Crona et al. 2013). The pleiotropic functional roles of heterochromatin-antagonizing proteins, in both euchromatin and heterochromatin, should constrain their evolution, making them less likely to be caught up in alternating selection for the sensitivity or the specificity of TE silencing. Nevertheless, the small sample size of heterochromatin-antagonizing genes (only 12 genes) could limit our ability to detect positive selection among these genes.

Our observed frequent positive selection acting on heterochromatin-related genes points out the high lability of the molecular feature of this functionally conserved genomic compartment. Instead of the usually assumed “arms race” with the changing repeatome, we proposed that the selective pressure on heterochromatin-related genes might mainly come from a need to maintain a delicate balance of its unique ability to “suppress” and to “spread,” which also consequently influence the evolutionary dynamics of TEs. Signatures of positive selection identified here could serve as an entry point to further investigate how the delicate balance of heterochromatin-mediated silencing may be conferred by vastly changing components between species (e.g. Rosin and Mellone 2016; Parhad et al. 2017), providing a fruitful opportunity to further dissect the molecular mechanisms shaping heterochromatin functions and evolution.

Materials and Methods

Evolutionary Analyses of Protein Sequences

Coding sequences and genome annotations for 16 studied species (Fig. 1b) were downloaded from NCBI, with GenBank ID listed in supplementary table S4, Supplementary Material online. Because we compiled the list of candidate genes based on D. melanogaster literature, we used orthologous information from OrthoDB (Kuznetsov et al. 2023; last accessed 12/10/2022) to retrieve one-to-one orthologs for D. melanogaster candidate genes. For genes without one-to-one ortholog according to OrthoDB, we performed BLAST search (see below). The retrieved coding sequences were translated to amino acid sequences, aligned using Clustal Omega (version 1.2.4; Sievers et al. 2011), and converted back to codon alignments using PAL2NAL (version 14; Suyama et al. 2006).

We performed unpolarized MK tests (McDonald and Kreitman 1991) using polymorphism within a D. melanogaster Zambian population (197 strains; Lack et al. 2015) and divergence between D. melanogaster and D. simulans. Following the recommendation of Lack et al. (2015), we treated genomic regions with non-African ancestry or identity-by-descent as missing data and only included sites with at least 100 nonmissing alleles. To count the number of nonsynonymous and synonymous changes, we used the mutational paths that minimize the number of amino acid changes. Codons with more than two alleles, considered both within species polymorphism and between species divergence, were excluded. Genes with fewer than 100 codons were excluded from the analysis due to insufficient statistical power. A gene is deemed under positive selection if the MK test, whose significance was assessed using Fisher's exact test, rejected the null hypothesis and the ratio of nonsynonymous to synonymous changes is greater for between-species substitutions than for within-species polymorphism. Sliding window MK tests were performed with windows of 100 codons and 10-codon steps.

We conducted PAML (Yang 2007) using 16 species to identify candidate genes experiencing positive selection over a long evolutionary time scale. We compared two site models, M8a (dN/dS ratio is beta-distributed and not greater than one) and M8 (dN/dS > 1), and determined the significance using likelihood ratio tests. The species tree reported in Suvorov et al. (2022) was used. Sites with >0.95 BEB posterior probability of coming from the site class with dN/dS > 1 (Yang et al. 2005) are considered under recurrent adaptive evolution.

Evolutionary Analysis for Gene Copy Number

To identify genes with varying gene copy numbers, we first used a genome-wide, high throughput search with liberal parameters to identify many potential candidates, followed by careful manual curations. We first used tblastn and reciprocal blastx (Camacho et al. 2009) to search for homologs and paralogs of candidate genes in studied species using D. melanogaster amino acid sequence as queries, with the following parameter: e < 10−2, amino acid identity > 20%, and blast score > 10. We required the best reciprocal blastx hit to be the original D. melanogaster query. Each potential CNV was manually validated using reciprocal best blast with more stringent criteria (e < 10−5), and orthologs and paralogs were distinguished using synteny of flanking genes (supplementary fig. S1, Supplementary Material online). We further examined the expression levels of candidate CNVs using RNA-seq exon coverage tracks of NCBI Data Viewer and removed those with no expression. Several duplicates identified in D. suzukii were filtered due to redundant contigs. For gene loss, we followed the procedures detailed in King et al. (2019) to confirm the true absence of a gene.

To examine the sex-biased expression, we deemed a gene male biased if the log2 fold change of the ratio of male and female expression (transcripts per million or fragments per kilobase of transcript per million fragments mapped) is >1, female biased if such value is <−1, and otherwise unbiased. For D. melanogaster, we used Insect Genome database (http://www.insect-genome.com/Sexdb/, last retrieved November 2023) and FlyAtlas2 RNA-seq data (Krause et al. 2022) for whole adult males and whole adult females. For the other 15 species, we mapped publicly available transcriptome data sets (supplementary table S5, Supplementary Material online) to publicly available genome assemblies with gene annotations using HiSAT2 (v2.2.1 with parameters –exon and –ss to specify the exon positions and splice sites; RRID:SCR_015530; Kim et al. 2019) and estimated the expression levels usingStringtie (v2.1.4 with parameters -dta -G to specify annotation files; RRID:SCR_016323; Kovaka et al. 2019). Candidate duplicates identified by manual curation but have no annotation or show no expression were excluded from the analysis.

Domain Annotations

We used UniProt (The UniProt Consortium 2023) to annotate known structured domains within the D. melanogaster allele of heterochromatin-related genes. We used flDPnn (Hu et al. 2021), which performed superior in the previous benchmark study (Necci et al. 2021) and the predicted binary index for IDRs (disorder propensity cutoff = 0.3) to annotate IDRs for all studied species. To detect the phylogenetic signal of the % of IDR, we computed Blomberg's K statistic using the phylosig function from phytools R package (Blomberg et al. 2003; Revell 2012). The tree structure and branch lengths were obtained from Suvorov et al. (2022) and generated by treeio R package (Wang et al. 2020).

Analyses of Repetitive Sequences and Their Coevolution with Candidate Genes

DNA for each studied species was extracted from 40 females using DNeasy Blood & Tissue Kit (Qiagen), following the manufacturer's protocol. To avoid PCR amplification bias during the preparation of sequencing libraries, which was found to skew the quantification of repeats (Wei et al. 2018), extracted DNA was prepared into Illumina sequencing library with PCR-free protocol and sequenced with 150-bp paired-end reads by Novogene (Sacramento, CA, USA).

We used Satellite Repeat Finder (Zhang et al. 2023) to estimate the abundance of total satellite repeat in the Illumina short-read sequencing data. Following the suggested procedures, we first counted K-mers (K = 21) in each sample using KMC (ver. 3.2.4; Kokot et al. 2017). Contigs of satellites were then generated and mapped back to the source Illumina sequences to estimate the abundance (in bases) of each satellite using minimap (ver. 2.24; Li 2018) and Satellite Repeat Finder. We also used K-seek (Wei et al. 2014) to estimate the abundance of simple repeats and reached similar conclusions (see Results).

TE abundance was estimated as the total TE reads from the Illumina sequencing data, with the assumption that TEs in the heterochromatin mainly originated from jumping events of euchromatic TEs and TE abundances in the heterochromatin and euchromatin are highly correlated. To minimize TE annotation bias across species, repetitive sequences from each of the 16 genomes (supplementary table S4, Supplementary Material online) were identified using RepeatMasker (version 4.1.0; http://www.repeatmasker.org/) with the Dfam database (Storer et al. 2021), using the command “RepeatMasker -q [genome sequence file] -species drosophila -e hmmer”. TE sequences annotated as LTR, LINE, DNA element, and unknown categories were obtained from all 16 genomes to create a master TE library. Illumina reads from each species were then mapped to the library with bwa-mem (version 0.7.17; Li and Durbin 2009) and viewed by samtools (version 1.15.1; Li 2011). The total number of reads mapped to the library, regardless of mapping quality, was considered the number of TE reads. It is worth noting that very few TEs were classified as unknown category, and inclusion/exclusion of such TEs did not change the results. To compare the abundance of satellite repeats and TEs across samples/species, these estimates were normalized. Illumina reads from each species were mapped to its repeat-masked genome using bwa-mem, with sites with mapping quality lower than 30 filtered (using samtools -q 30). The median read depth for the unmasked regions for each sample was then used to normalize the number of bases for satellite repeats and TEs.

We used coevolve (Lartillot and Poujol 2011) to estimate the correlation between dN/dS and repeat abundance (satellite repeats or TEs) given the tree structure. For each gene, the analysis was performed in duplicate to assess convergence (relative difference < 0.1) and a burn-in of 300 with at least 3,000 MCMC chains to get the final estimated correlation between dN/dS and repeat abundance. Because the number of Polycomb genes with significant PMAL tests is small, we only compared the coevolution results of candidate genes to those of 248 randomly chosen genes with accelerated rates of amino acid evolution (i.e. significant for PAML analysis).

Supplementary Material

msae181_Supplementary_Data

Acknowledgments

We thank the High Performance Cluster at UC Irvine for computational resources and Harsh Shukla for helping with the flDPnn software. We appreciate Aniek Janssen, Andrea Betancourt, and Adriana Ludwig for their helpful comments on the manuscript.

Supplementary Material

Supplementary material is available at Molecular Biology and Evolution online.

Funding

C.-H.C. was supported by Damon-Runyon Cancer Research Foundation postdoctoral fellowship DRG 2438-21 and National Institute of Health R01-GM74108 (to Hamit S. Malik). S.U.C. was supported by National Institute of Health R35GM139653 (to Gary H. Karpen). Y.H., L.L., J.M., and Y.C.G.L. were supported by National Institute of Health R35-GM14292 (to Y.C.G.L.).

Data Availability

PCR-free Illumina data have been deposited to SRA under the accession number PRJNA1113679. Scripts used in this study can be found at https://github.com/YuhengHuang87/HC_gene_evo and https://github.com/hilynano/Homology_detection.
==== Refs
References

Adrion  JR, Song  MJ, Schrider  DR, Hahn  MW, Schaack  S. Genome-wide estimates of transposable element insertion and deletion rates in Drosophila melanogaster. Genome Biol Evol. 2017:9 (5 ):1329–1340. 10.1093/gbe/evx050.28338986
Alekseyenko  AA, Gorchakov  AA, Zee  BM, Fuchs  SM, Kharchenko  PV, Kuroda  MI. Heterochromatin-associated interactions of Drosophila HP1a with dADD1, HIPP1, and repetitive RNAs. Genes Dev. 2014:28 (13 ):1445–1460. 10.1101/gad.241950.114.24990964
Allshire  RC, Madhani  HD. Ten principles of heterochromatin formation and function. Nat Rev Mol Cell Biol.  2018:19 (4 ):229–244. 10.1038/nrm.2017.119.29235574
Andersen  PR, Tirian  L, Vunjak  M, Brennecke  J. A heterochromatin-dependent transcription machinery drives piRNA expression. Nature. 2017:549 (7670 ):54–59. 10.1038/nature23482.28847004
Andreev  VI, Yu  C, Wang  J, Schnabl  J, Tirian  L, Gehre  M, Handler  D, Duchek  P, Novatchkova  M, Baumgartner  L, et al  Panoramix SUMOylation on chromatin connects the piRNA pathway to the cellular heterochromatin machinery. Nat Struct Mol Biol. 2022:29 (2 ):130–142. 10.1038/s41594-022-00721-x.35173350
Bannister  AJ, Kouzarides  T. Regulation of chromatin by histone modifications. Cell Res.  2011:21 (3 ):381–395. 10.1038/cr.2011.22.21321607
Bao  X, Deng  H, Johansen  J, Girton  J, Johansen  KM. Loss-of-function alleles of the JIL-1 histone H3S10 kinase enhance position-effect variegation at pericentric sites in Drosophila heterochromatin. Genetics. 2007:176 (2 ):1355–1358. 10.1534/genetics.107.073676.17435241
Barbash  DA, Siino  DF, Aaron  M, Roote  J. A rapidly evolving MYB-related protein causes species isolation in Drosophila. Proc Natl Acad Sci U S A.  2003:100 (9 ):5302–5307. 10.1073/pnas.0836927100.12695567
Bell  O, Burton  A, Dean  C, Gasser  SM, Torres-Padilla  M-E. Heterochromatin definition and function. Nat Rev Mol Cell Biol. 2023:24 (10 ):691–694. 10.1038/s41580-023-00599-7.37069331
Blomberg  SP, Garland  T  Jr. Tempo and mode in evolution: phylogenetic inertia, adaptation and comparative methods. J Evol Biol.  2002:15 (6 ):899–910. 10.1046/j.1420-9101.2002.00472.x.
Blomberg  SP, Garland  T, Ives  AR. Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution. 2003:57 (4 ):717–745. 10.1111/j.0014-3820.2003.tb00285.x.12778543
Blumenstiel  JP, Erwin  AA, Hemmer  LW. What drives positive selection in the Drosophila piRNA machinery? The genomic autoimmunity hypothesis. Yale J Biol Med. 2016:89 :499–512. https://pubmed.ncbi.nlm.nih.gov/28018141/.28018141
Brand  CL, Levine  MT. Functional diversification of chromatin on rapid evolutionary timescales. Annu Rev Genet.  2021:55 (1 ):401–425. 10.1146/annurev-genet-071719-020301.34813351
Brand  CL, Oliver  GT, Farkas  IZ, Buszczak  M, Levine  MT. Recurrent duplication and diversification of a vital DNA repair gene family across Drosophila. Mol Biol Evol.  2024:41 (6 ):msae113. 10.1093/molbev/msae113.38865490
Brown  CJ, Takayama  S, Campen  AM, Vise  P, Marshall  TW, Oldfield  CJ, Williams  CJ, Dunker  AK. Evolutionary rate heterogeneity in proteins with long disordered regions. J Mol Evol. 2002:55 (1 ):104–110. 10.1007/s00239-001-2309-6.12165847
Camacho  C, Coulouris  G, Avagyan  V, Ma  N, Papadopoulos  J, Bealer  K, Madden  TL. BLAST+: architecture and applications. BMC Bioinformatics. 2009:10 (1 ):421. 10.1186/1471-2105-10-421.20003500
Charlesworth  B, Langley  CH. The population genetics of Drosophila transposable elements. Annu Rev Genet.  1989:23 (1 ):251–287. 10.1146/annurev.ge.23.120189.001343.2559652
Chen  P, Kotov  AA, Godneeva  BK, Bazylev  SS, Olenina  LV, Aravin  AA. piRNA-mediated gene regulation and adaptation to sex-specific transposon expression in D. melanogaster male germline. Genes Dev. 2021:35 (11-12 ):914–935. 10.1101/gad.345041.120.33985970
Choi  JY, Lee  YCG. Double-edged sword: the evolutionary consequences of the epigenetic silencing of transposable elements. PLoS Genet.  2020:16 (7 ):e1008872. 10.1371/journal.pgen.1008872.32673310
Coronado-Zamora  M, González  J.  2023. The epigenetics effects of transposable elements are context dependent and not restricted to gene silencing:2023.11.27.568862. Available from: https://www.biorxiv.org/content/10.1101/2023.11.27.568862v1.
Cosby  RL, Chang  N-C, Feschotte  C. Host–transposon interactions: conflict, cooperation, and cooption. Genes Dev. 2019:33 (17-18 ):1098–1116. 10.1101/gad.327312.119.31481535
Crona  F, Dahlberg  O, Lundberg  LE, Larsson  J, Mannervik  M. Gene regulation by the lysine demethylase KDM4A in Drosophila. Dev Biol.  2013:373 (2 ):453–463. 10.1016/j.ydbio.2012.11.011.23195220
Czermin  B, Schotta  G, Hülsmann  BB, Brehm  A, Becker  PB, Reuter  G, Imhof  A. Physical and functional association of SU(VAR)3-9 and HDAC1 in Drosophila. EMBO Rep. 2001:2 (10 ):915–919. 10.1093/embo-reports/kve210.11571273
de Lima  LG, Ruiz-Ruano  FJ. In-depth satellitome analyses of 37 Drosophila species illuminate repetitive DNA evolution in the Drosophila genus. Genome Biol Evol.  2022:14 (5 ):evac064. 10.1093/gbe/evac064.35511582
Deng  H, Bao  X, Cai  W, Blacketer  MJ, Belmont  AS, Girton  J, Johansen  J, Johansen  KM. Ectopic histone H3S10 phosphorylation causes chromatin structure remodeling in Drosophila. Development. 2008:135 (4 ):699–705. 10.1242/dev.015362.18199578
Dimova  DK, Stevaux  O, Frolov  MV, Dyson  NJ. Cell cycle-dependent and cell cycle-independent control of transcription by the Drosophila E2F/RB pathway. Genes Dev. 2003:17 (18 ):2308–2320. 10.1101/gad.1116703.12975318
Eissenberg  JC, Elgin  SCR. HP1a: a structural chromosomal protein regulating transcription. Trends Genet.  2014:30 (3 ):103–110. 10.1016/j.tig.2014.01.002.24555990
Elgin  SCR, Reuter  G. Position-effect variegation, heterochromatin formation, and gene silencing in Drosophila. Cold Spring Harb Perspect Biol. 2013:5 (8 ):a017780. 10.1101/cshperspect.a017780.23906716
Feller  C, Forné  I, Imhof  A, Becker  PB. Global and specific responses of the histone acetylome to systematic perturbation. Mol Cell.  2015:57 (3 ):559–571. 10.1016/j.molcel.2014.12.008.25578876
Feng  W, Michaels  SD. Accessing the inaccessible: the organization, transcription, replication, and repair of heterochromatin in plants. Annu Rev Genet.  2015:49 (1 ):439–459. 10.1146/annurev-genet-112414-055048.26631514
Flores  HA, DuMont  VLB, Fatoo  A, Hubbard  D, Hijji  M, Barbash  DA, Aquadro  CF. Adaptive evolution of genes involved in the regulation of germline stem cells in Drosophila melanogaster and D. simulans. G3 (Bethesda). 2015:5 (4 ):583–592. 10.1534/g3.114.015875.25670770
Forman-Kay  JD, Mittag  T. From sequence and forces to structure, function, and evolution of intrinsically disordered proteins. Structure. 2013:21 (9 ):1492–1499. 10.1016/j.str.2013.08.001.24010708
Gallach  M, Betrán  E. Intralocus sexual conflict resolved through gene duplication. Trends Ecol Evol. 2011:26 (5 ):222–228. 10.1016/j.tree.2011.02.004.21397976
Garzino  V, Pereira  A, Laurenti  P, Graba  Y, Levis  RW, Le Parco  Y, Pradel  J. Cell lineage-specific expression of modulo, a dose-dependent modifier of variegation in Drosophila. EMBO J. 1992:11 (12 ):4471–4479. 10.1002/j.1460-2075.1992.tb05548.x.1425581
Ghildiyal  M, Seitz  H, Horwich  MD, Li  C, Du  T, Lee  S, Xu  J, Kittler  ELW, Zapp  ML, Weng  Z, et al  Endogenous siRNAs derived from transposons and mRNAs in Drosophila somatic cells. Science. 2008:320 (5879 ):1077–1081. 10.1126/science.1157396.18403677
Girton  JR, Johansen  KM.  2008. Chapter 1 chromatin structure and the regulation of gene expression: the lessons of PEV in Drosophila. In: Advances in genetics. Vol. 61 . Academic Press. p. 1–43. https://www.sciencedirect.com/science/article/pii/S0065266007000016?via%3Dihub.18282501
Graba  Y, Laurenti  P, Perrin  L, Aragnol  D, Pradel  J. The modifier of variegation modulo gene acts downstream of dorsoventral and HOM-C genes and is required for morphogenesis in Drosophila. Dev Biol. 1994:166 (2 ):704–715. 10.1006/dbio.1994.1349.7813788
Greil  F, de Wit  E, Bussemaker  HJ, van Steensel  B. HP1 controls genomic targeting of four novel heterochromatin proteins in Drosophila. EMBO J.  2007:26 (3 ):741–751. 10.1038/sj.emboj.7601527.17255947
Guthmann  M, Burton  A, Torres-Padilla  M. Expression and phase separation potential of heterochromatin proteins during early mouse development. EMBO Rep. 2019:20 (12 ):e47952. 10.15252/embr.201947952.31701657
Hahn  M, Dambacher  S, Schotta  G. Heterochromatin dysregulation in human diseases. J App Physiol. 2010:109 (1 ):232–242. 10.1152/japplphysiol.00053.2010.
Haig  D . Transposable elements: self-seekers of the germline, team-players of the soma. Bioessays. 2016:38 (11 ):1158–1166. 10.1002/bies.201600125.27604404
Hastings  P, Lupski  JR, Rosenberg  SM, Ira  G. Mechanisms of change in gene copy number. Nat Rev Genet. 2009:10 (8 ):551–564. 10.1038/nrg2593.19597530
Heitz  E . Das heterochromatin der moose. Jahrb Wiss Botanik. 1928:69 :762–818.
Helleu  Q, Levine  MT. Recurrent amplification of the heterochromatin protein 1 (HP1) gene family across Diptera. Mol Biol Evol. 2018:35 (10 ):2375–2389. 10.1093/molbev/msy128.29924345
Henikoff  S, Ahmad  K, Malik  HS. The centromere paradox: stable inheritance with rapidly evolving DNA. Science. 2001:293 (5532 ):1098–1102. 10.1126/science.1062939.11498581
Herz  H-M, Morgan  M, Gao  X, Jackson  J, Rickels  R, Swanson  SK, Florens  L, Washburn  MP, Eissenberg  JC, Shilatifard  A. Histone H3 lysine-to-methionine mutants as a paradigm to study chromatin signaling. Science. 2014:345 (6200 ):1065–1070. 10.1126/science.1255104.25170156
Hilfiker  A, Hilfiker-Kleiner  D, Pannuti  A, Lucchesi  JC. Mof, a putative acetyl transferase gene related to the Tip60 and MOZ human genes and to the SAS genes of yeast, is required for dosage compensation in Drosophila. EMBO J. 1997:16 (8 ):2054–2060. 10.1093/emboj/16.8.2054.9155031
Hoskins  RA, Carlson  JW, Kennedy  C, Acevedo  D, Evans-Holm  M, Frise  E, Wan  KH, Park  S, Mendez-Lago  M, Rossi  F, et al  Sequence finishing and mapping of Drosophila melanogaster heterochromatin. Science. 2007:316 (5831 ):1625–1628. 10.1126/science.1139816.17569867
Hoskins  RA, Carlson  JW, Wan  KH, Park  S, Mendez  I, Galle  SE, Booth  BW, Pfeiffer  BD, George  RA, Svirskas  R, et al  The release 6 reference sequence of the Drosophila melanogaster genome. Genome Res.  2015:25 (3 ):445–458. 10.1101/gr.185579.114.25589440
Hu  G, Katuwawala  A, Wang  K, Wu  Z, Ghadermarzi  S, Gao  J, Kurgan  L. flDPnn: accurate intrinsic disorder prediction with putative propensities of disorder functions. Nat Commun. 2021:12 (1 ):4438. 10.1038/s41467-021-24773-7.34290238
Huang  Y, Lee  YCG. Blessing or curse: how the epigenetic resolution of host-transposable element conflicts shapes their evolutionary dynamics. Proc Biol Sci  2024:291 (2020 ):20232775. 10.1098/rspb.2023.2775.38593848
Huang  Y, Shukla  H, Lee  YCG. Species-specific chromatin landscape determines how transposable elements shape genome evolution. eLife. 2022:11 :e81567. 10.7554/eLife.81567.35997258
Janssen  A, Colmenares  SU, Karpen  GH. Heterochromatin: guardian of the genome. Annu Rev Cell Dev Biol.  2018:34 (1 ):265–288. 10.1146/annurev-cellbio-100617-062653.30044650
Kamilar  JM, Cooper  N. Phylogenetic signal in primate behaviour, ecology and life history. Philos Trans R Soc Lond B Biol Sci. 2013:368 (1618 ):20120341. 10.1098/rstb.2012.0341.23569289
Kasinathan  B, Colmenares  SU  III, McConnell  H, Young  JM, Karpen  GH, Malik  HS. Innovation of heterochromatin functions drives rapid evolution of essential ZAD-ZNF genes in Drosophila. eLife. 2020:9 :e63368. 10.7554/eLife.63368.33169670
Kassis  JA, Kennison  JA, Tamkun  JW. Polycomb and Trithorax group genes in Drosophila. Genetics. 2017:206 (4 ):1699–1725. 10.1534/genetics.115.185116.28778878
Kelleher  ES . Protein–protein interactions shape genomic autoimmunity in the adaptively evolving rhino-deadlock-cutoff complex. Genome Biol Evol.  2021:13 (7 ):evab132. 10.1093/gbe/evab132.34115120
Kellner  WA, Ramos  E, Bortle  KV, Takenaka  N, Corces  VG. Genome-wide phosphoacetylation of histone H3 at Drosophila enhancers and promoters. Genome Res. 2012:22 (6 ):1081–1088. 10.1101/gr.136929.111.22508764
Kendek  A, Wensveen  MR, Janssen  A. The sound of silence: how silenced chromatin orchestrates the repair of double-strand breaks. Genes. 2021:12 (9 ):1415. 10.3390/genes12091415.34573397
Khan  T, Douglas  GM, Patel  P, Ba  N, Moses  AN, M  A. Polymorphism analysis reveals reduced negative selection and elevated rate of insertions and deletions in intrinsically disordered protein regions. Genome Biol Evol.  2015:7 (6 ):1815–1826. 10.1093/gbe/evv105.26047845
Kim  D, Paggi  JM, Park  C, Bennett  C, Salzberg  SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019:37 (8 ):907–915. 10.1038/s41587-019-0201-4.31375807
Kim  BY, Wang  JR, Miller  DE, Barmina  O, Delaney  E, Thompson  A, Comeault  AA, Peede  D, D’Agostino  ER, Pelaez  J, et al  Highly contiguous assemblies of 101 drosophilid genomes. eLife. 2021:10 :e66405. 10.7554/eLife.66405.34279216
King  TD, Leonard  CJ, Cooper  JC, Nguyen  S, Joyce  EF, Phadnis  N. Recurrent losses and rapid evolution of the condensin II complex in insects. Mol Biol Evol.  2019:36 (10 ):2195–2204. 10.1093/molbev/msz140.31270536
Klattenhoff  C, Xi  H, Li  C, Lee  S, Xu  J, Khurana  JS, Zhang  F, Schultz  N, Koppetsch  BS, Nowosielska  A, et al  The Drosophila HP1 homolog rhino is required for transposon silencing and piRNA production by dual-strand clusters. Cell. 2009:138 (6 ):1137–1149. 10.1016/j.cell.2009.07.014.19732946
Kokot  M, Długosz  M, Deorowicz  S. KMC 3: counting and manipulating k-mer statistics. Bioinformatics. 2017:33 (17 ):2759–2761. 10.1093/bioinformatics/btx304.28472236
Kovaka  S, Zimin  AV, Pertea  GM, Razaghi  R, Salzberg  SL, Pertea  M. Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol.  2019:20 (1 ):278. 10.1186/s13059-019-1910-1.31842956
Krause  SA, Overend  G, Dow  JAT, Leader  DP. FlyAtlas 2 in 2022: enhancements to the Drosophila melanogaster expression atlas. Nucleic Acids Res. 2022:50 (D1 ):D1010–D1015. 10.1093/nar/gkab971.34718735
Kuo  Y-M, Andrews  AJ. Quantitating the specificity and selectivity of Gcn5-mediated acetylation of histone H3. PLoS One. 2013:8 (2 ):e54896. 10.1371/journal.pone.0054896.23437046
Kuznetsov  D, Tegenfeldt  F, Manni  M, Seppey  M, Berkeley  M, Kriventseva  EV, Zdobnov  EM. OrthoDB v11: annotation of orthologs in the widest sampling of organismal diversity. Nucleic Acids Res.  2023:51 (D1 ):D445–D451. 10.1093/nar/gkac998.36350662
Lack  JB, Cardeno  CM, Crepeau  MW, Taylor  W, Corbett-Detig  RB, Stevens  KA, Langley  CH, Pool  JE. The Drosophila genome nexus: a population genomic resource of 623 Drosophila melanogaster genomes, including 197 from a single ancestral range population. Genetics. 2015:199 (4 ):1229–1241. 10.1534/genetics.115.174664.25631317
Larson  AG, Elnatan  D, Keenen  MM, Trnka  MJ, Johnston  JB, Burlingame  AL, Agard  DA, Redding  S, Narlikar  GJ. Liquid droplet formation by HP1α suggests a role for phase separation in heterochromatin. Nature. 2017:547 (7662 ):236–240. 10.1038/nature22822.28636604
Lartillot  N, Poujol  R. A phylogenetic model for investigating correlated evolution of substitution rates and continuous phenotypic characters. Mol Biol Evol. 2011:28 (1 ):729–744. 10.1093/molbev/msq244.20926596
Lee  YCG . The role of piRNA-mediated epigenetic silencing in the population dynamics of transposable elements in Drosophila melanogaster. PLoS Genet. 2015:11 (6 ):e1005269. 10.1371/journal.pgen.1005269.26042931
Lee  YCG, Karpen  GH. Pervasive epigenetic effects of Drosophila euchromatic transposable elements impact their evolution. eLife. 2017:6 :e25762. 10.7554/eLife.25762.28695823
Lee  J-H, Kim  EW, Croteau  DL, Bohr  VA. Heterochromatin: an epigenetic point of view in aging. Exp Mol Med. 2020:52 (9 ):1466–1474. 10.1038/s12276-020-00497-4.32887933
Lee  YCG, Leek  C, Levine  MT. Recurrent innovation at genes required for telomere integrity in Drosophila. Mol Biol Evol.  2017:34 (2 ):467–482. 10.1093/molbev/msw248.27836984
Lee  YCG, Ogiyama  Y, Martins  NMC, Beliveau  BJ, Acevedo  D, Wu  C-t, Cavalli  G, Karpen  GH. Pericentromeric heterochromatin is hierarchically organized and spatially contacts H3K9me2 islands in euchromatin. PLoS Genet.  2020:16 (3 ):e1008673. 10.1371/journal.pgen.1008673.32203508
Levine  MT, McCoy  C, Vermaak  D, Lee  YCG, Hiatt  MA, Matsen  FA, Malik  HS. Phylogenomic analysis reveals dynamic evolutionary history of the Drosophila heterochromatin protein 1 (HP1) gene family. PLoS Genet.  2012:8 (6 ):e1002729. 10.1371/journal.pgen.1002729.22737079
Lewis  SH, Salmela  H, Obbard  DJ. Duplication and diversification of dipteran argonaute genes, and the evolutionary divergence of piwi and aubergine. Genome Biol Evol. 2016:8 (3 ):507–518. 10.1093/gbe/evw018.26868596
Li  H . A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011:27 (21 ):2987–2993. 10.1093/bioinformatics/btr509.21903627
Li  H . Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018:34 (18 ):3094–3100. 10.1093/bioinformatics/bty191.29750242
Li  Y, Danzer  JR, Alvarez  P, Belmont  AS, Wallrath  LL. Effects of tethering HP1 to euchromatic regions of the Drosophila genome. Development. 2003:130 (9 ):1817–1824. 10.1242/dev.00405.12642487
Li  H, Durbin  R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009:25 (14 ):1754–1760. 10.1093/bioinformatics/btp324.19451168
Lin  Y, Currie  SL, Rosen  MK. Intrinsically disordered sequences enable modulation of protein phase separation through distributed tyrosine motifs. J Biol Chem. 2017:292 (46 ):19110–19120. 10.1074/jbc.M117.800466.28924037
Liu  J, Ali  M, Zhou  Q. Establishment and evolution of heterochromatin. Ann N Y Acad Sci. 2020:1476 (1 ):59–77. 10.1111/nyas.14303.32017156
Locke  J, Kotarski  MA, Tartof  KD. Dosage-dependent modifiers of position effect variegation in Drosophila and a mass action model that explains their effect. Genetics. 1988:120 (1 ):181–198. 10.1093/genetics/120.1.181.3146523
Malik  HS, Vermaak  D, Henikoff  S. Recurrent evolution of DNA-binding motifs in the Drosophila centromeric histone. Proc Natl Acad Sci U S A. 2002:99 (3 ):1449–1454. 10.1073/pnas.032664299.11805302
McDonald  JH, Kreitman  M. Adaptive protein evolution at the Adh locus in Drosophila. Nature. 1991:351 (6328 ):652–654. 10.1038/351652a0.1904993
Meignin  C, Davis  I. UAP56 RNA helicase is required for axis specification and cytoplasmic mRNA localization in Drosophila. Dev Biol.  2008:315 (1 ):89–98. 10.1016/j.ydbio.2007.12.004.18237727
Mérel  V, Gibert  P, Buch  I, Rodriguez Rada  V, Estoup  A, Gautier  M, Fablet  M, Boulesteix  M, Vieira  C. The worldwide invasion of Drosophila suzukii is accompanied by a large increase of transposable element load and a small number of putatively adaptive insertions. Mol Biol Evol.  2021:38 (10 ):4252–4267. 10.1093/molbev/msab155.34021759
Moesa  HA , Wakabayashi S, Nakai K, Patil A. Chemical composition is maintained in poorly conserved intrinsically disordered regions and suggests a means for their classification. Mol Biosyst. 2012:8 (12 ):3262–3273. 10.1039/c2mb25202c.23076520
Mohn  F, Sienski  G, Handler  D, Brennecke  J. The rhino-deadlock-cutoff complex licenses noncanonical transcription of dual-strand piRNA clusters in Drosophila. Cell. 2014:157 (6 ):1364–1379. 10.1016/j.cell.2014.04.031.24906153
Muller  HJ . Types of visible variations induced by X-rays in Drosophila. Journ. of Gen. 1930:22 (3 ):299–334. 10.1007/BF02984195.
Necci  M, Piovesan  D, Tosatto  SCE. Critical assessment of protein intrinsic disorder prediction. Nat Methods. 2021:18 (5 ):472–481. 10.1038/s41592-021-01117-3.33875885
Ninova  M, Chen  Y-CA, Godneeva  B, Rogers  AK, Luo  Y, Tóth  F, Aravin  K, Arvin  AA. Su(var)2-10 and the SUMO pathway link piRNA-guided target recognition to chromatin silencing. Mol Cell.  2020:77 (3 ):556–570.e6. 10.1016/j.molcel.2019.11.012.31901446
Nott  TJ, Petsalaki  E, Farber  P, Jervis  D, Fussner  E, Plochowietz  A, Craggs  TD, Bazett-Jones  DP, Pawson  T, Forman-Kay  JD, et al  Phase transition of a disordered nuage protein generates environmentally responsive membraneless organelles. Mol Cell. 2015:57 (5 ):936–947. 10.1016/j.molcel.2015.01.013.25747659
Obbard  DJ, Gordon  KHJ, Buck  AH, Jiggins  FM. The evolution of RNAi as a defence against viruses and transposable elements. Philos Trans R Soc Lond B Biol Sci.  2009:364 (1513 ):99–115. 10.1098/rstb.2008.0168.18926973
Obbard  DJ, Jiggins  FM, Halligan  DL, Little  TJ. Natural selection drives extremely rapid evolution in antiviral RNAi genes. Curr Biol.  2006:16 (6 ):580–585. 10.1016/j.cub.2006.01.065.16546082
Pak  CW, Kosno  M, Holehouse  AS, Padrick  SB, Mittal  A, Ali  R, Yunus  AA, Liu  DR, Pappu  RV, Rosen  MK. Sequence determinants of intracellular phase separation by complex coacervation of a disordered protein. Mol Cell. 2016:63 (1 ):72–85. 10.1016/j.molcel.2016.05.042.27392146
Parhad  SS, Tu  S, Weng  Z, Theurkauf  WE. Adaptive evolution leads to cross-species incompatibility in the piRNA transposon silencing machinery. Dev Cell.  2017:43 (1 ):60–70.e5. 10.1016/j.devcel.2017.08.012.28919205
Park  JI, Bell  GW, Yamashita  YM. Derepression of Y-linked multicopy protamine-like genes interferes with sperm nuclear compaction in D. melanogaster. Proc Natl Acad Sci U S A. 2023:120 (16 ):e2220576120. 10.1073/pnas.2220576120.37036962
Peacock  WJ, Lohe  AR, Gerlach  WL, Dunsmuir  P, Dennis  ES, Appels  R. Fine structure and evolution of DNA in heterochromatin. Cold Spring Harb Symp Quant Biol. 1978:42 (0 ):1121–1135. 10.1101/SQB.1978.042.01.113.98264
Regnard  C, Straub  T, Mitterweger  A, Dahlsveen  IK, Fabian  V, Becker  PB. Global analysis of the relationship between JIL-1 kinase and transcription. PLoS Genet.  2011:7 (3 ):e1001327. 10.1371/journal.pgen.1001327.21423663
Revell  LJ . Phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol.  2012:3 (2 ):217–223. 10.1111/j.2041-210X.2011.00169.x.
Rosin  L, Mellone  BG. Co-evolving CENP-A and CAL1 domains mediate centromeric CENP-A deposition across Drosophila species. Dev Cell. 2016:37 (2 ):136–147. 10.1016/j.devcel.2016.03.021.27093083
Ross  BD, Rosin  L, Thomae  AW, Hiatt  MA, Vermaak  D, de la Cruz  AFA, Imhof  A, Mellone  BG, Malik  HS. Stepwise evolution of essential centromere function in a Drosophila neogene. Science. 2013:340 (6137 ):1211–1214. 10.1126/science.1234393.23744945
Rudolph  T, Yonezawa  M, Lein  S, Heidrich  K, Kubicek  S, Schäfer  C, Phalke  S, Walther  M, Schmidt  A, Jenuwein  T, et al  Heterochromatin formation in Drosophila is initiated through active removal of H3K4 methylation by the LSD1 homolog SU(VAR)3-3. Mol Cell.  2007:26 (1 ):103–115. 10.1016/j.molcel.2007.02.025.17434130
Satyaki  PRV, Cuykendall  TN, Wei  KH-C, Brideau  NJ, Kwak  H, Aruna  S, Ferree  PM, Barbash  JS, A  D. The Hmr and Lhr hybrid incompatibility genes suppress a broad range of heterochromatic repeats. PLoS Genet.  2014:10 (3 ):e1004240. 10.1371/journal.pgen.1004240.24651406
Sawyer  IA, Sturgill  D, Dundr  M. Membraneless nuclear organelles and the search for phases within phases. Wiley Interdiscip Rev RNA.  2019:10 (2 ):e1514. 10.1002/wrna.1514.30362243
Schotta  G . Central role of Drosophila SU(VAR)3-9 in histone H3-K9 methylation and heterochromatic gene silencing. EMBO J.  2002:21 (5 ):1121–1131. 10.1093/emboj/21.5.1121.11867540
Schotta  G, Lachner  M, Sarma  K, Ebert  A, Sengupta  R, Reuter  G, Reinberg  D, Jenuwein  T. A silencing pathway to induce H3-K9 and H4-K20 trimethylation at constitutive heterochromatin. Genes Dev. 2004:18 (11 ):1251–1262. 10.1101/gad.300704.15145825
Shaffer  CD, Cenci  G, Thompson  B, Stephens  GE, Slawson  EE, Adu-Wusu  K, Gatti  M, Elgin  SCR. The large isoform of Drosophila melanogaster heterochromatin protein 2 plays a critical role in gene silencing and chromosome structure. Genetics. 2006:174 (3 ):1189–1204. 10.1534/genetics.106.057604.16980400
Sievers  F, Wilm  A, Dineen  D, Gibson  TJ, Karplus  K, Li  W, Lopez  R, McWilliam  H, Remmert  M, Söding  J, et al  Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol Syst Biol. 2011:7 (1 ):539. 10.1038/msb.2011.75.21988835
Smith  NGC, Eyre-Walker  A. Adaptive protein evolution in Drosophila. Nature. 2002:415 (6875 ):1022–1024. 10.1038/4151022a.11875568
Storer  J, Hubley  R, Rosen  J, Wheeler  TJ, Smit  AF. The Dfam community resource of transposable element families, sequence models, and genome annotations. Mob DNA. 2021:12 (1 ):2. 10.1186/s13100-020-00230-y.33436076
Strom  AR, Emelyanov  AV, Mir  M, Fyodorov  DV, Darzacq  X, Karpen  GH. Phase separation drives heterochromatin domain formation. Nature. 2017:547 (7662 ):241–245. 10.1038/nature22989.28636597
Suvorov  A, Kim  BY, Wang  J, Armstrong  EE, Peede  D, D’Agostino  ERR, Price  DK, Waddell  PJ, Lang  M, Courtier-Orgogozo  V, et al  Widespread introgression across a phylogeny of 155 Drosophila genomes. Curr Biol.  2022:32 (1 ):111–123.e5. 10.1016/j.cub.2021.10.052.34788634
Suyama  M, Torrents  D, Bork  P. PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucleic Acids Res.  2006:34 (Web Server ):W609–W612. 10.1093/nar/gkl315.16845082
Swanson  WJ, Nielsen  R, Yang  Q. Pervasive adaptive evolution in mammalian fertilization proteins. Mol Biol Evol.  2003:20 (1 ):18–20. 10.1093/oxfordjournals.molbev.a004233.12519901
Swenson  JM, Colmenares  SU, Strom  AR, Costes  SV, Karpen  GH. The composition and organization of Drosophila heterochromatin are heterogeneous and dynamic. eLife. 2016:5 :e16096. 10.7554/eLife.16096.27514026
Tatavosian  R, Kent  S, Brown  K, Yao  T, Duc  HN, Huynh  TN, Zhen  CY, Ma  B, Wang  H, Ren  X. Nuclear condensates of the Polycomb protein chromobox 2 (CBX2) assemble through phase separation. J Biol Chem. 2019:294 (5 ):1451–1463. 10.1074/jbc.RA118.006620.30514760
The UniProt Consortium . UniProt: the universal protein knowledgebase in 2023. Nucleic Acids Res.  2023:51 (D1 ):D523–D531. 10.1093/nar/gkac1052.36408920
Vermaak  D, Henikoff  S, Malik  HS. Positive selection drives the evolution of rhino, a member of the heterochromatin protein 1 family in Drosophila. PLoS Genet. 2005:1 (1 ):96–108. 10.1371/journal.pgen.0010009.16103923
Verschure  PJ, van der Kraan  I, de Leeuw  W, van der Vlag  J, Carpenter  AE, Belmont  AS, van Driel  R. In vivo HP1 targeting causes large-scale chromatin condensation and enhanced histone lysine methylation. Mol Cell Biol.  2005:25 (11 ):4552–4564. 10.1128/MCB.25.11.4552-4564.2005.15899859
Villeponteau  B . The heterochromatin loss model of aging. Exp Gerontol. 1997:32 (4-5 ):383–394. 10.1016/S0531-5565(96)00155-6.9315443
Wang  L-G, Lam  TT-Y, Xu  S, Dai  Z, Zhou  L, Feng  T, Guo  P, Dunn  CW, Jones  BR, Bradley  T, et al  Treeio: an R package for phylogenetic tree input and output with richly annotated and associated data. Mol Biol Evol.  2020:37 (2 ):599–603. 10.1093/molbev/msz240.31633786
Wang  Y, McNeil  P, Abdulazeez  R, Pascual  M, Johnston  SE, Keightley  PD, Obbard  D. Variation in mutation, recombination, and transposition rates in Drosophila melanogaster and Drosophila simulans. Genome Res. 2023:33 (4 ):587–598. 10.1101/gr.277383.122.37037625
Wei  KH-C, Grenier  JK, Barbash  DA, Clark  AG. Correlated variation and population differentiation in satellite DNA abundance among lines of Drosophila melanogaster. Proc Natl Acad Sci U S A.  2014:111 (52 ):18793–18798. 10.1073/pnas.1421951112.25512552
Wei  KH-C, Lower  SE, Caldas  IV, Sless  TJS, Barbash  DA, Clark  AG. Variable rates of simple satellite gains across the Drosophila phylogeny. Mol Biol Evol.  2018:35 (4 ):925–941. 10.1093/molbev/msy005.29361128
Wyman  MJ, Cutter  AD, Rowe  L. Gene duplication in the evolution of sexual dimorphism. Evolution. 2012:66 (5 ):1556–1566. 10.1111/j.1558-5646.2011.01525.x.22519790
Yang  Z . PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol.  2007:24 (8 ):1586–1591. 10.1093/molbev/msm088.17483113
Yang  Z, Wong  WSW, Nielsen  R. Bayes empirical Bayes inference of amino acid sites under positive selection. Mol Biol Evol. 2005:22 (4 ):1107–1118. 10.1093/molbev/msi097.15689528
Zarin  AA, Mark  B, Cardona  A, Litwin-Kumar  A, Doe  CQ. A multilayer circuit architecture for the generation of distinct locomotor behaviors in drosophila. eLife. 2019:8 (n.d. ):e51781. 10.7554/eLife.51781.31868582
Zarin  W, Veroniki  AA, Nincic  V, Vafaei  A, Reynen  E, Motiwala  SS, Antony  J, et al . Characteristics and knowledge synthesis approach for 456 network meta-analyses: a scoping review. BMC Med. 2017:15 :3 . 10.1186/s12916-016-0764-6.
Zhang  Y, Chu  J, Cheng  H, Li  H. De novo reconstruction of satellite repeat units from sequence data. Genome Res. 2023:33 (11 ):1994–2001. 10.1101/gr.278005.123.37918962
