GBE High Stability of the Epigenome in Drosophila Interspecific Hybrids n1, Marie Fablet2,3, Philippe Veber2, Cristina Vieira2, and Maria Pilar Garcıa Guerreiro1,* Alejandra Bodelo 1 mica, Bioinformatica i Biologia Evolutiva, Departament de Genètica i Microbiologia (Edifici C), Universitat Auto Grup de Geno noma de Barcelona, Spain 2 Laboratoire de Biometrie et Biologie Evolutive, UMR5558, Universite Claude Bernard Lyon 1, Villeurbanne, France 3 Institut universitaire de France, France *Corresponding author: E-mail: mariapilar.garcia.guerreiro@uab.es. Accepted: 2 February 2022 Abstract Interspecific hybridization is often seen as a genomic stress that may lead to new gene expression patterns and deregulation of transposable elements (TEs). The understanding of expression changes in hybrids compared with parental species is essential to disentangle their putative role in speciation processes. However, to date we ignore the detailed mechanisms involved in genomic deregulation in hybrids. We studied the ovarian transcriptome and epigenome of the Drosophila buzzatii and Drosophila koepferae species together with their F1 hybrid females. We found a trend toward underexpression of genes and TE families in hybrids. The epigenome in hybrids was highly similar to the parental epigenomes and showed intermediate histone enrichments between parental species in most cases. Differential gene expression in hybrids was often associated only with changes in H3K4me3 enrich- ments, whereas differential TE family expression in hybrids may be associated with changes in H3K4me3, H3K9me3, or H3K27me3 enrichments. We identified specific genes and TE families, which their differential expression in comparison with the parental species was explained by their differential chromatin mark combination enrichment. Finally, cis–trans compensatory regulation could also contribute in some way to the hybrid deregulation. This work provides the first study of histone content in Drosophila interspecific hybrids and their effect on gene and TE expression deregulation. Key words: epigenome, Drosophila, interspecific hybrids, deregulation, transposable elements, histone methylation. Significance The genomic stress caused by hybridization between different species is a cause of gene and transposable element (TE) deregulation. Histone modifications were often associated to genomic deregulation, but little is known to date about their role in Drosophila hybrid anomalies, and their elucidation is essential to unravel the role of these phenomena in speciation. Using transcriptomic and epigenomic analyses, we found that even though the epigenome of the parental species is highly conserved in hybrids, some changes detected in histone marks were associated to gene and TE deregulation. This work provides the first study of histone modification in Drosophila hybrids and its putative impact on gene a TE expression. Introduction Interspecific hybridizations are a source of phenotypic varia- tion (Soto et al. 2008; Zhu, Hu, et al. 2017), sterility (Naveira C The Author(s) 2022. Published by Oxford University Press on behalf of the Society for V and Fontdevila 1991; Michalak and Noor 2003; Moehring Molecular Biology and Evolution. et al. 2007), and genomic instability (Naveira and Fontdevila 1991; Vela et al. 2014, 2011) in hybrid offspring. These phe- This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unre- nomena are induced by the fusion of two different parental stricted reuse, distribution, and reproduction in any medium, provided the original work genomes and epigenomes, which act as a genomic stressful is properly cited. Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 1 n et al. Bodelo GBE factor in the hybrid genome (Garcıa Guerreiro 2012). This H3K27me3 are associated to gene repression (Yin et al. integration leads to deregulation of gene expression in 2011). H3K4me3 marks are abundant in euchromatin, and hybrids, which has been described in Drosophila (Michalak H3K9me3 and H3K27me3 take part of constitutive and fac- and Noor 2003; Ranz et al. 2004; Moehring et al. 2007; ultative heterochromatin, respectively (Ebert et al. 2006; Romero-Soriano et al. 2017; Gamez-Visairas et al. 2020) Schulze and Wallrath 2007). Moreover, it was also shown and other organisms (Malone and Michalak 2008; Zhu, Hu, that TEs are enriched in the H3K9me3 and H3K27me3 het- et al. 2017). For example, the misexpression of genes involved erochromatic marks, suggesting a putative colocalization of in spermatogenesis in Drosophila male hybrids (Michalak and these silencing marks in these sequences (Yin et al. 2011). In Noor 2003; Moehring et al. 2007) has been related with their plants, histone modification occurring during interspecific hy- sterility. Deregulation of gene expression was also observed in bridization has been associated to gene expression variation hybrid females of Drosophila buzzatii and Drosophila koep- (He et al. 2010; Zhu, Hu, et al. 2017). However, in Drosophila, ferae (Romero-Soriano et al. 2017; Gamez-Visairas et al. little is known about the role of histone modification on gene 2020) and in other Drosophila hybrids (Ranz et al. 2004). expression and their consequences on hybrid incompatibility Gene expression in Drosophila is controlled at different levels, and TE deregulation. Studies on D. melanogaster–D. simulans including epigenetic modifications (Yin et al. 2011) and cis- hybrids suggest that adaptive divergence of heterochromatin and trans-regulatory elements, among others. Divergence be- proteins is an important force driving the evolution of genes tween parental regulatory elements may play an important involved in hybrid incompatibility (Satyaki et al. 2014). Our role in gene expression deregulation in hybrids (Mack and previous work on D. buzzatii–D. koepferae hybrids suggested Nachman 2017). For example, differences in the cis and trans the existence of interacting phenomena, including incompat- parental regulatory elements were related with misexpression ibilities of the piRNA pathway, due to a functional divergence in Drosophila hybrids (Landry et al. 2005). between parental species, as one of the causes responsible of In addition to gene expression disruption, a growing num- TE mobilization (Romero-Soriano et al. 2017; review ber of evidence in Drosophila suggest that TEs could also play Fontdevila 2019). However, we cannot disregard other puta- an important role in the hybrid anomalies. For example, tive mechanisms in these species, including histone modifica- increases of transposition rates and/or misexpression were tions causing changes in TE expression. observed in interspecific hybrids between D. buzzatii and To unravel the causes of this hybrid incompatibility, we D. koepferae (Labrador et al. 1999; Vela et al. 2014; Garcıa investigated the expression profiles and the patterns of three Guerreiro 2015; Romero-Soriano and Garcia Guerreiro 2016) histone marks: H3K4me3, H3K9me3, and H3K27me3, using and other Drosophila species (Kelleher et al. 2012; Lopez- RNA-sequencing (RNAseq) and chromatin immunoprecipita- Maestre et al. 2017). In Drosophila, germline transposition is tion and deep sequencing (ChIPseq) in the genome of the controlled at two levels: transcriptional silencing, by hetero- parental species (D. buzzatii and D. koepferae) and their chromatinization processes, and post-transcriptional by piwi- hybrids. Drosophila buzzatii and D. koepferae are two cacto- interacting RNAs (piRNAs) involved in transposon transcript philic sibling species, belonging to repleta group (Garcıa degradation. In this sense, the causes of TE release in hybrids Guerreiro 2014), that diverged approximately 4–5 Ma have been the matter of recent research focused on the study (Go mez and Hasson 2003; Oliveira et al. 2012; Romero- of disruption of the germline mechanisms involved in TE reg- Soriano et al. 2017). Crosses between D. buzzatii males and ulation, basically those affecting piRNA production (Kelleher D. koepferae females result in F1 sterile males and fertile et al. 2012; Lopez-Maestre et al. 2017; Romero-Soriano et al. females (Marin et al. 1993), even though a few cases of partial 2017). Recent discoveries pointed to a dysfunction of piRNA fertility with atrophy in one of the ovaries is also observed pathway in Drosophila melanogaster–Drosophila simulans (Marin and Fontdevila 1998). Hybrid females can be back- hybrids (Kelleher et al. 2012), which would have a deficient crossed with D. buzzatii males (Naveira et al. 1986; Marin global piRNA production, whereas a similar or higher piRNA and Fontdevila 1998). In this system, we can only obtain production was observed in D. buzzatii-D. koeferae ovaries hybrids from crosses between D. koepferae females and (Romero-Soriano et al. 2017). Moreover, in the later study a D. buzzatii males, the reciprocal cross does not produce off- high proportion of the TEs overexpressed in hybrids did not spring (Marin et al. 1993). have associated piRNAs in parents or hybrids, which suggests We found that both genes and TE families detected as a more complex deregulation network, at least, in these differentially expressed in hybrid ovaries in comparison with hybrids. D. buzzatii and D. koepferae, tended to be mostly underex- Less is known about the behavior of epigenomes in inter- pressed. In contrast, we found a high conservation of the specific Drosophila hybrids, which is of fundamental interest parental chromatin mark patterns in hybrid genes and TEs, in evolutionary biology. Histone modification is a post- showing intermediate levels between the parental species. translational chemical modification of histone proteins. In Nevertheless, we could associate some changes in gene and Drosophila, the histone marks H3K4me3 are associated to TE family expression in hybrids with their corresponding his- activation of gene expression, whereas H3K9me3 and tone modifications versus parental species. 2 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE Results D. koepferae and D. buzzatii, respectively, whereas they de- creased to 1.62–1.72 (FC ¼ 0.33–0.30) for genes and to Deregulation of Gene and TE Expression in Hybrids 1.85–1.74 (FC ¼ 0.28–0.30) for TE families underexpressed We analyzed and compared the ovarian transcript amounts of in hybrids versus D. koepferae and D. buzzatii, respectively. D. buzzatii and D. koepferae and their F1 female hybrids. We Regarding the genes, most of the overexpressed shown in found that, out of a total of 13,621 protein coding genes and figure 1D, were involved in metabolism (small molecules or 658 TE families, 5.92% and 29.64%, respectively, were dif- protein), development, signaling and stimulus response, re- ferentially expressed between the parental species (supple- production and transportation, and cell organization (supple- mentary file 1: table S1, Supplementary Material online). mentary file 2, Supplementary Material online). Similarly, hybrids also showed a lower percentage of differ- An interesting example was the no hitter gene (nht), entially expressed genes, in comparison with the parental spe- highlighted in bold in figure 1D, which was overexpressed cies (4.57% vs. D. buzzatii and 3.99% vs. D. koepferae, in hybrids compared with both parental species. This gene is fig. 1A), than TE families (22.95% vs. D. buzzatii and involved in spermatogenesis and regulation of gene expres- 24.16% vs. D. koepferae, fig. 1B). Gene expression in hybrids sion and may be related to female hybrid sterility (Ranz et al. was more similar to the maternal species D. koepferae than to 2004). Moreover, as expected, most of the genes with higher the paternal species D. buzzatii (Z-test, P ¼ 0.018, supplemen- differences between hybrids and the parental species (gene tary file 1: table S2, Supplementary Material online), but the names labeled in fig. 1D) were at the same time detected as number of differentially expressed TE families in hybrids was the most different between the parental species (green and similar in comparison with both parental species (Z-test, blue gene names in supplementary file 3: fig. S1a, P ¼ 0.603, supplementary file 1: table S2, Supplementary Supplementary Material online). Material online). Regarding underexpressed genes, the Gene Ontology (GO) Globally, we could classify the differentially expressed analysis revealed a shared enrichment for 11 terms in hybrids genes and TEs in hybrids in three categories (fig. 1C): 1) versus both parental species (green GO in supplementary file D. koepferae- or D. buzzatii-like expression, which corre- 2, Supplementary Material online) and were mainly related to sponded to the majority of differentially expressed genes developmental processes, cell adhesion and reproduction (i.e., (fig. 1A) and TE families (fig. 1B). 2) Overexpressed or under- gonad development). Furthermore, the fact that most of expressed in hybrids in comparison with both parental spe- them had functions related with reproduction, reinforces cies, which were considered as deregulated in hybrids the idea of a putative role in hybrid fertility loss (Moehring (fig. 1C). We observed a deregulation trend toward under- et al. 2007). In the case of the overexpressed genes in hybrids, expression in hybrids compared with both parental species: only two enriched GO terms were shared with both parental 145 genes and 33 TE families underexpressed versus 29 genes species (green GO in supplementary file 2, Supplementary and 12 TE families overexpressed (fig. 1A and B). 3) Additive expression which included differentially expressed genes and Material online), and were linked to metabolic and cellular TE families in hybrids with an intermediate expression be- process. Additional GO terms, related to the same or other tween both parental species. Only a small number of genes processes (such as response to stimulus and biological regu- (29 and 40 genes) were included in this category (fig. 1A). lation, among others) were enriched in the differentially Regarding the TEs, the additive expressed TE families were expressed genes, but only in comparison with one of the pa- mostly underexpressed in comparison with D. buzzatii, and rental species (supplementary file 2, Supplementary Material overexpressed in comparison with D. koepferae (fig. 1B). Both online). results were consistent with the differences between the pa- To detect a putative location bias of the differentially rental species, where the mean TE family expression was expressed genes in hybrids in comparison with the parental higher in D. buzzatii than in D. koepferae (mean logarithmic species, previously reported in some Drosophila hybrids 2-fold change [log2FC] ¼ 0.27 or 1.21-fold increase in (Moehring et al. 2007; Wurmser et al. 2011), we studied their D. buzzatii vs. D. koepferae), and opposite for the genes distribution per chromosome. We found that the number of (mean log2FC¼ 0.23 or 0.85-fold decrease in D. buzzatii differentially expressed genes across chromosomes in hybrids vs. D. koepferae; supplementary file 1: table S1, compared with both parental species (supplementary file 3: Supplementary Material online). fig. S2a, Supplementary Material online), was different from We found that, even though there was a general bias to- the random expectation (chi-square test, P < 0.001, supple- ward underexpression of genes (fig. 1D) and TE families mentary file 1: table S3, Supplementary Material online). This (fig. 1E) in hybrids, those detected as overexpressed exhibit result is due to a higher number of differentially expressed the highest differences of expression compared with parental genes on the dot chromosome 6 (1.93% vs. D. buzzatii and species. In fact, the values of mean log2FC were 2.53–3.01 2.21% vs. D. koepferae) than the expected (0.54%). When (FC ¼ 5.78–8.06) for genes and 3.33–3.56 (FC ¼ 10.06– this chromosome was excluded from the test, the differences 11.79) for TE families overexpressed in hybrids versus become non-significant (chi-square test, D. buzzatii: Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 3 n et al. Bodelo GBE Fig. 1.—Differential expression analyses of genes and TEs: (A) Number of differentially expressed genes and (B) TE families in hybrids versus parental species of the total of 13,621 genes and 658 TE families. Colors indicate gene expression categories. (C) Expression categories in hybrids versus parental species. (D) Differential gene expression analyses in hybrids versus D. buzzatii (green) and D. koepferae (blue). Positive log2FC values correspond to genes more expressed in hybrids. The genes showing the 20 highest log2FC values and displaying an ortholog in D. melanogaster are shown. Genes common to both comparisons are in bold. (E) Violin plots representing the distribution of differentially expressed orders of TEs in hybrids versus parental species. Points indicate the log2FC of each TE family and for each comparison. The percentages of differentially expressed TE families per order are framed. The expected (total) TE family percentages are at the bottom. Red-dashed lines indicate the log2FC threshold for significancy (61). 4 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE P ¼ 0.320 and D. koepferae: P ¼ 0.443, supplementary file 1: model (RNA  K4 þ K9 þ K27 þ Input) described in Material table S3, Supplementary Material online). and Methods. We found that 62% (r 2 adjusted) of gene ex- Concerning the TEs, we studied the distribution and ex- pression variation was explained by the additive linear rela- pression of the differentially expressed TE families, per order, tionship with log-transformed histone mark enrichments, in hybrids in comparison with the parental species. We used with a P-value lower than 2.2  1016 in all species (supple- the Repbase classification (Jurka et al. 2005) where TEs were mentary file 1: table S4, Supplementary Material online). In divided in LTR, LINE, and DNA orders. Due to their particular the case of TE expression, the adjusted r 2 reached values from replication mechanism, we considered the RC/Helitron as a 61% in hybrids, to 75% in D. koepferae (supplementary file 1: group apart from the DNA order. As shown in figure 1E, the table S4, Supplementary Material online). number of differentially expressed TE families by order was Globally, we detected that the size effect, representing the similar to the random expectation considering their relative contribution of each chromatin mark to the gene expression proportions in the genome (chi-square test, P ¼ 0.183 for according to our model, was similar in hybrids and parental D. buzzatii and D. koepferae, supplementary file 1: table S3, species (fig. 2C), being H3K27me3 the chromatin mark with Supplementary Material online). The most extreme underex- the greatest contribution to gene expression. H3K4me3 was pression values were observed in LINE and LTR families in strongly positively associated with gene expression showing hybrids versus D. koepferae. In contrast, these families to- coefficient values from 1.59 in D. koepferae to 1.75 in gether with RC were the most overexpressed in hybrids versus hybrids, whereas H3K27me3 was strongly negatively associ- D. buzzatii. Last, the LTR order was the one where some ated: from 2.24 in D. koepferae, to 2.56 in hybrids (sup- families showed extreme overexpression and underexpression plementary file 1: table S4, Supplementary Material online). values in hybrids versus both parental species. When we went H3K9me3 was the one with the least positive contribution to deeper into the superfamilies of elements belonging to each gene expression (from 0.35 in D. koepferae to 0.5 in order (supplementary file 3: fig. S2b, Supplementary Material D. buzzatii). online), highly expressed elements in hybrids included mainly Regarding TEs, the contribution of H3K4me3 and Gypsy (LTR), Helitron (RC), and Jockey (LINE) elements in com- H3K9me3 to the expression was similar in hybrids and parental parison with D. buzzatii and, Gypsy (LTR), hAT (DNA), and species (fig. 2D). H3K4me3 was positively associated with TE TC1Mariner (DNA) elements compared with D. koepferae. expression, with coefficient values from 0.96 in D. koepferae Some of these TEs were detected as the most differentially to 1.58 in D. buzzatii (supplementary file 1: table S4, expressed between the parental species (green and blue TE Supplementary Material online), whereas H3K9me3 was neg- names in supplementary file 3: fig. S1b, Supplementary atively associated (1.04 in D. buzzatii to 1.16 in hybrids), as Material online). Consistently with what was found in previ- expected. Contrary to the pattern observed in genes, ous works (Romero-Soriano et al. 2017), the most highly over- H3K27me3 was the chromatin mark with the lowest contri- expressed TEs in hybrids versus both parental species bution to TE expression and its contribution depended on the belonged to Gypsy subfamily (supplementary file 3: fig. S2b, species: positively associated in D. buzzatii and hybrids (0.38 Supplementary Material online). and 0.71 respectively) and negatively in D. kopferae (0.33). Epigenetic Mark Landscapes Are Conserved and Maintenance of the Parental Enrichment of Epigenetic Contribute to Gene and TE Expression Marks in the Hybrid Genome In Drosophila, the constitutive heterochromatin is enriched in We computed read counts corresponding to enrichments in H3K9me3 epigenetic mark, whereas H3K27me3 and H3K4me3, H3K9me3, and H3K27me3, in the gene bodies H3K4me3 are associated to facultative heterochromatin and and TE families from hybrids and their parental species. We euchromatin, respectively (Boros 2012). We investigated the found that only a small percentage of all genes (1.40–2.83%) landscape of epigenetic marks in the parental species, un- were differentially enriched in hybrids compared with any pa- known to date, and in their hybrids. Analyses of the distribu- rental species and chromatin mark (fig. 3A). A lower number tion of histone marks in genes and their surrounding regions of differentially H3K4me3 and H3K27me3 enriched genes (fig. 2A and B) showed that H3K4me3 euchromatic mark is was detected in hybrids when they were compared with enriched around the start codon (SC) and throughout the D. buzzattii than to D. koepferae (Z-test, P < 0.001, supple- coding sequence in hybrids and their parental species, which mentary file 1: table S5, Supplementary Material online). is similar to the well-described D. melanogaster species. Regarding TE families, the hybrid versus parents comparisons H3K9me3 and H3K27me3 epigenetic marks were mainly de- revealed very similar patterns for H3K9me3 enrichments pleted in the gene body in all species (fig. 2A and B). (1.52–1.37% of differentially enriched TE families), whereas We then studied the contribution of these histone marks to they were the most contrasted for H3K27me3 enrichments the gene and TE family expression in the parental species (9.27–6.53% differentially enriched TE families, fig. 3B). D. buzzatii, D. koepferae, and their hybrids, using the linear However, the number of differentially enriched TEs in hybrids Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 5 n et al. Bodelo GBE Fig. 2.—Global epigenetic patterns in D. buzzatii, D.koepferae, and hybrids and their contribution to expression: (A) Average levels of the histone marks H3K4me3, H3K9me3, and H3K27me3 over all annotated gene sequences (exons and introns) and the 62 kb surrounding regions. (B) Heatmaps showing the density scores of each histone mark in D. buzzatii, D. koepferae, and their hybrids. White regions in the heatmaps indicate missing data. SC: start codon, EC: end codon. (C, D) Mean linear effects of H3K4me3, H3K9me3, and H3K27me3 enrichments on RNA read counts (log transformed) in (C) genes and (D) TE families. Colors indicate the species. The 95% confidence intervals are indicated. was similar to one or another parental species (Z-test, Supplementary Material online). On the contrary, in TE fam- P ¼ 0.817 H3K4me3 and H3K9me3; P ¼ 0.198 H3K27me3; ilies H3K27me3 was the mark with the largest differential supplementary file 1: table S5, Supplementary Material on- enrichment in hybrids versus the parental species, showing line). Moreover, H3K9me3 was the chromatin mark with the largest contrast between the parental species (supplemen- more similar patterns between hybrids and the parental spe- tary file 1: table S6, Supplementary Material online). cies for both genes and TE families, which was consistent with In general, most of the differentially enriched genes and the similar patterns of this chromatin mark when the parental TEs in hybrids have D. buzzatii or D. koepferae-like chromatin species were compared (supplementary file 1: table S6, mark levels. However, we also observe a high number of 6 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE Fig. 3.—Comparison of epigenetic marks in hybrids and parental species: (A) Number of differentially enriched genes and (B) TEs of the total of 13,621 genes and 658 TE families in H3K4me3, H3K9me3, or H3K27me3 chromatin marks in hybrids versus the parental species D. buzzatii and D. koepferae, respectively. Colors indicate the different categories of enrichment. (C) Violin plots representing the distribution of differentially H3K4me3, H3K9me3, and H3K27me3 enriched TE orders in hybrids versus parental species. Points indicate the size effect. The percentages of differentially enriched TE families per order are framed. The expected (total) TE family percentages are at the bottom. Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 7 n et al. Bodelo GBE genes displaying histone mark enrichments outside of the P  0.001: comparison with both parental species, supple- range of parental values, especially for H3K4me3—48 less mentary file 1: table S8, Supplementary Material online). enriched and 62 more enriched in hybrids versus both paren- Indeed, genes that were underexpressed in hybrids compared tal species (fig. 3A). In the same way, a high number of TE with D. buzzatii also showed a reduced H3K4me3 enrich- families show intermediate patterns of enrichment between ment. However, and unexpectedly, genes that were under- the parental species (fig. 3B). expressed in hybrids compared with D. koepferae also We studied the distribution across chromosomes of the displayed an increased H3K4me3 enrichment. This apparent differentially enriched genes for each epigenetic mark in inconsistency could be explained by the differences in enrich- hybrids in comparison with each parental species and we ment of this epigenetic mark between parental species: the did not find significant differences in any case (chi-square genes enriched in H3K4me3 in hybrids versus D. koepferae test, P > 0.05 in all cases, supplementary file 1: table S7, often corresponded to those more enriched in D. buzzatii Supplementary Material online). Globally, we observed that than in D. kopferae (blue color, fig. 4A), the opposite was the extreme changes in chromatin mark levels were toward also observed. an increase in hybrids in comparison with the parental species TE families, as genes, showed an intermediate inheritance (supplementary file 3: fig. S3, Supplementary Material online). of the chromatin marks in hybrids: TEs more enriched in a For example, some extreme changes of the euchromatic mark chromatin mark in hybrids in comparison with one parental H3K4me3 were observed in chromosome 5 in comparison species tended to be less enriched when they were compared with D. koepferae, and of the euchromatic mark with the other parental species (fig. 4G, I and K). In addition, H3K27me3 in comparison with D. koepferae (chromosomes we found an association between TE expression changes in 4 and 5) and D. buzzatii (chromosome 2). hybrids in comparison with the parental species and the cor- When we studied the orders of the differentially enriched responding changes in the epigenome. For example, an im- TE families in hybrids versus parental species, we did not find poverishment of the heterochromatic mark H3K9me3 was differences (chi-square test, P > 0.05 in all cases, supplemen- found in many underexpressed TEs in hybrids compared tary file 1: table S7, Supplementary Material online) for any with both parental species (Fisher’s exact test, P < 0.001, sup- chromatin mark (fig. 3C). If we focus on the specific enrich- plementary file 1: table S8, Supplementary Material online). In ment of epigenetic marks, three TEs belonging to DNA, LINE, the same way, a high TE number showing a decrease of and LTR families showed the most extreme values of H3K4me3 or H3K27me3 chromatin marks were underex- H3K4me3 in hybrids versus D. koepferae (ranging from 6 to presed in hybrids in comparison with D. buzzatii (Fisher’s exact 9). However, RC and unknown TEs showed values of this test, H3K4me3: P < 0.001 and H3K27me3: P ¼ 0.002, sup- chromatin mark similar to both parental species, being only plementary file 1: table S8, Supplementary Material online). a small percentage differentially enriched. We also found a These results were opposite to those found in hybrids com- high increment of H3K27me3 in one LTR element in compar- pared with D. koepferae, where an increase of H3K27me3 ison with D. koepferae (6), whereas the changes were less was observed in underexpressed elements (P < 0.001, supple- extreme in comparison with D. buzzatii. As reported for gene mentary file 1: table S8, Supplementary Material online). To enrichment, only small differences of H3K9me3 amounts determine if there was an association between chromatin were observed in hybrids versus parental species. marks in TEs, we checked the correlation between the eu- chromatic mark H3K4me3 and the other two, in hybrids in comparison with the parental species (supplementary file 3: Small Changes in the Hybrid Epigenome Affect Their Gene fig. S4, Supplementary Material online). We found that the and TE Expression increase of H3K4me3 in hybrids was associated with an in- We analyzed the association between gene and TE expression crease of H3K9me3 and H3K27me3 (Linear model, P  1.669 changes in hybrids in comparison with the parental species  1015 in all cases, supplementary file 1: table S9, and the corresponding changes in chromatin marks (fig. 4A– Supplementary Material online) and the opposite, with L, and supplementary file 1: table S8, Supplementary Material some TEs highly expressed and extremely enriched in online) using the whole set of genes and TE families (signifi- H3K4me3 in hybrids (red dots in the upper part right of sup- cant and nonsignificant). Regarding the genes, we found that plementary file 3: fig. S4, Supplementary Material online). changes in the enrichment of the heterochromatic marks When we focus on the genes considered as differentially H3K9me3 and H3K27me3 were not associated with expres- expressed and differentially enriched in hybrids, we observe sion changes in hybrids in comparison with the parental spe- that the epigenetic mark associated with the most extreme cies (Fisher’s exact test, D. buzzatii: P ¼ 0.763 both chromatin expression changes in hybrids (absolute log2FC values ranging marks; D. koepferae: P ¼ 0.970 and P ¼ 0.155, respectively, from 4 to 6) in comparison with both parental species was supplementary file 1: table S8, Supplementary Material on- H3K4me3 (fig. 4M). Changes in the remaining chromatin line). However, changes in H3K4me3 in hybrids seemed to be marks were related with less extreme expression changes (ab- associated with gene expression changes (Fisher’s exact test, solute log2FC values from 1 to 3). Regarding the TEs showing 8 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE Fig. 4.—Association of expression changes with the corresponding chromatin mark changes: (A–L) Differential expression (log2FC) and differential enrichment (size effect) of H3K4me3, H3K9me3, and H3K27me3 in all genes (B, D, F) and TE families (H, J, L) (significant and nonsignificant) in hybrids versus D. buzzatii and D. koepferae, respectively. Top plots represent the distribution of the differential enrichment between parental species in genes (A, C, E) and TE families (G, I, K). Colors represent the results of the differential enrichment between parental species, and values that cannot be computed are included in the NA category. Fisher’s exact test P-values are shown in the right corner. Gene outliers are not included in (A–F). (M, N) Log2FC and size effect representation of significant differentially expressed and enriched genes (M) and TEs (N) in hybrids versus the parental species. Significant genes (M) in at least two comparisons are named using their D. melanogaster ortholog or, otherwise, a number. The TE families (N) common to at least two comparisons are marked, indicating their category and a subscript. Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 9 n et al. Bodelo GBE significant differences in enrichment and expression in hybrids studied how the regulatory divergence influences expression versus parental species (fig. 4N), we observed two TE families differences between parental species. We found more differ- whose increase in H3K4me3 was associated to an extreme entially expressed genes between parental species with trans- overexpression in hybrids in comparison with D. buzzatii regulatory divergence (Fisher’s exact test, P < 0.01, supple- (log2FC from 6 to 11). Additionally, increases in the number mentary file 1: table S11, Supplementary Material online) of H3K27me3 marks in hybrids in comparison with than expected, whereas no differences to what was expected D. kopferae were likely related to extreme overexpression in were observed in genes with cis-regulatory divergence hybrids (log2FC from 3 to 6). Chromatin mark changes in the (Fisher’s exact test, P ¼ 0.128 for cis-regulatory, supplemen- remaining TEs were related with less extreme expression tary file 1: table S11, Supplementary Material online). Finally, changes (absolute log2FC values from 1 to 3). we studied how the regulatory divergence observed in If we focus on the expression and chromatin enrichment of parents influenced the inheritance of gene expression in specific genes (fig. 4M) in hybrids, we observe some examples hybrids. We analyzed genes with regulatory divergence be- whose expression is explained by the combination of different tween D. koepferae and D. buzzatii, which were differentially chromatin marks. For example, the CG3902 gene, involved in expressed in hybrids versus parental species, using the expres- oxidation–reduction processes, was underexpressed in hybrids sion categories described previously (fig. 1C). We found a in comparison with D. buzzatii and less enriched in H3K4me3 higher number of differentially expressed genes with trans- but more in H3K9me3. The gene ptc, involved in different regulatory divergence than with other regulatory classes in all processes including gonad development, was underexpressed expression categories, except those of genes deregulated in in hybrids in comparison with D. koepferae and more hybrids. In this latter category, we also observed a significantly enriched in H3K4me3 and H3K27me3. In other cases, the higher number of genes with compensatory regulation (cis- expression changes in genes (fig. 4M) and TEs (fig. 4N) could and trans-regulatory differences compensate each other) than be explained by the content of a unique epigenetic mark in the other expression categories (additive or D. buzzatii/ when hybrids are compared with one or another parental koepferae-like) (Fisher’s exact test, additive expression: species. For example, the gene Nepl4 (fig. 4M), involved in P ¼ 0.039, D. buzzatii and koepferae-like expression: protein processing, had a low H3K4me3 content in hybrids P < 0.001, supplementary file 1: table S12, Supplementary versus both parental species and was underexpressed in Material online). hybrids. On the other hand, the gene named 3 (fig. 4M) To gain insight into other mechanisms affecting TE dereg- and the unknown1 TE (fig. 4N) were detected as underex- ulation in hybrids, we examined the expression of piRNA path- pressed and less enriched in H3K4me3 in hybrids in compar- way genes considering their role in germline development, ison with D. buzzatii, but overexpressed and more enriched in epigenetic regulation, and TE silencing. We observed that au- this chromatin mark in comparison with D. koepferae. Finally, bergine (aub) and Sister of yellow body (SoYb) were deregu- LTR2 and DNA3 were enriched in H3K27me3 in comparison lated toward underexpression in hybrids in comparison with with one parental species, but the opposite in comparison both parental species (supplementary file 3: fig. S5, with the other (fig. 4N). Supplementary Material online). The chromatin mark levels were similar to the parental species, with only a few genes differentially enriched, mostly in comparison with only one Additional Mechanisms Influencing Hybrid Deregulation parental species. Globally, differentially enriched genes de- The divergence in cis- and trans- regulatory elements was creased H3K4me3 levels in hybrids in comparison with the proposed as one of the causes of expression differences be- parental species. Yellow body (Yb), which was detected as tween Drosophila species, as well as of gene deregulation in overexpressed in hybrids versus both parental species, was their hybrids (Wittkopp et al. 2004; Landry et al. 2005; intermediate enriched in H3K27me3 in hybrids versus both McManus et al. 2010). To detect regulatory divergence be- parental species, whereas vreteno (vret), which had an addi- tween D. buzzatii and D. koepferae and its putative effects on tive expression in hybrids versus parental species, was less hybrid expression, we compared the relative allele expression enriched in H3K4me3 in hybrids versus both parental species. in hybrids (H), the differential expression between parental Brother of Yellow body (BoYb) and armitage (armi), which species (P), and the ratio between these two metrics (T). We were detected as overexpressed in hybrids in comparison categorized the genes in different classes, as in (McManus with one parental species, were more enriched in H3K4me3 et al. 2010) and described in Material and Methods. As shown and less enriched in H3K9me3, respectively, in comparison in figure 5A, most of the genes (57.65%), categorized as with the same parental species. However, most differentially conserved, do not have regulatory divergence between pa- enriched genes were not detected as differentially expressed rental species. However, we found slightly more genes show- (supplementary file 3: fig. S5, Supplementary Material online). ing evidence of trans-regulatory divergence (6.83%) than cis- We next studied the association of piRNA amounts and TE regulatory divergence (6.24%) (Z-test, P ¼ 0.049, supplemen- expression in hybrids and parental species. We observed that, tary file 1: table S10, Supplementary Material online). We also on average, parental species showed differences in piRNA 10 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE Fig. 5.—Additional factors affecting gene and TE expression, respectively: (A) cis- and trans- regulatory divergence between parental species. The pie chart represents the percentage of total genes categorized by regulatory divergence class. (B) Differential expression (log2FC) values of TEs and piRNA amounts (log2FC) in hybrids versus D. buzzatii and D. koepferae, respectively. Colors represent the results of the piRNA amount changes between parental species. Values that cannot be computed are included in the NA category. Fisher’s exact test P-values are shown in the right corner. amounts: TE families were associated with more piRNAs in families were overexpressed in comparison with both parental D. koepferae than in D. buzzatii (mean log2FC of 2.14 in species (no data provided about underexpression). In contrast, D. buzzatii vs. D. koepferae comparison, fig. 5B), as observed 78% of differentially expressed genes of D. melanogaster– in our previous work (Romero-Soriano et al. 2017). Hybrids D. simulans female hybrid versus their parental species were exhibited an additive pattern of piRNA amounts between pa- found in other studies (Ranz et al. 2004). These discrepancies rental species: less piRNA amounts than D. koepferae and could be due to the different methodological approaches more than D. buzzatii. Additionally, as shown in figure 5B, used along with the different genetic background of the the differences in piRNA amounts were not associated to the Drosophila stocks. changes in TE transcript amounts in hybrids in comparison On the other hand, in our hybrid females most of the dif- with the parental species (Fisher’s exact test, P ¼ 0.328 for ferentially expressed genes tended to be underexpressed, as D. buzzatii and P ¼ 0.086 for D. koepferae, supplementary observed in previous studies in plants (Wang et al. 2006) and file 1: table S13, Supplementary Material online). There in Drosophila females (Ranz et al. 2004). The underexpression were also no differences when TE classes were analyzed sep- seems also to be the rule in Drosophila hybrid males between arately (Fisher’s exact test, P ¼ 0.732 for Class I and P ¼ 0.564 species of the melanogaster group (Michalak and Noor 2003; for Class II, supplementary file 1: table S13, Supplementary Moehring et al. 2007), which have been associated to male Material online). sterility. In our case, even if most F1 hybrid females are fertile, some are partially fertile (Marin and Fontdevila 1998), which Discussion could explain that GO terms of our underexpressed genes are associated to developmental processes, cell adhesion, and re- Bias to Underexpression of Genes and TEs in Hybrid production. Instead, only a few genes were extremely over- Females expressed in hybrids, for example the nht, which is involved in The analysis of the total ovarian transcriptome of parental spermatogenesis processes. The overexpression of male- species showed that 5.92% and 29.64% of the protein- specific genes in female hybrids is not new, and has been coding genes and TE families, respectively, were differentially attributed to a failure in the mechanisms controlling the ex- expressed between D. buzzatii and D. koepferae parental spe- pression of these genes in females (Ranz et al. 2004). cies. These differences are similar to those obtained in other If we focus on the resemblance of gene expression in studies with other species of the repleta group (Lopez- hybrids versus parental species, we found a bias toward genes Maestre et al. 2017). When the transcriptome of F1 hybrids that resemble more to one of the parental species, being the and their parental species were compared, around 4% of number of genes having an additive expression between both genes and 23% of TEs from hybrids were differentially parental species low, as observed in previous studies with expressed compared with any parental species. A greater de- D. melanogaster–D. simulans hybrids (Ranz et al. 2004). regulation of TEs than of genes was already observed in a Moreover, the number of genes in hybrids sharing similar previous work with D. melanogaster–D. simulans artificial transcript amounts with D. koepferae was higher than with hybrids (Kelleher et al. 2012), but the percentages were dif- D. buzzatii. Maternal effects were pointed out in Arabidopsis ferent: 0.7% of genes were deregulated and 12% of TE lyrata (Videvall et al. 2016), Xenopus (Malone and Michalak Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 11 n et al. Bodelo GBE 2008), and in coral (Chan et al. 2021) intraspecific hybrids. The Epigenomic Landscapes of Parental Species and Because only one direction of cross can be performed in our Hybrids Are Similar to Other Drosophila Species case, it is difficult to attribute these results unequivocally to Gene and TE deregulation in D. buzzatii and koepferae maternal effects. hybrids were previously described in the literature (Soto When the distribution of the derepressed genes in F1 hy- et al. 2008; Vela et al. 2014; Garcıa Guerreiro 2015; brid females along chromosomes was considered, we did not Romero-Soriano and Garcia Guerreiro 2016; Romero- find an overrepresentation of deregulated genes in any chro- Soriano et al. 2017; Gamez-Visairas et al. 2020), as well as mosome, except in chromosome 6 (dot). This bias could be in other hybrids of Drosophila (Moehring et al. 2007; Ortız- explained if this chromosome was different between our pa- Barrientos et al. 2007; Satyaki et al. 2014) and of other organ- rental species, which is reinforced by the highest rate of mo- isms (De Araujo et al. 2005; He et al. 2013; Zhu, Greaves, lecular evolution found in this dot chromosome in the closely et al. 2017). To get insight about these genomic deregulation related species D. mojavensis (Allan and Matzkin 2019). No mechanisms, we studied the epigenomic landscape of a eu- bias was found in a previous Drosophila study (Michalak and chromatic mark (H3K4me3) and two heterochromatic marks Noor 2003). However, an X-chromosome bias of differentially (H3K9me3 and H3K27me3) in D. buzzatii, D. koepferae, and expressed genes was found in other studies performed in their F1 hybrids, and we found a high similarity across these Drosophila hybrid males (Moehring et al. 2007; Wurmser species. Our results do not globally differ from those of the et al. 2011). The faster evolution of X-linked genes (Coyne D. melanogaster species epigenome described in the literature and Orr 1989) proposed to explain the differential expression and are consistent with the reported high conservation of the between X and autosomal genes has, however, been ques- active chromatin epigenome landscapes across Drosophila tioned by other authors (Hollocher and Wu 1996) who con- species (Brown and Bachtrog 2014). As reported in other sider that the homozygous autosome effects in reproductive studies (modENCODE Consortium et al. 2011; Yin et al. isolation, approach those of X chromosome. 2011; Boros 2012; Parey et al. 2019), we found the We showed that 6.84% of TE families were completely H3K4me3 active chromatin mark located at the 50 ends of deregulated in F1 hybrid females, in comparison with both actively transcribed genes and a depletion of H3K9me3 and parental species, and had a trend toward underexpression. H3K27me3 in the gene body. H3K9me3 was previously de- These results, obtained using a different approach (normali- scribed to be enriched in promoters but depleted in 50 tran- zation of TE counts using gene counts) and updated analysis scribed regions of active genes (Yin et al. 2011; Boros 2012), tools, contrast with our previous results where the number of whereas the H3K27me3 was underrepresented in the gene TE families upregulated slightly exceeded that of downregu- body (Filion et al. 2010; Yin et al. 2011; Boros 2012). lated (Romero-Soriano et al. 2017). However, a few TE fam- However, both were reported to be enriched in transposons ilies belonging to the Gypsy superfamily showed values of and repetitive sequences in Drosophila (Yin et al. 2011; Parey expression very high compared with parental species, concor- et al. 2019) and other organisms (Walter et al. 2016; dantly with previous results in these species where an increase Montgomery et al. 2020). of transposition and expression of the Osvaldo element was We next investigated whether the genomic expression was shown (Labrador et al. 1999; Vela et al. 2014; Garcıa explained by the epigenetic marks, finding that globally the Guerreiro 2015). Results on TE expression in hybrids reported gene expression was positively correlated with the active mark in the literature are very heterogeneous, finding cases of over- and was negatively correlated with the H3K27me3 repressive expression of specific elements by RT-qPCR (Carnelossi et al. mark, concordantly with previous studies in Drosophila (Yin 2014; Garcıa Guerreiro 2015) or by RNAseq (Kelleher et al. et al. 2011) and plants (He et al. 2010). Correlations between 2012). However, examples of underexpression in hybrids, af- the epigenomic landscape and gene expression in Drosophila fecting most (Renaut et al. 2014) or a high percentage of TEs have been previously described in other studies (Schübeler were also reported (Dion-Côte et al. 2014), but usually results et al. 2004; Yin et al. 2011). According to our model, 62% of underexpression, if they exist, are poorly discussed. Finally, of gene expression was explained by epigenetic marks, being no evidence of TE reactivation was found in natural hybrid H3K27me3 the one that most influence gene expression, lineages of Saccharomyces, suggesting that other factors like followed by H3K4me3 as observed in other Drosophila studies population structure and hybrid genotype are major determi- (Yin et al. 2011). H3K4me3 was positively associated with nants of TE content (Henault et al. 2020). The finding of gene expression, consistent with its enrichment in actively underexpressed TE families in our hybrids highlights that reg- transcribed genes found in other works (modENCODE ulation of some TEs exits in a way. Because cases of overex- Consortium et al. 2011; Yin et al. 2011; Boros 2012; Parey pressed TE families were also observed, we suggest that et al. 2019). H3K27me3 was found to be negatively corre- deregulation processes could be closely related to the TE fam- lated to mRNA levels, consistent with its reported function of ily and the genetic background of species involved in the hy- binding target for Polycomb repressive complex 1 (PRC1) bridization processes. (Riddle et al. 2011; Yin et al. 2011; Boros 2012; Parey et al. 12 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE 2019). Surprisingly, the epigenetic mark H3K9me3, usually Because in our study, we cannot distinguish individual TE cop- associated to heterochromatin regions (Yin et al. 2011), was ies, they were analyzed at a family level (658 TE families), positively associated with the gene expression, even if its con- meaning that changes in expression or chromatin mark tribution is very low. Increases of this epigenetic mark, along amounts are the result of the addition of the different copies with other euchromatic ones, have been described in intronic per family. Indeed, intraspecific and interspecific variations in and 30 end regions of some active heterochromatic genes in histone marks were observed in Drosophila TE copies (Rebollo some stages of D. melanogaster (Riddle et al. 2011; Saha et al. et al. 2012), indicating that specific epigenetic modifications 2019). in TE individual copies in our hybrids could go unnoticed. The epigenetic marks found on TEs explained more than Further inspection showed that increases of H3K4me3 were 60% of their expression, being H3K9me3 and H3K4me3 the associated with increases of H3K9me3 and H3K27me3 in most influential, with a negative and positive association to TE hybrids. Changes in the epigenetic status in individuals sub- expression, respectively. H3K9me3 together with H3K27me3, mitted to other stress were already observed in Bari-Jheh TE: a are known to be abundant in TEs and are associated to their H3K9me3 dominant pattern turned to increases in H3K4me3, silencing (Walter et al. 2016; Parey et al. 2019) and H3K4me3 H3K9me3, and H3K27me3 enrichments (Guio et al. 2018). to their activation in Drosophila (Riddle et al. 2011; Yin et al. 2011) and other organisms (Walter et al. 2016). The different Other Mechanisms May Have a Role in Hybrid Genome association of H3K27me3 with some species (positively cor- Deregulation related with TE expression in D. buzzatii and hybrids) could be explained by the different chromatin marks associated to each Even though in hybrids we found a relationship between epi- TE copy inside the same family. genetic status and expression changes in genes and TEs, the histone modification patterns did not account for the whole Hybrids Exhibited Limited Changes in Histone Marks genome deregulation. To get more insight into this aspect, we studied the regulatory divergence between D. buzzatii and We studied differences in the enrichment of three chromatin D. koepferae, and its contribution to gene expression differ- marks (H3K4me3, H3K9me3, and H3K27me3) in genes and ences between them and their hybrids. In general, we found TEs between the parental species D. buzzatii, D. koepferae, that most genes did not have regulatory divergence between and their hybrids, constituting a pioneer study in the parental species, which is consistent with the small percent- Drosophila genus. We found that most of the genomic age (5.92%) of differentially expressed genes found between regions in hybrids show similar histone modification patterns D. buzzatii and D. koepferae. We also detected more trans- compared with the parental species. This is consistent with regulatory than cis-regulatory divergence between these spe- the maintenance of the parental histone modification pat- cies, as well as a higher number of differentially expressed terns, as well as their additive inheritance in intraspecific trans-regulated genes. These results are in concordance (Moghaddam et al. 2011; Dong et al. 2012; Zhu, Greaves, with the high trans-regulatory divergence reported between et al. 2017) and interspecific (Zhu, Hu, et al. 2017) plant D. melanogaster and D. sechellia species (McManus et al. hybrids described in the literature. However, we also observe that 1.40–2.83% of genes and 1.37–9.27% of TE families are 2010) and the association of this class with gene differential differentially enriched in hybrids compared with any parental expression between parental species, but opposite to other species. H3K4me3 seems to be related with most gene ex- studies in D. melanogaster and D. simulans (Wittkopp et al. pression changes, consistent with the findings in rice and 2004; Graze et al. 2009). When we studied how the regula- maize hybrids (He et al. 2010, 2013) where gene expression tory divergence affects the gene inheritance in hybrids, we changes were correlated to the enrichment of this epigenetic found that most of the differentially expressed genes in mark. hybrids versus parental species had trans-regulatory diver- In the case of TEs, TE family expression changes in hybrids gence. These results are in contrast to what was reported in in comparison with D. buzzatii were related to their corre- the additive expressed genes, where cis-regulatory changes sponding enrichment changes in H3K4me3, H3K9me3, and were more frequent than trans-regulatory changes (Lemos H3K27me3, whereas they were related only to H3K9me3 and et al. 2008). Finally, we also found a high number of deregu- H3K27me3 enrichment changes when hybrids were com- lated genes in hybrids with cis–trans compensatory evolution, pared with D. koepferae. These results are in disagreement which was also previously reported in D. melanogaster and with previous results reported in interspecific Arabidopsis D. simulans hybrids and considered an important cause of hybrids, where changes in H3K9me2 and H3K27me3 do hybrid deregulation (Landry et al. 2005). The differences not coincide with the TEs having their expression changed reported, both between previous works and our own, could (Göbel et al. 2018). These differences could be explained by be due to use of different Drosophila species, different meth- the small number of differentially expressed TEs found in their odologies and the use, in this work, of ovarian tissues versus hybrids and the use of a different organism and methodology. whole adults. Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 13 n et al. Bodelo GBE We as well examined the expression of piRNA pathway genomic deregulation in hybrids. Nevertheless, we cannot genes for their role in germline TE silencing and we found discard additional mechanisms, resulting from the incompat- that aub and SoYb were both underexpressed, whereas Yb ibility of the two different paternal genomes in the hybrids, was overexpressed in hybrids versus both parental species. which could interact forming a complex gene network and Results of expression of these genes showed the same general contribute to the deregulation patterns observed. This and the trend as observed in a previous study with these species using fact that this study is limited to hybrid females, makes that a different analysis approach of RNAseq data (Romero- additional studies are necessary to go deeper into a better Soriano et al. 2017). Nevertheless, overexpression of three knowledge of the regulatory mechanisms and the factors in- additional piwi pathway genes in the same hybrids was volved in hybrid male sterility. detected by RT-qPCR in our previous work (Gamez-Visairas et al. 2020), suggesting that the activation of these genes Materials and Methods could be a primary response to the hybridization stress. Discrepancies found between RT-qPCR and RNAseq results Drosophila Stocks and Crosses highlight the different sensibilities of these two techniques. We used D. buzzatii Bu 28 and D. koepferae Ko2 inbred In addition, we found that the chromatin mark levels of the strains described in our previous works (Vela et al. 2014; piRNA pathway genes was similar to the parental species, Garcıa Guerreiro 2015; Romero-Soriano et al. 2016; with only differences in a few genes, which is consistent to Romero-Soriano and Garcia Guerreiro 2016). Both strains the general trend observed in other genes. In addition, we were maintained by brother–sister mating for at least a de- found that the piRNA amount changes were not associated cade and then by mass culturing. We performed ten different with TE expression changes in hybrids, because we found an interspecific crosses of ten D. buzzatii males to ten intermediate inheritance of piRNAs between the parental spe- D. koepferae females, due to the scarce offspring obtained cies in the hybrid, being D. koepferae the parental species in interspecific crosses. All stocks and crosses were reared at with the highest amount, as observed in our previous work 25  C in a standard Drosophila medium supplemented with (Romero-Soriano et al. 2017). yeast. We suggest that changes in transcript amounts in hybrids are either the result of the enrichment/impoverishment of a Chromatin Preparation specific mark or the disproportion in the active/repressive mark content, together with other mechanisms such as cis– Ovaries of 5-day-old Drosophila females were dissected in PBT trans compensatory regulation. Three histone marks, consid- (1 phosphate-buffered saline and 0.2% Tween 20). We ered as relevant to the expression in Drosophila, were consid- performed chromatin extraction from two biological repli- ered in this work, but we cannot discard the effect of other cates of 50 ovaries per sample of parental species and 70 epigenetic marks like H3K9ac highly positively correlated to ovaries for hybrids. Samples were resuspended in Buffer A1 gene expression (Yin et al. 2011). Other factor that could (HEPES 15 mM, sodium butyrate 10 mM, KCl 60 mM, Triton affect gene and TE deregulation is asynapsis, frequently ob- x100 0.5%, NaCl 15 mM) plus 1.8% formaldehyde, homog- served in our hybrids and reported in other repleta hybrids enized with a Dounce tissue grinder (15 times) and incubated (Naveira and Fontdevila 1991). Asynapsis, is known for for 10 min at room temperature. The crosslink was then influencing the trans regulation of Ultrabithorax gene alleles stopped with glycine to a final concentration of 125 mM. in Drosophila (Goldsborough and Kornberg 1996), for con- Samples were subsequently incubated 3 min, kept on ice tributing to mouse intersubspecific hybrids infertility and washed three times with Buffer A1. We then added (Bhattacharyya et al. 2013) and to the female meiotic losses 0.2 ml of lysis buffer (HEPES 15 mM, EDTA 1 mM, EGTA in mammals (Cloutier et al. 2015). 0.5 mM, sodium butyrate 10 mM, SDS 0.5%, sodium deoxy- cholate 0.1%, N-lauroylsarcosine 0.5%, Triton x100 1%, and NaCl 140 mM) and incubated 1–2 h at 4  C. After the lysis Conclusions process, we sonicated the samples using Biorruptor pico son- The hybridization between D. buzzatii and D. koepferae spe- ication device from Diagenode: 32 cycles of 30 s ON/30 s OFF, cies promotes a higher deregulation in TE families than in for parental species, and 35 cycles of 30 s ON/30 s OFF for genes, both toward underexpression. The epigenome of the hybrids. The sheared crosslinked chromatin was recovered parental species is in general highly preserved in the hybrids, from the pellet after a spin step at 10,000 g 4 min at 4  C. but some changes of the parental chromatin landscape are To check DNA size, samples were previously de-crosslinked also observed in hybrids and are associated with their new with NaCl 5M, boiled 15 min and treated with 1 ml of 10 mg/ gene and TE family expression patterns. Finally, cis–trans com- ml RNAse A (37  C, 30 min). They were purified with phenol- pensatory regulation could also be involved in expression de- chloroform, precipitated in absolute ethanol, washed in 70% regulation of some genes. The present study has contributed ethanol, resuspended in 20 ml of DEPC and run in a 1.5%- to a better understanding of the mechanism affecting agarose gel. 14 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE For the immunoprecipitation step, Magna ChIPTM G chro- values were computed using computeMatrix and visualized matin immunoprecipation Kit (Millipore) was used together with plotHeatmap. with antibodies against H3K4me3 (Abcam; ab8580), H3K9me3 (Abcam; ab8898), and H3K27me3 (Abcam; Gene Alignments ab6002). We separated 20 ml of chromatin for input and RNAseq reads from Romero-Soriano et al. (2017) were the remaining 180 ml were distributed in three aliquots where treated the same way as ChIPseq reads to ensure that results 1 ml of each antibody plus 20 ml of magnetic beads and dilu- were comparable. Gene sequences (only body region: from tion buffer up to 530 ml were added. Samples were incubated the SC to the EC) from D. buzzatii were obtained with get- overnight at 4  C in an agitation wheel. After the beads were fasta of Bedtools v2.29.2 (Quinlan and Hall 2010) using the washed with low salt buffer, high salt buffer, LiCl buffer, and genome sequence and its gene annotations (Guillen et al. TE buffer. Separation of chromatin from the beads was per- 2014). Gene sequences were masked using RepeatMasker formed using 0.5 ml of CHIP elution buffer. One microliter of v4.1.1. (http://www.repeatmasker.org) and the Repeat Proteinase K was then added to each sample and incubated at Masker from the D. buzzatii browser (http://dbuz.uab.cat, 62  C for 2 h in a shaker at 300 rpm. Samples were then last accessed January 7, 2015). A total of 37 genes were purified with the columns provided by the kit and stored at completely masked and excluded from further analysis, and 20  C. ChIP enrichment was quantified by real-time PCR of a total of 13,621 protein coding genes were included in the a well-known genomic region enriched for the different his- reference sequences. tone marks studied: rp49 for H3K4me3, kkv for H3K9me3, To ensure that there was no bias due to the use of only the and light for H3K27me3. The following gene-specific primers D. buzzatii reference genome, a de novo transcriptome for were used: rp49-forward: 50 -GTCGTCGCTTCAAGGGCCAA each parental species was created using Trinity v2.9 (Grabherr T-30 , rp49-reverse: 50 -ATGGGCGATCTCACCGCAGTA-30 , et al. 2011) and the corresponding RNAseq. Drosophila buz- kkv-forward: 50 -TAATCCAGCCACGCCCATTT-30 , kkv- zatii and D. koepferae transcriptomes were then aligned reverse: 5 -CCCAACGTTTGCATTGCTGA-30 , light-forward: 0 against all the ChIPseq inputs (supplementary file 4: 50 -CGAGTACAAAATGAATAGCTCCG-30 , light-reverse: 5’- Alignment_inputChIPSeq, Supplementary Material online). GCGGTTCTCCTCAATGAT-30 . We also randomly selected 40 genes amongst Trinity outputs, and computed nucleotide divergence between parental spe- Chromatin Sequencing cies using the Mega software (Tamura et al. 2021) and the Duplicate Truseq ChIP libraries, corresponding to two biolog- Jukes-Cantor model (Jukes and Cantor 1969) (supplementary ical replicates per sample, were performed by Macrogen Inc., file 4: Nucleotide_Divergence, Supplementary Material on- Seoul, Korea. Sequencing was carried out using an Illumina line). Both, the similar alignment percentages and the low Hiseq4000. We obtained 22–34 millions of paired-end reads average divergence, indicate that bias of using D. buzzatii as for each sample, resulting in a total of 659.9 millions of reference genome, if any, is marginal. Nevertheless, we de- paired-end reads. cided to be conservative and use only the protein coding genes. Trimmed RNAseq and ChIPseq reads were aligned to ChIPseq Visualization Patterns the D. buzzatii masked gene sequences using Bowtie2 v2.3.5.1 (Langmead and Salzberg 2012) and the parameters ChIPseq sequenced reads were trimmed using Trimmomatic explained above. eXpress v1.5.1 (Roberts and Pachter 2013) software v0.39 (Bolger et al. 2014) and aligned to the with the default options and then an additional online EM D. buzzatii genome downloaded from the Drosophila buzzatii round (to increase the accuracy) was used to quantify read Genome Project web page (http://dbuz.uab.cat, last accessed counts for each gene. All isoforms of a gene were considered January 7, 2015) using Bowtie2 v2.3.5.1 (Langmead and together. We used rounded effective counts for the following Salzberg 2012). For the alignment, the default parameters steps. of the –very-sensitive-local mode with two extra- modifications to increase the sensibility (–local -D 20 -R 3 -L 20 -N 1 -p 12 –gbar 3 –mp 5,1 –rdg 4,2 –rfg 4,2) were used to TE Alignment reach the highest percentage of aligned reads with both pa- TE RNAseq (Romero-Soriano et al. 2017) and ChIPseq read rental species and their F1 hybrids. Reads with map quality counts were analyzed using the TEcount module of the score lower than 30 and unmapped reads were filtered using TEtools pipeline (Lerat et al. 2017) (available at https:// Samtools v1.10 (Li et al. 2009) and excluded from further github-com/l-modolo/TEtools). First, the manually constructed analysis. Deeptools v3.3.2 (Ramırez et al. 2016) was used to TE library containing consensus TEs from both D. buzzatii and visualize the enrichment of each chromatin mark around the D. koepferae, described in a previous work (Romero-Soriano SC and the end codon (EC). First, bamCompare was used to et al. 2017), was aligned with the RNAseq and ChIPseq reads, normalize the ChIPseq samples by depth using the RPKM using Bowtie 2 v2.2.4 (Langmead and Salzberg 2012) with method and by the input (control). Finally, the read density the most sensitive option and keeping a single alignment for Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 15 n et al. Bodelo GBE reads mapping to multiple positions (–very-sensitive). Read TE families with the lowest P-values were considered as dif- counts were computed per TE family, adding all reads ferentially enriched. mapped on copies from the same family. Count tables corre- The significant genes were assigned to their respective sponding to genes and TEs were concatenated and were used chromosomes following the scaffold assignation to chromo- for the differential expression and enrichment analyses. Genes somes obtained from a previous work (Guillen et al. 2014) to counts were used to normalize TE counts, following the detect potential chromosome biases of differentially guidelines of TEtools pipeline (Lerat et al. 2017). expressed genes. Additionally, the order of the significant TE families was analyzed to detect possible order biases in Differential Expression and Enrichment Analyses our significant TEs. Within each genome, we quantified the contribution of The statistical analyses were performed using R v4.0 (R Core histone mark enrichments on transcript amounts using the Team 2020). The visualization of the results was performed following linear model on log-transformed read counts for using the R package ggplot2 v3.3.2 (Wickham 2016). The each gene and TE family: RNA  K4 þ K9 þ K27 þ Input. In DESeq2 function from the R Bioconductor package DESeq2 addition, we tested whether changes in RNAseq counts in v1.28.1 (Love et al. 2014) was used to normalize read counts, hybrids were associated with changes in ChIPseq counts using using the default method, and to model the read counts using Fisher’s exact test on 2  2 matrices: genes and TE families a negative binomial distribution. were classified as displaying a positive or negative log2FC of For RNAseq, the DESeq2 package was used to identify RNAseq counts in hybrids versus parents and as displaying an differentially expressed genes and TE families while perform- increase or decrease in histone mark enrichment in hybrids ing a Wald test (Love et al. 2014). The P-values were adjusted versus parents. for multiple testing using the procedure of Benjamin and Hochberg (Benjamini et al. 2001) with an FDR cutoff of 0.05, and were obtained using the results function from the Allele-Specific Expression Analysis DESeq2 package. The log2FC was shrunken using the default We created a de novo transcriptome with the RNAseq of both and recommended apeglm algorithm (Zhu et al. 2019) of the parental species, together with the hybrids using Trinity v2.9 lfcShrink function. Genes with an adjusted P-value lower than (Grabherr et al. 2011). A SuperTranscript (Davidson et al. 0.05 and at least double of expression in one species above 2017) was then created as a general reference transcriptome. the other (absolute shrunken log2FC > 1) were considered as The GATK pipeline for variant calling (Van der Auwera and differentially expressed. GO term enrichment analyses of bio- O’Connor 2020) was used to detect differences between the logical processes were performed for the underexpressed and general reference transcriptome and each parental species. overexpressed significant genes using the topGO R package The VCF files were filtered following the GATK guidelines, v.2.42.0. (Alexa and Rahnenfuhrer 2020) (“weight01” algo- including a coverage depth of at least 10, and the exclusion rithm and Fisher’s statistic) and the ortholog genes in of variants only present in one replicate, which were consid- D. melanogaster obtained from the D. buzzatii Genome ered assembly errors. These variants were replaced in the Project web page (http://dbuz.uab.cat, last accessed January general reference transcriptome using the 7, 2015). The obtained Fisher’s exact test P-values were not FastaAlternateReferenceMaker GATK tool (Van der Auwera adjusted but, from the total gene ontologies with P < 0.05, and O’Connor 2020) to create a reference transcriptome for only the top 13 with the lowest P-values were considered as each parental species. enriched. HISAT2 (Kim et al. 2015) with the “no-softclip” option was Regarding the ChIPseq data, a handmade script was used used to align the RNAseq from hybrids to both parental ref- to perform the differential enrichment analysis between erence transcriptome. Then, CompareBams (Lindenbaum hybrids and both parental species. Briefly, the “Regularized 2015) of Jvarkit was used to compare the alignments and log” transformed counts obtained using DESeq2 for each FilterSamReads from Picard (Broad Institute 2018) to filter chromatin mark were considered and analyzed using the fol- out the reads aligning in different position when the data lowing linear model for each gene: RldKm  RldKi þ species were aligned to each parental species. Samtools v1.10 (Li (RldKm: log-transformed counts for the considered histone et al. 2009) was used to remove multimapped and unmapped mark, RldKi: log-transformed counts for the input). We con- reads and BamTools v.2.4.0 (Barnett et al. 2011) to keep only sidered contrasts between pairs of species modalities reads that align without mismatches. Bedtools multicov (D. buzzatii, D. koepferae, hybrid). Because of the high num- (Quinlan and Hall 2010) and a manually updated GTF were ber of tests and the lack of extreme P-values, the Benjamin used to count reads aligning to each gene. and Hochberg (Benjamini et al. 2001) adjustment of the The reference transcriptomes were annotated using BLAT P-values had a high effect in our results. For this reason, the v.35x1 (Kent 2002) with the parameters –minIdentity ¼ 80 P-values were not adjusted and, from the genes and TE and –maxIntron ¼ 75,000 and the gene sequences of families with P < 0.05, the rounded 13 of the top genes and D. buzzatii, keeping the best match with an overlap of at least 16 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE 50%. Subsequently, statistical analyses were performed using (Langmead et al. 2009) and keeping a single alignment for R v4.0 (R Core Team 2020) and DESeq2. The DESeq2 function reads mapping to multiple positions. These reads were con- was used to normalize read counts, using the default method. sidered to correspond to microRNAs and counts were com- First, the results function was used to compute the log2FC puted using eXpress v1.5.1 (Roberts and Pachter 2013). TE and the adjusted P-values using the Benjamin and Hochberg counts among piRNAs were computed using the TEcount (Benjamini et al. 2001) method. A FDR cutoff of 0.05 was module of TEtools (Lerat et al. 2017). piRNA counts were used to identify differentially expressed genes between paren- then normalized so that the sum of microRNA counts is the tal species (P). Second, the same procedure and cutoff was constant across samples. used to identify genes with a different abundance of the pa- Regarding TE sequences, RNAseq and CHIPseq data were rental and the maternal allele in hybrids (H), and were con- integrated following the same approaches and processes as sidered as genes with cis-regulatory divergence. Finally, for genes. To get insight into the chromatin mark combina- significant genes in either P or H were analyzed for trans- tion, we study the association between the euchromatic mark regulatory effects (T) by comparing the P and H ratios with H3K4me3 and the heterochromatic marks (H3K9me3 and the same results function. Genes were then categorized in the H3K27me3) in hybrids versus each parental species. following groups as reported in a previous work (McManus et al. 2010): Statistical Tests • Cis- only: Significant differential expression in P and H, but Three main statistical tests were used in the article and were no significant T. performed using the R v4.0 (R Core Team 2020) program: • Trans- only: Significant differential expression in P, and T, but not H. • Two proportion Z-test was used to compare distributions of • Cis- þ trans-: Significant differential expression in P, H, and significant sequences (genes and TE families) across com- T. Cis- and trans- regulatory differences favor expression of parisons and cis- and trans-regulatory classes. the same allele. • Chi-square test under equal assumption was used to detect • Cis-  trans-: Significant differential expression in P, H, and chromosome-biases and TE-category-biases. T. Cis- and trans- regulatory differences favor expression of • Fisher’s exact test under independence assumption was opposite alleles. computed using a 2  2 contingency table to detect asso- • Compensatory: Significant differential expression in H, and ciations in: gene and TE expression (log2FC) and chromatin T, but not P. Cis- and trans-regulatory differences compen- mark enrichment (size effect); TE euchromatic mark sate each other, resulting in no expression differences be- (H3Kme3) and heterochromatic marks (H3K9me3 and tween parental species. H3K27me3); trans- and cis-regulatory divergence and dif- • Conserved: No significant differential expression in H, P, or ferences of expression between parental species; compen- T. Conserved regulation. satory and remaining classes and gene deregulation in • Ambiguous: Other patterns with no clear biological hybrids, and TE expression and piRNA quantities, in hybrids interpretation. versus each parental species independently. The results were corrected for multiple testing using the Piwi Pathway Genes and piRNA Amount Study Benjamin and Hochberg (Benjamini et al. 2001) method. Additionally, reciprocal tblast v2.10.1. (Camacho et al. 2009) of the D. melanogaster piRNA pathway proteins from UniProt (The UniProt Consortium 2021) and the reference gene se- Supplementary Material quence list were performed. A total of 31 proteins were Supplementary data are available at Genome Biology and detected and associated with D. buzzatti genes. The Evolution online. Argonaute 3 gene (ago3), not included in the reference genes, was manually included in the list from our previous Acknowledgments results (Gamez-Visairas et al. 2020). Small RNAseq raw reads from our previous work (Romero- The authors wish to thank Victor Gamez and Judit Salces- Soriano et al. 2017) were used to study the piRNA regulatory Ortiz for their contribution to the experimental part of this data in the TE expression study. Using PRINSEQ lite (Schmieder work. They want to thank as well to Daniel Siqueira de and Edwards 2011), we isolated 23–30 nt-long reads and Oliveira for his contribution in the DESeq pipeline and considered them as piRNAs. For normalization purposes, we Mariana Galvao Ferrarini for her advises in GO analyses. also isolated 20–23 nt-long reads and searched for microRNA Finally, they wish to thank Lars Ootes for his advices with sequences: low-quality reads were removed using UrQt the allele-specific expression analyses and two anonymous (Modolo and Lerat 2015), then trimmed 20–23 nt reads reviewers for their valuable comments on the manuscript. were aligned to the masked genes using Bowtie V.1.3.0 This work was supported by grants CGL2017-89160P from Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 17 n et al. Bodelo GBE Ministerio de Economıa y Competitividad (Spain; cofinanced Chan WY, Chung J, Peplow LM, Hoffmann AA, van Oppen MJH. 2021. with the European Union FEDER funds), 2017SGR 1379 from Maternal effects in gene expression of interspecific coral hybrids. Mol Ecol. 30(2):517–527. Generalitat de Catalunya, and ANR ExHyb from the French Cloutier JM, et al. 2015. Histone H2AFX links meiotic chromosome asyn- National Research Agency. A.B. was supported by a PIF pre- apsis to prophase I oocyte loss in mammals. PLoS Genet. doctoral fellowship from the Universitat Auto noma de 11(10):e1005462. Barcelona (Spain). Coyne JA, Orr HA. 1989. Two rules of speciation. In: Otte D, Endler JA, editors. Speciation and its consequences. Sunderland (MA): Sinauer Associates. p. 180–207. Author Contributions Davidson NM, Hawkins ADK, Oshlack A. 2017. SuperTranscripts: a data driven reference for analysis and visualisation of transcriptomes. A.B.: data analysis, interpretation, and article writing; Genome Biol. 18(1):148. M.P.G.G. and C.V.: conception, design, interpretation, and De Araujo PG, et al. 2005. Transcriptionally active transposable elements in article writing; M.F.: data analysis, interpretation, and article recent hybrid sugarcane. Plant J. 44(5):707–717. Dion-Côt e AM, Renaut S, Normandeau E, Bernatchez L. 2014. RNA-seq writing; PV: data analysis. All authors read and approved the reveals transcriptomic shock involving transposable elements reactiva- final manuscript. tion in hybrids of young lake whitefish species. Mol Biol Evol. 31(5):1188–1199. Dong X, et al. 2012. Natural variation of H3K27me3 distribution between Data Availability two Arabidopsis accessions and its association with flanking transpos- The data are included in the article and in its Supplementary able elements. Genome Biol. 13(12):R117. Ebert A, Lein S, Schotta G, Reuter G. 2006. Histone modification and the Material online. Raw ChIPseq data generated from this article control of heterochromatic gene silencing in Drosophila. Chromosome have been deposited at the NCBI Sequence Read Archive Res. 14(4):377–392. (SRA) under the BioProject accession PRJNA796032: accession Filion GJ, et al. 2010. Systematic protein location mapping reveals five numbers from SRR17535990 to SRR17536013. principal chromatin types in Drosophila cells. Cell 143(2):212–224. Fontdevila A. 2019. Hybrid genome evolution by transposition: an update. J Hered. 110(1):124–136. Literature Cited Gamez-Visairas V, Romero-Soriano V, Martı-Carreras J, Segarra-Carrillo E, Garcıa Guerreiro MP. 2020. Drosophila interspecific hybridization Alexa A, Rahnenfuhrer J. 2020. topGO: enrichment analysis for gene on- causes a deregulation of the piRNA pathway genes. Genes 11(2):215. tology. R package version 2.42.0. Available from: https://bioconduc- Garcıa Guerreiro MP. 2015. Changes of Osvaldo expression patterns in tor.org/packages/release/bioc/html/topGO.html germline of male hybrids between the species Drosophila buzzatii and Allan CW, Matzkin LM. 2019. Genomic analysis of the four ecologically Drosophila koepferae. Mol Genet Genomics. 290(4):1471–1483. distinct cactus host populations of Drosophila mojavensis. BMC Garcıa Guerreiro MP. 2012. What makes transposable elements move in Genomics. 20(1):732. the Drosophila genome. Heredity 108(5):461–468. Barnett DW, Garrison EK, Quinlan AR, Strömberg MP, Marth GT. 2011. Garcıa Guerreiro MP. 2014. Interspecific hybridization as a genomic BamTools: a Cþþ API and toolkit for analyzing and managing BAM stressor inducing mobilization of transposable elements in files. Bioinformatics 27(12):1691–1692. Drosophila. Mob Genet Elements. 4:e34394. Benjamini Y, Drai D, Elmer G, Kafkafi N, Golani I. 2001. Controlling the Göbel U, et al. 2018. Robustness of transposable element regulation but false discovery rate in behavior genetics research. Behav Brain Res. no genomic shock observed in interspecific Arabidopsis hybrids. 125(1–2):279–284. Genome Biol Evol. 10(6):1403–1415. Bhattacharyya T, et al. 2013. Mechanistic basis of infertility of mouse Goldsborough AS, Kornberg TB. 1996. Reduction of transcription by ho- intersubspecific hybrids. Proc Natl Acad Sci USA. 110:E468–E477. mologue asynapsis in Drosophila imaginal discs. Nature Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: a flexible trimmer for 381(6585):807–810. Illumina sequence data. Bioinformatics 30(15):2114–2120. Go  mez GA, Hasson E. 2003. Transpecific polymorphisms in an inversion Boros IM. 2012. Histone modification in drosophila. Brief Funct Genomics. linked esterase locus in Drosophila buzzatii. Mol Biol Evol. 11(4):319–331. 20(3):410–423. Broad Institute. 2018. Picard toolkit. Broad Institute, GitHub Repos. Grabherr MG, et al. 2011. Full-length transcriptome assembly from RNA- Available from: https://broadinstitute.github.io/picard/. Accessed Seq data without a reference genome. Nat Biotechnol. February 21, 2018. 29(7):644–652. Brown EJ, Bachtrog D. 2014. The chromatin landscape of Drosophila: Graze RM, McIntyre LM, Main BJ, Wayne ML, Nuzhdin SV. 2009. comparisons between species, sexes, and chromosomes. Genome Regulatory divergence in Drosophila melanogaster and D. simulans, Res. 24(7):1125–1137. a genomewide analysis of allele-specific expression. Genetics Camacho C, et al. 2009. BLASTþ: architecture and applications. BMC 183(2):547–561. Bioinformatics. 10(1):1–9. Guillen Y, et al. 2014. Genomics of ecological adaptation in cactophilic Carnelossi EAG, et al. 2014. Specific activation of an I-like element in Drosophila. Genome Biol Evol. 7(1):349–366. Drosophila interspecific hybrids. Genome Biol Evol. Guio L, Vieira C, Gonz alez J. 2018. Stress affects the epigenetic marks 6(7):1806–1817. added by natural transposable element insertions in Drosophila mela- Center for Research in Agricultural Genomics, National Center for nogaster. Sci Rep. 8:1–10. Genomic Analysis, Genomics Bioinformatics and Evolution Research He G, et al. 2010. Global epigenetic and transcriptional trends among two Group of the Universitat Auto  noma de Barcelona, Barcelona rice subspecies and their reciprocal hybrids. Plant Cell. 22(1):17–33. Supercomputing Center. Drosophila buzzatii Genome Project. He G, et al. 2013. Conservation and divergence of transcriptomic and Available from: http://dbuz.uab.cat (accessed November 26, 2020). epigenomic variation in maize hybrids. Genome Biol. 14(6):R57. 18 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 High Stability of the Epigenome in Drosophila Interspecific Hybrids GBE Henault M, Marsit S, Charron G, Landry CR. 2020. The effect of hybrid- Modolo L, Lerat E. 2015. UrQt: an efficient software for the unsupervised ization on transposable element accumulation in an undomesticated quality trimming of NGS data. BMC Bioinformatics. 16(1):1–8. fungal species. eLife 9:e60474. Moehring AJ, Teeter KC, Noor MAF. 2007. Genome-wide patterns of Hollocher H, Wu CI. 1996. The genetics of reproductive isolation in the expression in Drosophila pure species and hybrid males. II. Drosophila simulans clade: x vs. autosomal effects and male vs. female Examination of multiple-species hybridizations, platforms, and life cy- effects. Genetics 143(3):1243–1255. cle stages. Mol Biol Evol. 24(1):137–145. Jukes T, Cantor C. 1969. Evolution of protein molecules. New York: Moghaddam AMB, et al. 2011. Additive inheritance of histone modifica- Academic press. tions in Arabidopsis thaliana intra-specific hybrids. Plant J. Jurka J, et al. 2005. Repbase update, a database of eukaryotic repetitive 67(4):691–700. elements. Cytogenet Genome Res. 110(1–4):462–467. Montgomery SA, et al. 2020. Chromatin organization in early land plants Kelleher ES, Edelman NB, Barbash DA. 2012. Drosophila interspecific reveals an ancestral association between H3K27me3. Curr Biol. hybrids phenocopy piRNA-pathway mutants. PLoS Biol. 30(4):573–588. 10(11):e1001428. Naveira H, Fontdevila A. 1991. The evolutionary history of Drosophila Kent WJ. 2002. BLAT–the BLAST-like alignment tool. Genome Res. buzzatii. XXI. Cumulative action of multiple sterility factors on sperma- 12(4):656–664. togenesis in hybrids of D. buzzatii and D. koepferae. Heredity Kim D, Langmead B, Salzberg SL. 2015. HISAT: a fast spliced aligner with 67(1):57–72. low memory requirements. Nat Methods. 12(4):357–360. Naveira H, Pla C, Fontdevila A. 1986. The evolutionary history of Labrador M, Farre M, Utzet F, Fontdevila A. 1999. Interspecific hybridiza- Drosophila buzzatii XI. A new method for cytogenetic localization tion increases transposition rates of Osvaldo. Mol Biol Evol. based on asynapsis of polytene chromosomes in interspecific hybrids 16(7):931–937. of Drosophila. Genetica 71(3):199–212. Landry CR, et al. 2005. Compensatory cis-trans evolution and the dysre- Oliveira DCSG, et al. 2012. Monophyly, divergence times, and evolution of gulation of gene expression in interspecific hybrids of Drosophila. host plant use inferred from a revised phylogeny of the Drosophila Genetics 171(4):1813–1822. repleta species group. Mol Phylogenet Evol. 64(3):533–544. Langmead B, Salzberg SL. 2012. Fast gapped-read alignment with Bowtie Ortız-Barrientos D, Counterman BA, Noor MAF. 2007. Gene expression 2. Nat Methods. 9(4):357–359. divergence and the origin of hybrid dysfunctions. Genetica Langmead B, Trapnell C, Pop M, Salzberg SL. 2009. Ultrafast and memory- 129(1):71–81. efficient alignment of short DNA sequences to the human genome. Parey E, Crombach A, Alba M. 2019. Evolution of the Drosophila mela- Genome Biol. 10(3):r25. nogaster chromatin landscape and its associated proteins. Genome Lemos B, Araripe LO, Fontanillas P, Hartl DL. 2008. Dominance and the Biol Evol. 11(3):660–677. evolutionary accumulation of cis- and trans-effects on gene expres- Quinlan AR, Hall IM. 2010. BEDTools: a flexible suite of utilities for com- sion. Proc Natl Acad Sci USA. 105(38):14471–14476. paring genomic features. Bioinformatics 26(6):841–842. Lerat E, Fablet M, Modolo L, Lopez-Maestre H, Vieira C. 2017. TEtools R Core Team. 2020. R: a language and environment for statistical com- facilitates big data expression analysis of transposable elements and puting. Vienna (Austria): R Foundation for Statistical Computing. reveals an antagonism between their activity and that of piRNA genes. Available from: https://www.r-project.org/. Nucleic Acids Res. 45(4):e17. Ramırez F, et al. 2016. deepTools2: a next generation web server for deep- Li H, et al. 2009. The sequence alignment/map format and SAMtools. sequencing data analysis. Nucleic Acids Res. 44(W1):W160–W165. Bioinformatics 25(16):2078–2079. Ranz JM, Namgyal K, Gibson G, Hartl DL. 2004. Anomalies in the expres- Lindenbaum P. 2015. JVarkit: java-based utilities for Bioinformatics. fig- sion profile of interspecific hybrids of Drosophila melanogaster and share. doi: https://doi.org/10.6084/m9.figshare.1425030. Drosophila simulans. Genome Res. 14(3):373–379. Lopez-Maestre H, et al. 2017. Identification of misexpressed genetic ele- Rebollo R, et al. 2012. A snapshot of histone modifications within trans- ments in hybrids between Drosophila-related species. Sci Rep. posable elements in Drosophila wild type strains. PLoS One. 7:40618–40613. 7(9):e44253. Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change Renaut S, Rowe HC, Ungerer MC, Rieseberg LH. 2014. Genomics of and dispersion for RNA-seq data with DESeq2. Genome Biol. homoploid hybrid speciation: diversity and transcriptional activity of 15(12):550–521. long terminal repeat retrotransposons in hybrid sunflowers. Philos Mack KL, Nachman MW. 2017. Gene regulation and speciation. Trends Trans R Soc Lond B Biol Sci. 369 (1648):20130345. Genet. 33(1):68–80. Riddle NC, et al. 2011. Plasticity in patterns of histone modifications and Malone JH, Michalak P. 2008. Gene expression analysis of the ovary of chromosomal proteins in Drosophila heterochromatin. Genome Res. hybrid females of Xenopus laevis and X. muelleri. BMC Evol Biol. 21(2):147–163. 8(1):82. Roberts A, Pachter L. 2013. Streaming fragment assignment for real-time Marin I, Fontdevila A. 1998. Stable Drosophila buzzatii–Drosophila koep- analysis of sequencing experiments. Nat Methods. 10(1):71–73. ferae hybrids. J Hered. 89(4):336–339. Romero-Soriano V, et al. 2016. Drosophila females undergo genome ex- Marin I, Ruiz A, Pla C, Fontdevila A. 1993. Reproductive relationships pansion after interspecific hybridization. Genome Biol Evol. among ten species of the Drosophila repleta group from South 8(3):556–561. America and the West Indies. Evolution 47(5):1616. Romero-Soriano V, et al. 2017. Transposable element misregulation is McManus CJ, et al. 2010. Regulatory divergence in Drosophila revealed by linked to the divergence between parental piRNA pathways in mRNA-seq. Genome Res. 20(6):816–825. Drosophila hybrids. Genome Biol Evol. 9(6):1450–1470. Michalak P, Noor MAF. 2003. Genome-wide patterns of expression in Romero-Soriano V, Garcia Guerreiro MP. 2016. Expression of the retro- Drosophila pure species and hybrid males. Mol Biol Evol. transposon helena reveals a complex pattern of TE deregulation in 20(7):1070–1076. Drosophila hybrids. PLoS One. 11(1):e0147903. modENCODE Consortium, et al. 2011. Identification of functional ele- Saha P, Sowpati DT, Mishra RK. 2019. Epigenomic and genomic landscape ments and regulatory circuits by Drosophila modENCODE. Science of Drosophila melanogaster heterochromatic genes. Genomics 330:1787–1797. 111(2):177–185. Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022 19 n et al. Bodelo GBE Satyaki PRV, et al. 2014. The Hmr and Lhr hybrid incompatibility genes Videvall E, Sletvold N, Hagenblad J, Ågren J, Hansson B. 2016. Strong suppress a broad range of heterochromatic repeats. PLoS Genet. maternal effects on gene expression in Arabidopsis lyrata hybrids. 10(3):e1004240. Mol Biol Evol. 33(4):984–994. Schmieder R, Edwards R. 2011. Quality control and preprocessing of Walter M, Teissandier A, P erez-Palacios R, Bourc’his D. 2016. An epige- metagenomic datasets. Bioinformatics 27(6):863–864. netic switch ensures transposon repression upon dynamic loss of DNA Schübeler D, et al. 2004. The histone modification pattern of active genes methylation in embryonic stem cells. eLife 5:1–30. revealed through genome-wide chromatin analysis of a higher eukary- Wang J, et al. 2006. Genomewide nonadditive gene regulation in ote. Genes Dev. 18(11):1263–1271. Arabidopsis allotetraploids. Genetics 172(1):507–517. Schulze SR, Wallrath LL. 2007. Gene regulation by chromatin structure: Wickham H. 2016. ggplot2: elegant graphics for data analysis. New paradigms established in Drosophila melanogaster. Annu Rev York: Springer-Verlag. Available from: https://ggplot2.tidyverse. Entomol. 52:171–192. org. Soto EM, Soto IM, Carreira VP, Fanara JJ, Hasson E. 2008. Host-related life Wittkopp PJ, Haerum BK, Clark AG. 2004. Evolutionary changes in cis and history traits in interspecific hybrids of cactophilic Drosophila. Entomol trans gene regulation. Nature 430(6995):85–88. Exp Appl. 126:18–27. Wurmser F, et al. 2011. Population transcriptomics: insights from Tamura K, Stecher G, Kumar S. 2021. MEGA11: molecular evolu- Drosophila simulans, Drosophila sechellia and their hybrids. Genetica tionary genetics analysis version 11. Mol Biol Evol. 38(7): 139(4):465–477. 3022–3027. Yin H, Sweeney S, Raha D, Snyder M, Lin H. 2011. A high-resolution The UniProt Consortium. 2021. UniProt: the universal protein knowledge- whole-genome map of key chromatin modifications in the adult base in 2021. Nucleic Acids Res. 49:D480–D489. Drosophila melanogaster. PLoS Genet. 7(12):e1002380. Van der Auwera G, O’Connor B. 2020. GATK, and WDL in terra Zhu A, Greaves IK, Dennis ES, Peacock WJ. 2017. Genome-wide analyses genomics in the cloud: using Docker. 1st ed. Sebastopol (CA): of four major histone modifications in Arabidopsis hybrids at the ger- O’Reilly Media, Inc. minating seed stage. BMC Genomics. 18(1):1–9. Vela D, Fontdevila A, Vieira C, Garcıa Guerreiro MP. 2014. A genome- Zhu A, Ibrahim JG, Love MI. 2019. Heavy-tailed prior distributions for se- wide survey of genetic instability by transposition in Drosophila quence count data: removing the noise and preserving large differ- hybrids. PLoS One 9(2):e88992. ences. Bioinformatics 35(12):2084–2092. Vela D, Garcıa Guerreiro MP, Fontdevila A. 2011. Adaptation of Zhu W, Hu B, et al. 2017. Altered chromatin compaction and histone the AFLP technique as a new tool to detect genetic instability methylation drive non-additive gene expression in an interspecific and transposition in interspecific hybrids. Biotechniques 50(4): Arabidopsis hybrid. Genome Biol. 18(1):1–16. 247–250. Associate editor: Andrea Betancourt 20 Genome Biol. Evol. 14(2) https://doi.org/10.1093/gbe/evac024 Advance Access publication 10 February 2022