Skip to main content

Regulatory and evolutionary impact of DNA methylation in two songbird species and their naturally occurring F1 hybrids



Regulation of transcription by DNA methylation in 5’-CpG-3’ context is a widespread mechanism allowing differential expression of genetically identical cells to persist throughout development. Consequently, differences in DNA methylation can reinforce variation in gene expression among cells, tissues, populations, and species. Despite a surge in studies on DNA methylation, we know little about the importance of DNA methylation in population differentiation and speciation. Here we investigate the regulatory and evolutionary impact of DNA methylation in five tissues of two Ficedula flycatcher species and their naturally occurring F1 hybrids.


We show that the density of CpG in the promoters of genes determines the strength of the association between DNA methylation and gene expression. The impact of DNA methylation on gene expression varies among tissues with the brain showing unique patterns. Differentially expressed genes between parental species are predicted by genetic and methylation differentiation in CpG-rich promoters. However, both these factors fail to predict hybrid misexpression suggesting that promoter mismethylation is not a main determinant of hybrid misexpression in Ficedula flycatchers. Using allele-specific methylation estimates in hybrids, we also determine the genome-wide contribution of cis- and trans effects in DNA methylation differentiation. These distinct mechanisms are roughly balanced in all tissues except the brain, where trans differences predominate.


Overall, this study provides insight on the regulatory and evolutionary impact of DNA methylation in songbirds.


Genetic lineages that diverge and form new species acquire both genetic and phenotypic changes in the process. As a by-product, intrinsic reproductive isolation in the form of hybrid dysfunction may evolve. The genetic basis of hybrid dysfunction is usually interactions between incompatible alleles that have never been tested in the same cellular environment [1]. These so-called Bateson–Dobzhansky–Muller incompatibilities (BDMIs) [2,3,4] are believed to be important for the completion of speciation since they are not dependent on the environment [5]. While some progress have been made in mapping BDMI loci [6], the molecular mechanisms causing hybrid dysfunction is less clear [7]. Incompatible interactions between transcription factors and their binding sites is a class of BDMIs that may cause low fitness in hybrids by the misexpression of genes [8]. Misexpression in hybrids could also be caused by other gene regulatory mechanisms, e.g., aberrant levels of different epigenetic marks such as DNA methylation [9]. Given the large number of genes, transcription factors, and cis-regulatory elements in a vertebrate genome, gene regulatory mechanisms are likely to be an important albeit understudied source of BDMIs [10, 11]. With large genomic data sets, we can now both characterize genetic and epigenetic sources of gene regulatory variation and thus gain insight not only on the degree of hybrid misexpression but also the molecular basis of misregulation [12,13,14]. In this study, we specifically assess the role of DNA methylation in gene misexpression of naturally occurring hybrids.

In vertebrates, DNA methylation is most frequently occurring at cytosines in the 5’-CpG-3’ dinucleotide context and is often associated with transcriptional repression [15]. How DNA methylation regulates gene expression is explained by the “molecular lock” model [16]. Following de novo DNA methylation, a gene is locked in a silent state preventing further transcription until the DNA methylation marks are removed, either passively through lack of DNA methylation maintenance during replication or actively using enzymatic activity. In addition to silencing genes, DNA methylation in vertebrates is also used to silence the expression of transposable elements (TEs) [17,18,19]. Keeping genes (and TEs) repressed through DNA methylation is thought to occur through three main mechanisms that are not mutually exclusive: (1) DNA methylation prevents binding of transcription factors to target sequences [20], (2) DNA methylation acts as substrate for proteins mediating repression such as MeCP2 [21], or (3) de novo methylation alters the chromatin structure to a more compact inactive state [22]. Common to these mechanisms is that they are all predicted to yield a negative relationship between the level of DNA methylation in the promoter region and expression level. Gene body methylation on the other hand shows (if any) a positive or quadratic relationship with expression level, but the function of gene body methylation is debated [23,24,25,26,27].

While most of the genome is methylated in adult vertebrates, some sequence regions escape the wave of de novo methylation occurring during development [28]. These so-called CpG islands (CGI) are characterized by a high density of CpG dinucleotides compared to the level predicted by their GC content [29]. Typically, unmethylated CpG dinucleotides contribute to an open chromatin state permissible for transcriptional initiation [30]. Roughly 70% of all promoters in the human genome have CGIs [31]. In humans, most CpG-dense promoters (CGI promoters) associated with housekeeping genes remain unmethylated in adult tissues, while other CGI promoters are dynamically regulated [32]. In contrast, questions remain to what extent CpG-deficient promoters are regulated by DNA methylation [32, 33].

Differences in methylation caused by genetic changes among alleles, populations, and species can arise from two mechanisms: either cis-regulatory changes specific to the DNA sequence at a locus or trans-regulatory changes because of divergence in structure or regulation of interacting factors with potential to change methylation level. In model organisms, this relationship has been investigated using data from hundreds or thousands of individuals to determine quantitative trait loci affecting methylation level [34]. These regulatory mechanisms may also be distinguished using F1 hybrid individuals, which was pioneered in studies of gene expression divergence in mice and fruit flies [35, 36], but also recently applied to other molecular phenotypes [37, 38]. The hybrid represents a trans environment in which allele-specific measures of molecular phenotypes can be measured and contrasted with parental species. Here, we extend this approach to DNA methylation data and develop a statistical framework to infer the molecular mechanism of DNA methylation differentiation.

Specifically, in this study, we investigate the regulatory and evolutionary impact of DNA methylation in two Ficedula flycatcher species and their naturally occurring F1 hybrids. The pied flycatcher (Ficedula hypoleuca) and the collared flycatcher (F. albicollis) are two species of songbirds that diverged approximately 1 MYA [39]. Interspecific mating occurs in sympatry, e.g., on the island of Öland in the Baltic Sea [40]. F1 hybrids of both sexes are infertile [41, 42] and show reduced viability [43], and males have an increased metabolic rate [44, 45]. A number of studies has investigated DNA methylation in birds (e.g., [46,47,48,49,50,51,52,53,54]). In general, methylation patterns in birds are similar to other vertebrates. However, compared to mammals, much less is known about genome-wide DNA methylation patterns across tissues. Also, the interplay between DNA methylation, gene expression, and genetic variation, during speciation and divergence (i.e., the evolutionary impact of DNA methylation), remains largely unexplored. For this reason, we first provide a detailed examination of the association between DNA methylation and gene expression across five different tissues. We then identify differentially methylated regions (DMRs) on a genome-wide scale between tissues and species and between parental species and F1 hybrids. We also investigate molecular mechanisms that could explain the differentiation in DNA methylation as well as inheritance patterns of DNA methylation. Finally, we explore the role of genetic vs. epigenetic change in differential gene expression among species as well as in hybrid misexpression.


Study system and sequence data

We performed whole-genome bisulfite sequencing from samples of five tissues: the brain, heart, kidney, liver, and testis in 14 wild-caught male flycatchers belonging to two Ficedula flycatcher species and their naturally occurring F1 hybrids, sampled on the island of Öland in the Baltic Sea (Fig. 1). In total, we sequenced 6 collared flycatchers, 5 pied flycatchers, and 3 F1 hybrids (41.5 billion reads in total) [55]. All F1 hybrids (HYB) were offspring from crosses between female pied (PIE) flycatchers and male collared (COL) flycatchers [56]. Each of the 70 samples was sequenced using 2–9 technical replicates, corresponding to 372 technical samples in total (Additional file 1: Table S1). On average, 592 M reads were obtained per biological sample. We then collated the whole-genome bisulfite sequencing data with previously sequenced RNA-seq data from the same 70 tissue samples [56, 57], thus forming matched genomic data sets.

Fig. 1
figure 1

Distribution map of pied flycatcher (PIE; F. hypoleuca) and collared flycatcher (COL; F. albicollis) in the southern Baltic Sea region. Pied flycatchers are found all over this region (gray) and occur in sympatry (purple) with collared flycatchers in parts of continental Europe as well as two main islands in the Baltic Sea, including Öland where sampling was done for this study

Association between DNA methylation and gene expression across tissues

We assessed the average DNA methylation profile across a set of 8563 genes with annotated 5’ and 3’ UTRs, i.e., defined transcription start sites (TSSs) and transcription termination sites (TTSs) from an updated gene annotation (Additional file 2: Supplementary Results 1). For this purpose, we a priori split the genes into two sets based on presence or absence of CGI annotation in their promoters, referred to as CGI and Other promoters, respectively. Following Mugal et al. [56], we defined promoter regions as the 2-kb upstream region of the TSS. CGI promoters made up 59% of the promoter set. We further split the genes into three categories of gene expression level: low (L: 20% of gene with the lowest gene expression level), high (H: 20% of gene with the highest gene expression level), and medium (M: 60% remaining genes with an intermediate gene expression level).

We computed the DNA methylation profile from 5 kb upstream of the TSS to 5 kb downstream of the TTS separately for tissues and the different sets of genes. Among tissues and irrespective of promoter type, the brain showed the highest average methylation level and the testis the lowest (Fig. 2A–J and Additional file 2: Fig. S1). On average, gene bodies showed higher methylation levels than the 5-kb up- and downstream regions. Splitting the gene body into exons and introns revealed a higher methylation level in exons (Additional file 2: Fig. S2). Genes with CGI promoters showed a drop in DNA methylation levels especially around the TSS (Fig. 2A–E). Genes with low expression generally showed higher promoter methylation, in particular for CGI promoters, but in most tissues lower gene body methylation than the other categories, regardless of promoter type. The exception to the latter pattern was the brain which showed the lowest gene body methylation for genes with the highest amount of expression.

Fig. 2
figure 2

Relationship between DNA methylation and gene expression. Average methylation level was plotted in up- and downstream as well as gene bodies for genes with CGI (AE) and Other (FJ) promoters. Correlation profile between DNA methylation and gene expression across genes (K and L) and across tissue (M). DNA methylation profiles were computed for n = 8563 genes with annotated 5’ and 3’ UTRs

To further assess the relationship between DNA methylation and gene expression, we computed a gene profile of their correlation across genes (Fig. 2K, L). Here, all tissues except the brain showed a weak positive correlation between gene body methylation and gene expression regardless of promoter type. For the brain, the relationship was negative. Correlating gene expression and DNA methylation across genes also revealed a clear difference between promoter types. While higher methylation levels in CGI promoters were consistent with lower expression (Fig. 2K), methylation levels in Other promoters only showed a negative relationship with gene expression for the brain (Fig. 2L).

We also correlated gene expression and DNA methylation for each gene across tissues within groups of samples, which tests whether variation in DNA methylation across tissues is associated with variation in gene expression across tissues (Fig. 2M). Similar to the correlations among genes within each tissue, DNA methylation at CGI promoters were more negatively correlated with gene expression than in Other promoters, but the overall effect was weaker than correlations across genes. A weak negative correlation was observed across tissues for gene body methylation and expression. For these general patterns of gene body and gene-proximal DNA methylation, the difference between COL, hybrids, and PIE was marginal. To conclude, these results show that promoter type is an important predictor for the association between promoter methylation and gene expression and that CGI promoter methylation is associated with tissue-specific gene expression.

Tissue-specific patterns of DNA methylation

We identified regions with tissue-specific patterns of DNA methylation on a genome-wide scale. For this purpose, we called differentially methylated regions (DMRs) between tissues using the BSmooth method [58], separately per species (COL, HYB, and PIE). We then called regions with tissue-specific methylation (tsDMRs) by identifying DMRs unique for a certain focal tissue that is significantly differentially methylated in comparisons to all other tissues. For all sample groups, the testis had the highest number of tsDMRs, followed by the heart in COL and the brain in PIE and hybrids (Fig. 3A–E). In the heart, kidney, and liver, most tsDMRs were hypomethylated while the reverse was true in the testis (Additional file 2: Table S2). We here use hyper- and hypomethylated as relative terms of lower and higher methylation in a comparison, following Hansen et al. [58]. Our findings therefore indicate that tissue-specific methylation can on average either have permissive or repressive functions, dependent on tissue.

Fig. 3
figure 3

Differentially methylated regions between tissues and species. Tissue-specific DMRs (tsDMRs) were shared between sample groups to a much greater extent than expected by chance (odds ratio >1 in all cases; AE). F Enrichment analysis of tsDMRs and annotation features. tsDMRs are especially enriched in CGI promoters and generally (but not always) depleted in transposable elements. G Enrichment analysis of between-species DMRs (spDMRs) and annotation features. spDMRs are also enriched at CGI promoters but generally depleted in CNEEs and introns. Significance in all cases was calculated using 1000 Monte Carlo replicates with a family wise error rate (FWER) of 0.1. In FG, “fixed differences” refer to the 100 ± bp vicinity of fixed differences

It is possible that functionally important tissue-specific patterns are shared between PIE and COL and potentially also present in viable F1 hybrids. For all tissues, more than half of tsDMRs shared between COL and PIE were also found in hybrids (Monte Carlo p-value ≈ 0; Fig. 3A–E). This indicates that the identified tsDMRs represent functionally important genomic regions with constrained methylation levels. Next, we investigated the overlap between tsDMRs and regions with functional annotation (Fig. 3F). A general pattern emerged with tsDMRs enriched in regulatory regions such as 5’ UTRs and promoter regions (especially CGI promoters). Also, conserved non-exonic elements (CNEEs), which are suggested to function as tissue-specific regulatory elements [59], were enriched in tsDMRs in most tissues. Transposable elements (TEs) were underrepresented in all tissues except in the testis, which also showed the greatest bias towards hypermethylation and thus likely repression among their tsDMRs.

We overlapped tsDMRs shared by COL, PIE, and hybrids with genes and their ± 5-kb up- and downstream regions to investigate if tsDMRs show gene ontology (GO) terms consistent with tissue-specific function (Additional files 3–6: Tables S3-6). This revealed that brain-specific DMRs were enriched for genes involved in ion transport, while kidney-specific DMRs were enriched in HOX genes with developmental function. Heart-specific DMRs showed no significant enrichment while testis-specific DMRs were enriched in, e.g., autophagy-related terms but also cardiac cell functions in concordance with most of testis-DMRs being hypermethylated (Additional file 2: Table S2). These findings highlight that the inferred tsDMRs show signatures of tissue-specific repression as well as permission [60]. We next investigated the relationship between tsDMRs overlapping with genes and tissue-specific expression and found evidence for an association in the brain and testis (Additional file 2: Supplementary Results 2). One interesting example of tissue-specific regulation is that of KDM2A, a gene that represses transcription and is involved in pericentromeric chromatin formation through binding unmethylated DNA [61], possibly an important function in the demethylated germline genome. Three conserved tsDMRs are found in the testis inside intron 11 of KDM2A. One of these testis-specific DMRs completely overlaps a CGI, which may act as an alternative promoter for a shorter 5’-truncated forms of KDM2A, which are known from other species [62].

Differential methylation between species is enriched in regulatory regions

We identified DMRs between species (spDMRs) by pairwise comparison between COL, PIE, and hybrids in order to investigate general patterns of DNA methylation differences between the two parental species and their viable F1 hybrids. The brain showed the highest number of spDMRs (11,330) between species while the testis had the lowest (7430) (Table 1). For all tissues, COL were significantly more often hypermethylated compared to PIE (p-value < 2.2 × 10–16; binomial test). In hybrids, the pattern varied among tissues, with lower methylation in the heart and higher methylation in the testis than both parentals. For all tissues, there were more DMRs in the comparison between PIE and hybrids than between COL and hybrids. This means that F1 hybrids have a methylation level closer to COL than PIE. Since all crosses had a COL sire, this could either indicate that the paternal methylation pattern has a greater influence on methylation levels in male offspring or collared-dominant inheritance.

Table 1 Number of spDMRs and the frequency of hypermethylation per comparison. Values were rounded to two decimal points. Differentially methylated regions were assessed genome-wide with COL n = 6, PIE n = 5, and HYB n = 3 samples (except for the kidney were only n = 2 samples were available for HYB)

We next investigated the overlap between genome features and spDMRs. This revealed a significant enrichment of spDMRs in regulatory regions such as promoters and UTRs (Fig. 3G; adjusted Monte Carlo p-value < 0.05, for a sample size of 1000). In contrast, other putative regulatory regions such as CNEEs were significantly underrepresented among spDMRs in most tissues except the testis, in all pairwise comparisons. This indicates that CpG methylation at CNEEs may be under functional constraint. Transposons, in particular, long terminal repeat (LTR) retrotransposons were enriched in spDMRs.

DNA methylation patterns vary more among tissues than species

Our results showed that tsDMRs and spDMRs had similar enrichment patterns in promoters and UTRs but distinct enrichment patterns in CNEEs and TEs (Fig. 3). In order to assess the relative contribution of tissue and species to methylation variation within annotation sets quantitatively, we performed between-groups principal component analysis (BCA) and tested the significance of the amount of explained variation (R 2) using permutations [63]. Overall tissue differences in methylation level were much greater than evolutionary difference between the species (hybrids excluded; Table 2), where TEs showed the lowest between-tissue R 2 (51%) and highest R 2 between species (11%) among annotated features. The ± 100 bp region around fixed differences showed the greatest between-species effect (16%), indicating a cis-genetic coupling between genetic and DNA methylation differences. Annotations with higher between-tissue R 2 had in general lower between-species R 2 even after controlling for between-tissue variation.

Table 2 Between-groups principal component analysis. R 2 is the proportion of variance explained. Values were rounded to two decimal points

Genetic differentiation is correlated with differences in methylation between species

We tested the association between absolute methylation differentiation (M diff) and genetic differentiation (F ST) between COL and PIE across genes (Fig. 4A–C). In general, we observed a clear positive correlation between F ST and M diff of ~ 0.15–0.25 for F ST based on CpG sites in the reconstructed ancestor (CpG F ST) and ~ 0.05–0.1 for F ST based on other sites (non-CpG F ST) (Fig. 4D). There was some variation between tissues, with a stronger correlation coefficient in the kidney and brain and lowest in the liver (Additional file 2: Fig. S3). For CGI promoters, we observed a strong reduction in correlation in the promoter region using both CpG F ST and non-CpG F ST. This region is close to the peak of annotated CGIs along the gene profile (Additional file 2: Fig. S4). Assuming a causal effect of genetic change on DNA methylation differences among species, this could mean that the DNA methylation level at CGI promoters is under purifying selection and that only genetic changes that do not disrupt methylation level are tolerated to segregate at appreciable frequencies. In support of this hypothesis, we do not observe this reduction in correlation for genes with Other promoters, which we previously determined had a weaker relationship with gene expression level and thus epigenetic changes may have less of an effect on expression (Fig. 2L). Furthermore, M diff is three times lower in the promoter region of CGI promoters compared to Other promoters (Fig. 4A). Both non-CpG F ST and CpG F ST are higher in promoters compared to the gene body (Fig. 4B, C), especially the exons, which show lower F ST compatible with purifying selection [64]. This indicates that genetic differentiation in promoter regions is on average tolerated in both CGI and Other promoters but less so when affecting CpG sites in CGI promoters.

Fig. 4
figure 4

Association between genetic and methylation differentiation. A Correlation between F ST and M diff. B Gene profile for M diff. Gene profile for non-CpG (C) and CpG F ST (D). Different lines in (A and D) are different tissues. Different line types represent the average for different promoter types. See Additional file 2: Fig. S3 for (A and D) colored by tissue. Lines in plots represent loess curves with shaded region representing the 95% confidence interval of the local regression. The legend for “Annotation” refers to panels (A–C), while F ST refers to panel (D), and “Promoter type” refers to all panels. For M diff: n = 6 COL and n = 5 PIE samples were used per tissue (generated within this study). For F ST: n = 19 COL and n = 19 PIE individuals were used [65]. M diff and F ST were averaged in bins for n = 8563 genes

Molecular mechanism of DNA methylation differentiation

Changes in DNA methylation that are caused by genetic differentiation between species can either be due to substitutions that affect the local genomic region (cis) or elsewhere in the genome (trans). These molecular mechanism may be distinguished through contrasting allele-specific effects in F1 hybrids with parental differences [36]. For this purpose, we developed a statistical framework applicable to DNA methylation data inspired by previous categorization systems for gene expression [36, 66]. We classified loci into categories based on differences in DNA methylation between the two parental species as well as between the two alleles in the hybrids using a beta regression model and an FDR of 0.1 (Additional file 2: Table S7). To distinguish between alleles in the hybrids, we polarized bisulfite-seq reads based on fixed differences between COL and PIE and calculated allele-specific methylation. In total, around 37,000 fixed differences per tissue were used as markers (Additional file 2: Table S8). We measured the allele-specific methylation in the ± 100 bp region from a fixed difference. Most loci either had no CpG site or too low coverage to be considered. This restricted the dataset to 3251, 5939, 6120, and 7701 loci in the brain, heart, liver, and testis, respectively.

In our analysis, we focused on the relative contribution of cis and trans effects on DNA methylation and if this varies with tissues and sequence divergence. A prediction based on gene expression differences in flies and yeast is that the proportion of cis-changes should increase with sequence divergence [67, 68]. Since the Z sex chromosome is more genetically differentiated between flycatchers than autosomes [69], we performed the analysis separately for Z and autosomes. In total, 1.6% of all callable loci across tissues showed a statistically detectable difference in at least one of the pairwise comparisons COL/PIE or the two alleles in the hybrids (Fig. 5A). Distribution of regulatory patterns were significantly different among tissues (Fig. 5; Fisher’s exact test, p ≈ 0.0002). The brain showed an excess of conserved loci (Fig. 5B). Among divergent (non-conserved) loci, there was also a difference between tissues (p ≈ 0.01). An excess of trans differences was found on brain autosomes (Binomial test, p ≈ 0.001). For other tissues and chromosome types, an even contribution of cis and trans differences could not be rejected. Overall, distribution of divergence categories was not significantly different between autosomes and the Z chromosome (p > 0.05). The exception was the testis (p ≈ 0.001), which showed an excess of non-conserved loci on the Z compared to autosomes. A trend with a larger proportion of cis relative to trans on Z compared to autosomes were observed (Fig. 5B) but was not significant for any tissue (Fisher’s exact test, p > 0.05).

Fig. 5
figure 5

Molecular mechanism of DNA methylation differentiation. A Proportion of different classes of molecular mechanisms of DNA methylation differentiation at callable loci. B Proportion of divergent (non-conserved) classes per tissue and chromosome type. C–F Difference between COL and PIE alleles in their native parental and hybrid cellular environment. Loci with cis-only effects are clustered close to the origin as predicted when the methylation level is independent of cellular environment. In contrast, loci with trans effects show greater variation but the general trend is a spread along the additive inheritance axis. Bold black lines in (C–F) are major axis regression lines. Values of p in (C–F) are rounded to two decimals. Around 37,000 fixed differences were used as markers with the exact number dependent on tissue (Additional file 2: Table S8). At each marker loci, the sample size was n = 3 for each sample group

We next determined the inheritance pattern of methylation at these loci and investigated its relationship with allele-specific methylation in hybrids. This revealed that the inheritance pattern is tightly related to the mechanism of molecular divergence (Fig. 5C–F). Major axis regressions between parental and hybrid differences were negative for all tissues (p < 0.05), indicating an overall pattern of additivity (Fig. 5C–F). Subtle differences in spread are visible among tissues, with the liver having the largest variation between PIE in parental and in hybrids indicating more COL dominance (Fig. 5E). The results illustrate two different mechanisms of additive inheritance of a molecular phenotype. Loci with methylation difference due only to cis effects showed essentially no difference in methylation level whether in parentals or in hybrids which indicate that they are strictly additively inherited. In other words, cis loci are unaffected by the hybrid cellular environment. Other loci are affected by the hybrid cellular environment unveiling a trans effect. That trans difference here generally makes the alleles in hybrids more similar to each other which creates a different route to additive inheritance since the overall state is in-between parentals. These two types of additive inheritance are only distinguishable when measuring allele-specific methylation.

Tissue-specific association of differential expression with genetic- and epigenetic change between collared- and pied flycatchers

Much of the interest in DNA methylation lies in its ability to regulate gene expression. Nevertheless, it is not clear to what extent changes in DNA methylation are involved in the evolution of differential expression (DE) between species. Besides DNA methylation patterns, also genetic changes in promoter regions are expected to affect transcriptional regulation and could ultimately result in differential expression. We sought to test the relative importance of changes in promoter methylation and genetic change in the evolution of differential expression between COL and PIE. Here, changes in promoter methylation are used to assess epigenetic change. We used non-CpG F ST to assess genetic changes to more clearly separate genetic and epigenetic effects. First, we explored the patterns of differentiation around well-annotated genes (having both 5’ and 3’ UTR). The number of differentially expressed (FDR = 0.1) genes between COL and PIE were lowest in the brain and highest in the testis (Additional file 2: Table S9). As reported above (Fig. 3G), between-species DMRs are enriched in the CGI promoter region, with an elevated enrichment in DE compared to non-DE genes (Table 3, Fig. 6A). Notably, this same effect is not observed in Other promoters (Table 3, Fig. 6C). In addition, we observed no clear difference in DMR frequency in the bodies of DE vs. non-DE genes for most tissues. Genetic differentiation (non-CpG F ST) is higher for DE genes in CGI promoters though only significant in the brain and the testis after multiple-test correction (Table 3, Fig. 6B). All tissue and promoter type combinations had significantly higher non-CpG F ST in the bodies of DE vs. non-DE genes (Table 3, Figs. 4C and 6B), revealing a strong association between local genetic differentiation and differential expression in this system.

Table 3 Determinants in cis of DE between COL and PIE. Average DMR frequency differences and non-CpG F ST in the 2 k upstream promoter region and throughout the gene body were compared between DE and non-DE genes using paired t-tests. The table displays p-values of those tests
Fig. 6
figure 6

Patterns of genetic and epigenetic change at differentially expressed genes between COL and PIE. DE genes had more DMRs in CGI promoters (A) but not in the Other promoters (C). DE genes also had higher non-CpG F ST (B) across the gene bodies in all tissues and for both promoter types (B and D). Lines in plots represent loess curves and shaded regions around lines are the 95% confidence intervals. Differential expression was assessed at n = 8418 of the 8563 well-annotated genes. For each tissue, the sample size was n = 6, n = 5, and n = 3 for COL, PIE, and hybrids, respectively

Cis-genetic and -epigenetic changes do not predict misexpression in hybrids

We also investigated whether genetic differentiation between COL and PIE as well as hybrid-specific methylation changes could predict misexpression in hybrids, i.e., DE between hybrids and both parental species in the same direction. In other words, an overdominant or underdominant expression pattern. Overall inheritance of methylation patterns for both promoter types were mainly additive (Additional file 2: Supplementary Results 3), in concordance with results from methylation differentiation around fixed differences (Fig. 5). We observed less mismethylation in CGI promoters compared to Other promoters (Additional file 2: Supplementary Results 3). In general, neither DMRs between each parental species and hybrids nor non-CpG F ST in promoter regions were higher for DE genes between hybrids and both parental species (Additional file 2: Fig. S5). In gene bodies, non-CpG F ST was still higher in misexpressed genes in some cases, while significantly lower in the heart (Additional file 2: Table S10). These results indicate that both cis-genetic and cis-methylation differences play a relatively minor part in F1 hybrid misexpression patterns and are perhaps dwarfed in importance by trans effects or distal cis effects.


In this study, we investigated both the regulatory and evolutionary impact of DNA methylation in two species of Ficedula flycatchers and their F1 hybrids. Using these two layers of analyses, we could evaluate both the functional impact of DNA methylation in two wild bird species and its relation to genetic divergence and differential gene expression. In addition, by using F1 hybrids, we could investigate the role of methylation in hybrid misexpression, thus probing deeper into the regulatory underpinnings of hybrid dysfunction.

Genome-wide we observed much greater differentiation of DNA methylation between tissues than species, as expected. A similar pattern was observed in a study of three primates including humans [70]. Among the tissues studied here, the brain showed the most unique methylation profile. Only the brain tissue had a negative correlation between gene expression and gene body methylation. This distinctive pattern could be driven by the special role of the MeCP2 gene, which is expressed at remarkable levels of more than > 16 M proteins per nucleus in nerve cells of mammals [71]. MeCP2 binds to methylated C’s and recruits the NCoR1/2 co-repressor complex causing transcriptional repression possibly through deacetylating histone tails [21, 72]. This mechanism could potentially make DNA methylation marks at gene bodies a target of chromatin compaction, which could explain why gene body methylation was negatively correlated with expression in this case.

We observed that the relationship between DNA methylation at promoters and gene expression differed among promoters with or without CGIs, which we suggest affect their different evolutionary constraints. DNA methylation level at CGI promoters showed a stronger relationship with gene expression, in congruence with a recent study of fibroblasts from six mammals and chicken [54]. We cannot rule out that some genes with Other promoters are also regulated by DNA methylation as have been observed in some systems [23, 33]. Our study highlights the benefits of whole-genome bisulfite sequencing for gaining a complete view of the importance of DNA methylation, in contrast to the popular reduced-representation bisulfite sequencing method which may bias analysis to CpG-rich genes and promoters where variation in DNA methylation may be more impactful [73, 74]. For example, we observed that CGI promoter methylation to a higher degree was associated with tissue-specific expression compared to methylation patterns at Other promoters, a conclusion which would be difficult to draw using reduced-representation data generated with CG-specific digestion enzymes.

Beyond promoters and gene bodies, less is known of the regulatory impact of differential methylation. What is for example the impact of DNA methylation levels on enhancer or other kinds of cis-regulatory sequences? We observed an enrichment of tsDMRs and a lack of spDMRs in CNEEs, a subset of which putatively acts as tissue-specific cis-regulatory elements. DNA methylation at cis-regulatory elements impacts binding affinities both positively and negatively, but the analysis is complicated by the fact that transcription factor binding itself seems to induce active demethylation possibly through attracting TET demethylases [20, 75]. Still, we can conclude that flycatcher CNEEs are likely to be involved in tissue-specific regulation of expression potentially through tissue-specific methylation, though experimental studies would be needed for robust confirmation.

While DNA methylation was constrained among species at CNEEs, the reverse was observed for TEs. They showed of a lack of tsDMRs and an enrichment of spDMRs. More pronounced interspecific differences in DNA methylation at TEs than at CNEEs have previously been observed in primates [76]. TEs could be showing relatively more variation among species due to a higher turnover of CpG sites resulting from high methylation levels or relaxed purifying selection on methylation level. In addition, active copies of LTR retrotransposons carry regulatory elements which could potentially have a functional role in regulating gene expression [77]. TEs may affect divergence and speciation of lineages in more ways than through their associated regulatory elements affecting expression levels of proximal genes. In hybrids, TEs can be misregulated, which could lead to hybrid dysfunctions or sterility as has been observed for the P-element in Drosophila [78]. While beyond the scope of this paper, the presented dataset could be used to test this hypothesis. With that said, in flycatchers, it is possible that genic effects are causing hybrid male sterility since dysregulation of meiotic genes have been observed through single-cell RNA sequencing of male testes [79].

Despite having a relationship with gene expression and thus stronger functional constraints, we did observe an enrichment of spDMRs in CGI promoters. When examining the gene profile of DMR frequency, we find that the peak enrichment of spDMRs in the promoter is in most cases distal to the TSS within the promoter. Therefore, the enrichment of spDMRs in promoters is not coinciding with the minima of correlation between F ST and M diff which is more proximal to the TSS, rich in CGI annotation (c.f. Figs. 4D and 6A, and Additional file 2: Fig. S4). Nevertheless, since we observed an enrichment of spDMRs within CGI promoters but a marked drop in cis-genetic correlation, it is likely that some spDMRs at CGI promoters in this system are driven by trans differences in methylation. A non-exclusive explanation would be that DMRs are concentrated to the shores of CGIs, which have been shown to be especially prone to differential methylation [80]. While we could survey genome-wide patterns of cis and trans differences in methylation using fixed differences as marker loci, we lack power to distinguish their relative contribution to differentiation at promoters. Ideally, the F1 hybrid method of determining the molecular mechanism of DNA methylation differentiation is complemented by sequencing of trios (parents and hybrid offspring), which would enable much more and less biased marker loci to be used. Even so, we did observe a stronger species-effect in DNA methylation at fixed differences compared to other annotation sets further indicating that genetic and DNA methylation change are coupled [34]. In spite of this, tissue was a stronger determinant of DNA methylation level than species also at fixed differences, highlighting the importance of studying several organs/tissues or cell types for understanding the evolution of DNA methylation [81].

For loci with fixed sequence differences between the two species, we observed an equal amount of methylation differentiation attributable to cis or trans effects for all tissues except the brain where trans differences dominated. In a cross between red-jungle fowl and domestic white leghorn chickens, a greater share of cis differences was observed in a quantitative trait locus analysis of hypothalamus tissue [51]. This could be due to tissue- or species-specific patterns or the ~ 3 × greater genetic differentiation between red-jungle fowl and white leghorn compared to COL and PIE. A greater share of cis difference in gene expression with greater sequence divergence has been observed in fruit flies and yeast [67, 68]. Though insignificant, we here observed the same trend for DNA methylation in the contrast between Z and autosomes, where the more differentiated Z has a trend of higher share of cis. By contrasting the allele effect in hybrids with parental species, we also observed that trans differences generally where additively inherited but through a distinct mechanism where both HYB alleles converged to the midparental value of DNA methylation. This could for example be caused by the relatively unexplored molecular process of transvection [82, 83], in which alleles affect each other’s phenotype, which have previously been observed to affect DNA methylation patterns during meiosis in mice [82]. However, since we only used F1 hybrids from one direction of the cross, inheritance patterns could to some extent be confounded with parental effects.

In the last two decades, many studies have presented evidence of misexpression in F1 hybrids of a wide range of species [14, 56, 84,85,86]. While some patterns are emerging, the regulatory mechanisms underlying hybrid misexpression remains relatively obscure. Dependent on tissue, both higher DMR frequency and larger F ST in promoters were associated with DE between COL and PIE, while the same variables failed to predict DE between parentals and hybrids. If hybrid misexpression was driven by evolution at proximal cis-regulatory elements of many genes, we would have expected to find similar patterns in both COL-PIE and parental-hybrid comparisons. It has been suggested that incompatibilities between divergent cis-regulatory elements and trans-acting factors in hybrids result in misregulation of genes [11, 13]. Misregulation evolves quickest when interacting cis- and trans factors diverge under positive selection [8], which may be the case for a subset of the genes differentially expressed between COL and PIE [87]. However, expression is likely in many cases to be under stabilizing selection [88, 89], which may also cause BDMIs if different compensatory mutations fix in diverging lineages [84]. We did observe compensatory evolution of DNA methylation at fixed differences, though most were additively inherited and consequently not mismethylated in hybrids. In addition, we did not observe greater frequencies of DMRs at promoters of misexpressed DE genes which would be expected in a model where hybrid over- or underdominance in gene expression is caused by loss or gain of promoter methylation. One caveat is that such mismethylation might be so deleterious that we miss it when sampling adult yearlings that survived migration to Africa and back and thus are likely to be in better condition than the average natural F1 hybrid hatched [43]. In theory, genetic or epigenetic misregulation at a single or few two-locus interactions could potentially cascade throughout the expression network and cause hybrid misexpression at many genes. Misregulation of upstream cis-trans interactions could overshadow the cis effects of methylation and genetic difference associated with DE in the parental comparison. In this model, most misexpressed genes do not have incompatibilities in their promoter regions, instead they are symptoms of a few rare cascading interactions. It is possible that the hybrid misexpression observed here are metabolic responses possibly related to the transgressive metabolic rate observed in F1 hybrid flycatchers [44, 45] and somewhat in line with results in a copepod with known mitochondrial dysfunction in F2 hybrids [90].

Our results showed that there was a weak but pervasive correlation between non-CpG genetic differentiation and DNA methylation except proximal to the TSS (for CGI promoters). While it is possible that a genetic change affects, e.g., the binding of a transcription factor leading indirectly to a change in methylation, it is harder to conceive of a molecular mechanism supporting the other direction of causation though it cannot be entirely disregarded. Recent empirical and theoretical studies suggest that epigenetic variation could act as a first substrate for divergent selection and may promote or slow down speciation [91,92,93]. In birds, epigenetic effects independent of genetic effects may have a limited impact on speciation processes due to weak evidence for genomic imprinting [94, 95]. However, transgenerational effects have been observed in several species [reviewed in 95] and thus research is needed to understand the molecular mechanism of transgenerational effects in birds and potential associations with reproductive isolation. With this in mind, a limited transgenerational inheritance or extent of imprinting does not mean that epigenetic mechanisms are unimportant for speciation in birds and other vertebrates. Epigenetic mechanisms are fundamental in conserving and reshaping transcriptional states, as illustrated by the prominent role of DNA methylation in eye degeneration of the cave morph of the Mexican tetra [96]. In other words, epigenetic mechanisms can play important parts without being the ultimate cause in cases of adaptive differentiation and hybrid dysfunction, both of which are important aspects of speciation.


In this study, we investigated the relationship between genetic differentiation, DNA methylation, and gene expression in two songbird species and their F1 hybrids. The hybrids were used both as a tool to investigate the relative contribution of cis- and trans-factors in DNA methylation differentiation between the parental species and to investigate the role of DNA methylation in hybrid misexpression. We showed that DNA methylation matters more for both general gene expression patterns and differential expression in genes with CGI promoters. In general, genetic differentiation predicted methylation differentiation but this relationship broke down close to the transcription start site where methylation patterns matter the most for gene regulation. This indicates that DNA methylation levels at CGI promoters are under strong purifying selection and remains conserved among species. While these results suggest that DNA methylation could be important in hybrid misexpression, we did not find evidence for this, indicating that trans factors are more important than cis factors for hybrid misexpression. Overall, this study provides a deeper understanding of the evolution of DNA methylation patterns in songbirds and their role in gene regulation.


Sampling scheme and tissue collection

Male flycatchers were collected at the Baltic Island of Öland (57°10’N, 16°56’E) during the breeding season of 2014 [56]. Collected samples included six collared (Ficedula albicollis), five pied (F. hypoleuca), and three F1 hybrid flycatchers (♀PIE x ♂COL). Initially, four individuals were classified as hybrids using plumage score [40]. Later three of those individuals were confirmed genetically as hybrids [56], while a fourth was identified as a collared flycatcher (in this study included within the six collared flycatcher individuals). Tissue collection is described in more detail in Mugal et al. [56]. All sampling procedures were approved by the Swedish Board of Agriculture (Jordbruksverket—DNR 21-11). The following organs/tissues were used for RNA sequencing [56] and whole-genome bisulfite sequencing (WGBS; within the present study): brain (caudal region of the telencephalon) heart, kidney, liver, and testis. In addition, we included RNA sequencing data for spleen for the collared flycatcher individuals.

Nucleic acid extractions and sequencing

Samples were homogenized using a bead beater with ceramic beads and aliquots were used to extract DNA and RNA. For details on RNA extraction and sequencing, see [56]. DNA was extracted using the phenol-chloroform method. Library preparation and whole-genome bisulfite sequencing were performed by the SciLife SNP&SEQ Technology Platform in Uppsala, Sweden. Sequencing libraries were prepared from 100 ng of DNA using the TruSeq (EpiGnome) Methylation kit (Illumina Inc., EGMK91324) according to the manufacturers’ protocol (#15066014). Samples were multiplexed and split into several lanes as well as separate sequencing runs (Additional file 1: Table S1). Sequencing was performed in two separate sequencing efforts, split into a pilot study consisting of brain samples, and a second effort consisting of the remaining samples. Brain samples were sequenced using v4 sequencing chemistry and 125 bp paired-end reads on the Illumina HiSeq2500. In total, 14 lanes were used with two technical replicates and one biological sample per lane (Additional file 1: Table S1). The rest of samples were sequenced using v2.5 sequencing chemistry HiSeqX with 150 bp paired-end reads. In total, 87 lanes were used with 6–9 technical replicates per biological sample. One kidney sample from one F1 hybrid flycatcher (HYB02) was discarded due to allelic imbalance (data not shown).

Processing of bisulfite sequence reads and methylation calls

Quality control, filtering, and mapping of reads as well as methylation calls were performed using the reproducible Nextflow workflow v20.10.0 nf-core Methylseq v1.5 [97]. Raw reads were quality controlled using FastQC [98]. Adapter sequences were removed using Trim Galore! v0.6.4_dev, a wrapper for Cutadapt v2.9 [99]. Further clipping was performed according to the Epignome profile (8 bp from both 5’ and 3’ ends of both reads in a pair). Reads were aligned to the chromosome version of the collared flycatcher reference genome, FicAlb1.5 [100] using Bismark v0.22.3 with Bowtie2 as alignment tool [101]. Alignment quality was assessed using Qualimap v2.2.2-dev [102] and visualized using MultiQC v1.8 [103]. The percentage of uniquely mapped reads ranged from 50.2 to 70.1% (Additional file 1: Table S1). Fixed differences between collared and pied flycatchers—determined using 19 individuals of each species previously sampled on Öland [65]—were masked prior to read-mapping. Bismark v0.22.3 was then used to deduplicate alignments and extract methylation calls for CpG sites. After read mapping and deduplication, the median coverage ranged from 7 × to 39 × .

Methylation level

Since the tissue samples we used in this study consisted of a population of cells each with a possibility of either having a methyl mark or not at a certain CpG position, we assessed the methylation status at each CpG position as the proportion of methylated reads. To measure methylation level, we summarized the number of methylated reads (x mCpG) mapping to both strands of a reference CpG dinucleotide and divided by the total number of reads: x uCpG + x mCpG. As a measure of methylation level over n CpG dinucleotides, within a defined region, we used the average proportion across the individual dinucleotides (i),

$$Methylation \,level= \frac{\sum_{i=1}^{n}\frac{{x}_{mCpG}}{{x}_{mCpG}+{x}_{uCpG}}}{n}.$$

We only included dinucleotides with at least 6 and at most 200 mapped reads unless otherwise stated. The lower limit is included to reduce bias in estimating methylation level caused by lack of data. The upper limit is included to reduce bias induced by collapsed repetitive regions. Other methods and filtering thresholds were used for the smoothed methylation values produced by BSmooth for DMR analysis and BiSeq for the analysis of molecular mechanism of DNA methylation differentiation respectively (see below).

Transcriptome assembly

To build the transcriptome assembly, we used RNA-seq reads across all five tissues plus spleen from the six collared flycatchers described above (in total 36 samples). We concatenated all samples and generated four separate and one consensus de novo transcriptome assembly using the Oyster-River protocol v2.3.1 [104]. Adapters were removed and bases with a Phred score lower than 36 were trimmed from the ends of reads using Trimmomatic v0.39-1 [105]. Reads were further error-corrected using Rcorrector v1.0.3 [106]. Two transcriptome assemblies were created using Spades RNAseq assembler v3.14 with k-mer size 55 and 75, respectively [107]. Further transcriptome assemblies were created with Trans-ABySS v2.0.1 [108] and Trinity v2.9.1 [109]. Assemblies were merged using OrthoFuse [104]. First the separate assemblies were concatenated and groups of transcripts were identified using a modified version of OrthoFinder [110]. The best transcript in each group was then found based on contig score using a modified version of TransRate [111]. Transcripts with less than 1 TPM in the concatenated set of RNA samples and no hit to the Swissprot database were removed from the consensus assembly. The consensus assembly had a TransRate score of 0.1441 and 96.6% of all BUSCO v3.0.2 genes using the aves_odb9 lineage dataset [112]. Optimal TransRate score for the assembly was 0.2341 but filtering away lower quality contigs reduced the BUSCO coverage to 94%. Since we here used the transcriptome assembly as a basis for a gene annotation update, we decided to use the more complete assembly of slightly lower contig quality.

Gene annotation update

With a main goal of improving the annotation of transcription start sites (TSS) and transcription termination sites (TTS) of genes by annotation of their untranslated regions (UTR), we updated the gene annotation for the collared flycatcher genome assembly based on a de novo transcriptome assembly (see above). For building gene models, we used MAKER v2.31.10 [113]. We configured MAKER to update the gene models from the collared flycatcher Ensembl annotation (v96), by setting this annotation as pred_gff in the maker_opts.ctl file. Collared flycatcher RNA-seq evidence for the Ensembl annotation consisted of eight adult organs/tissues as well as embryo [114]. As additional evidence in the annotation, we used the Oyster-River protocol transcriptome assembly and proteins from the chicken (Gallus sea, Ensembl annotation v98) and zebra finch (Taeniopygia guttata, Ensembl annotation v102). We configured MAKER to allow gene models to be built directly from transcripts and protein homology. We also included genes found using a pipeline designed to identify so-called missing genes, i.e., genes that have proven difficult to annotate in the bird genome because of repetitiveness or extreme base composition such as high GC content [115]. Using candidate proteins (n = 2454) from the Chinese softshell turtle (Pelodiscus sinensis, Ensembl annotation v98), we found collared flycatcher candidate hits (n = 1389) using tBLASTn v2.7.1+ to an earlier version of the transcriptome assembly [116]. Candidates were converted to gff3 file format and included as predictions in MAKER (pred_gff option). We used 20,000 bp as the expected max intron size for evidence alignments and conservatively did not consider single exon transcript evidence when generating annotations. We also used an updated repeat annotation consisting of Aves repeats from RepeatMasker v4.0.7_Perl5.24.1, which was mainly repeats curated from chicken and zebra finch, as well as repeats from collared flycatcher [117], hooded crow [118], blue-capped cordon-bleu [119], paradise crow (Lycocorax pyrrhopterus, [120], Huon astrapia (Astrapia rothschildi), and paradise riflebird (Ptiloris paradiseus) [121]. Including consensus sequences derived from TEs in related species has been shown to improve detection of repeats missed by species-specific repeat libraries [119] and should consequently decrease the risk of including TE genes in the gene annotation.

RNA sequence read analysis

RNA-seq reads were mapped to the collared flycatcher reference genome, FicAlb1.5 [100], and the updated gene annotation with fixed differences between collared and pied flycatchers masked. All steps from quality control, read mapping to differential expression analysis were performed using Nextflow v21.02.0.edge nf-core rnaseq v3.0 pipeline [97]. Quality of raw reads were assessed with FastQC v0.11.9 [98] and adapters were removed using Trim Galore! v0.6.6, a wrapper for Cutadapt v2.10 [99]. Bases with a Phred score < 20 were removed from the 3’ end of reads. Trimmed reads were aligned to the reference genome (see above) using STAR v2.6.1d [122]. Transcript quantification was performed using Salmon v1.4.0 [123]. Quality control was done using RSeQC v3.0.1 [124], SAMtools v1.1.0 [125], and visualized with MultiQC v1.9 [103]. Differential expression analysis was done with DeSeq2 v1.28 [126] with Salmon count data imported using tximport [127]. We considered genes with an FDR-adjusted p-value < 0.1 as differentially expressed.

Additional annotation tracks

We defined promoters as the 2-kb upstream region of the TSS for the set of 9597 genes with at least one transcript with annotated 5’ UTR. CpG islands (CGIs) were inferred using CpGCLUSTER v1.0, with default parameter settings and a minimum length of at least 50 bp [128]. Promoters intersecting CGIs were classified as CGI promoters using BEDtools v2.29.2 [129]. Promoters without any overlapping CGIs were defined as having other types of promoters (Other). Phylogenetic conserved elements (CEs) based on PhastCons [130] and a whole-genome alignment of 23 sauropsids were retrieved from [131]. Conserved non-exonic elements (CNEEs) were defined as CEs which did not overlap exons.

Gene profile

We investigated the average methylation patterns at genes and their up- and downstream regions (gene profile). To construct a gene profile, we filtered for genes having at least one transcript with 5’ and one with a 3’UTR such that TSS and TTS were defined. Upstream and downstream 5-kb windows of each gene were split into 100 bp nonoverlapping segments. Genes were also split into 100-bp segments, then averaged per 99 ranks across the gene length. Variables of interest, such as methylation level, were averaged across genes in each of these 100-bp segments. For correlation gene profiles, we calculated the Spearman rank correlation coefficient between two selected variables in a segment across all genes using the cor.test function in R v4.0.4 [132].

Identification of differentially methylated regions

Differentially methylated regions (DMRs) between samples were identified using the BSmooth method [58]. We defined DMRs as regions in the 0.01 and 0.99 quantiles of methylation difference, keeping only CpG sites where at least two samples per group (of the pairwise comparison) have a coverage ≥ 2. In addition, DMRs needed to span at least 3 CpGs with a mean methylation difference equal to or larger than 0.1. Significant deviation from random (none of the compared groups had more DMRs with higher methylation, i.e., the random expectation is 0.5) in the either direction (hypo- or hypermethylation) for a set of DMRs was determined using binomial tests in R v4.0.4 [132].

Classification of DMRs

We called DMRs between both tissues and species (spDMRs). We defined tissue-specific DMRs (tsDMRs) as regions with differential methylation in a focal tissue compared to all other tissues using BEDtools v2.29.2 intersect requiring at least 25% reciprocal overlap [70]. Furthermore, to fulfill the criterion of tissue specificity, the same region could not be classified as DMR in any other tissue comparison.

A method for enrichment analysis between two sets of genomic ranges using resampling

Enrichment analyses of various classes of DMRs in genomic annotation tracks were performed using a custom Bash script employing BEDtools v2.29.2 intersect and calculation of empirical p-values using a Monte Carlo randomization procedure. For a certain overlap between an annotation track (e.g., introns) and a set of DMRs (e.g., between COL and PIE heart samples), we shuffled the DMRs 1000 times across the genome and calculated the total number of base pairs that overlapped per resampling replicate. To calculate the empirical p-value, we compared the overlap in the resampling replicates with that of the actual data and calculated the p-value as r/n, where r is the number of replicates with an overlap greater than or equal to the overlap for actual data [133]. p-values were corrected using the Bonferroni method to a family-wise error rate of 0.1. Enrichment was defined as the following odds ratio:

$$\textit{Odds ratio}=\frac{\textit{Overlap}(\textit{DMRs and annotation track})}{\textit{Total length of DMRs}}/\frac{\textit{Total length of annotation track}}{\textit{Genome length}}.$$

Gene ontology analysis

We performed GO analysis of tsDMRs shared among COL, PIE, and hybrids overlapping genes and their ± 5-kb neighboring regions, using ShinyGO v0.77 [134] with GO biological process as database. We used an FDR of 0.1 and collared flycatcher as reference annotation.

Tissue specificity of methylation and gene expression

Tissue-specific gene expression was calculated using the preferential expression measure (PEM, [135]) based on the five tissues represented in the study. Genes with no expression in any of the five tissues were excluded. PEM is a relative measure which ranges from − ∞ to 1, with values below and above 0 representing relative under- and overexpression, respectively. To check for an association between tissue-specific methylation and expression, we selected tsDMRs overlapping annotated promoters. We then tested for a significant deviation in PEM rank for genes with a tsDMR in a focal tissue using χ 2 test of independence (chisq.test in R). If there is a relationship between tsDMRs and PEM rank, we expect to see the focal tissue with tsDMR in promoter being either the most or least expressed gene among tissues. For tissues and species combinations (e.g., COL × testis) that rejected the null hypothesis of the first test, we also proceeded with a secondary test. We tested for a difference in mean PEM between genes with a promoter tsDMR in a focal tissue versus genes with no promoter tsDMR in any tissue (reference set) using the non-parametric Wilcoxon test. Here, we tested for differences between overexpressed (PEM > 0) and underexpressed (PEM < 0) genes separately. This test investigates the relative importance of DNA methylation in tissue-specific expression since the reference set consists both of genes with ubiquitous expression, as well as tissue-specific genes controlled by mechanisms other than tissue-specific DNA methylation.

Classification framework for hybrid inheritance patterns at promoters

We classified the methylation patterns of F1 hybrids using a cutoff strategy into the following categories of inheritance pattern: conserved, additive, collared-dominant, pied-dominant, and the mismethylation categories overdominant and underdominant. Additive means a hybrid methylation level between parentals. Collared-dominant corresponds to a methylation level in hybrids close to the value of COL and vice-versa for pied-dominant. Over- and underdominant means hybrid methylation level above and below the value of both COL and PIE, respectively. We used a cutoff strategy to determine inheritance pattern instead of statistical tests because different number of significant tests would be needed for different categories. For example, no significant difference between F1 and both parentals for conserved but significantly greater than both parentals for overdominant [86]. However, a simple cutoff strategy is also biased. Consider that we measure the phenotype P of two parental species A and B and their F1 hybrids H. If we assume that the phenotypic value P is pairwise-independent between A, B, and H then,

$$\textit{cov}\left(P_{\mathit H}-P_{\mathit A},P_{\mathit H}-P_{\mathit B}\right)=\textit{var}\left(P_{\mathit H}\right).$$

Since this variance is almost certainly not zero in all practical cases, applying a cutoff will inflate the relative proportion of over- and underdominance. To deal with this artifact, for each promoter, we randomly picked a hybrid sample for P H in P HP A and another for P H in P HP B. In effect, this reduces the correlation in error between X and Y axis [136]. Ideally, with a large sample size, many samples would be picked for each group. Simulations showed that this gives a roughly rhombic error profile around the origin, regardless of whether the measured phenotype is a uniform or Poisson random variable (data not shown). To make further classifications as fair as possible, we split the X, Y field into 8 areas defined by slices of π/8 radians. To accommodate this fairer classification scheme, we then defined a circular cutoff of 0.1 to classify promoters as either conserved (< 0.1) or not (> 0.1).

Ancestral genome reconstruction and estimates of genetic differentiation

To reconstruct an ancestral genome for the black-and-white flycatchers, we used all-sites genotype data mapped to the COL genome assembly from one red-breasted flycatcher (Ficedula parva) and one snowy browed flycatcher (F. hyperythra) previously sequenced [137]. Sites for which red-breasted and snowy browed flycatcher were monomorphic were considered callable with the ancestral state being that allele. Sites polymorphic within or between these two species were ignored. To estimate genetic differentiation, we used SNPs from 19 COL and 19 PIE flycatchers previously sampled on Öland [65]. Genetic differentiation between COL and PIE was estimated separately for ancestral CpG sites (CpG) and other contexts (non-CpG) using the fixation index F ST [138] implemented in vcftools v0.1.16 [139].

Allele-specific methylation estimation

We estimated allele-specific methylation in F1 hybrids using fixed differences (see above) between collared and pied flycatchers as markers within the bisulfite sequencing reads. First, deduplicated.bam files of bisulfite sequence reads produced by bismark v0.22.3 were split according to parent-of-origin allele using SNPsplit v0.3.2 [140]. Fixed differences C and T for forward strand alignments and G and A for reverse strand alignments were ignored since they are indistinguishable from the pattern produced by the bisulfite treatment. Allelic imbalance in the number of allele-specific reads were determined with SNPsplit and the kidney sample for HYB02 was excluded due to extreme allelic imbalance in bisulfite-seq reads but not RNA-seq reads (which were analyzed to determine whether the imbalance of the bisulfite-seq reads in HYB02 kidney was a biological effect or a technical artifact). To ensure a sample size of at least three, the kidney was not considered further in allele-specific analyses. Allele-specific methylation was called using bismark v0.22.1 methylation extractor and the methylation level was measured in 200-bp windows centered at the fixed difference. To measure methylation difference between samples, we used the R package BiSeq v1.28.0 [141]. In F1 hybrids, read coverage was limited to the 0.9 quantile. Parental species for each tissue were randomly downsampled to 3 individuals and coverage was capped at the 0.45 quantile to mimic the sample size and read coverage of hybrids. Difference in methylation between groups of samples for each 200-bp locus were determined using the beta regression model in BiSeq v1.28.0 [141].

Statistical framework for molecular mechanism of DNA methylation differentiation

We classified the molecular mechanism of DNA methylation differentiation by comparing the methylation difference between the parental species and between the parental alleles in the F1 hybrid environment [36, 66]. To ensure that the same number of null hypotheses needed to be rejected for calling cis and trans differences, we also compared the methylation difference of parental alleles and hybrid alleles of the same origin (e.g., PIE: parental—hybrid; Additional file 2: Table S7). This constitutes a trans effect test since by definition, we expect the same allele to have different methylation levels between the parent and hybrid environment, if such an effect occurs. By comparing both COL and PIE methylation differences in hybrids and the trans effect test, we tested both row and column null hypotheses of a 2 × 2 matrix with COL, PIE, parental, and hybrid as column and row names, respectively. Due to the dependence of tests in this approach which is a feature it shares with the original framework [136], column and row tests should ideally be done using different samples, though that would require a sample size of hybrids of at least 6, in this case. If any sample group lacked methylation read information for a locus, then that locus was classified as ambiguous. We classified a locus as conserved if there was no difference between COL and PIE in either the parental or hybrid comparison. For a cis-only change, the COL and PIE allele at a locus had to be significantly different in the same direction both within hybrid and between parentals and no significant trans effect. Three different outcomes were possible if both a cis and a trans effect were acting at a locus: (1) cis + trans, significant difference between parental and hybrid allele for either COL or PIE but with a trans effect in the same direction as the cis effect, (2) cis x trans, cis, and trans effect in opposing directions, and (3) compensatory, no significant difference in the parental comparison while significant difference in hybrids and a trans effect. For trans-only, there needed to be a significant difference between parentals but not between alleles within F1 hybrids and a significant trans effect. To estimate inheritance patterns at fixed difference marker loci, we performed major axis regression (a type of model II regression that does not assume that the predictor variable is fixed) using COL:Parental-Hybrid as response variable and PIE:Parental-Hybrid as predictor variable.

Availability of data and materials

All data generated or analyzed during this study are included in this published article, its supplementary information files, and publicly available repositories. Bisulfite sequencing data generated for this study is available at the European Nucleotide Archive ( under accession number PRJEB71458. RNA-seq data from spleen of the six F. albicollis individuals are also available under the same accession number (PRJEB71458). Previously generated RNA-seq data is available at NCBI BioProject database ( under accession number PRJNA551584. Scripts are available at the GitHub repository:



Between-groups principal component analysis


Bateson–Dobzhansky–Muller incompatibility


CpG island


Conserved non-exonic element


Collared flycatcher


CG dinucleotide


Chicken repeat 1 transposon


Differential expression


Differentially methylated region

F ST :

Genetic differentiation


Gene ontology


F1 hybrid between female PIE and male COL


Long terminal repeat transposon

M diff :

Absolute methylation differentiation


Pied flycatcher


DMR between species, and between species and hybrids


Transposable element


Tissue-specific DMR


Transcription start site


Transcription termination site


Untranslated region


  1. Reifová R, Ament-Velásquez SL, Bourgeois Y, Coughlan J, Kulmuni J, Lipinska AP, et al. Mechanisms of intrinsic postzygotic isolation: from traditional genic and chromosomal views to genomic and epigenetic perspectives. Cold Spring Harb Perspect Biol. 2023;15(10):a041607.

    Article  PubMed  Google Scholar 

  2. Bateson W. Heredity and variation in modern lights. In: Seward AC, editor. Darwin and modern science. Cambridge: Cambridge University Press; 1909. p. 85–101.

    Google Scholar 

  3. Dobzhansky T. Studies on hybrid sterility – I. Spermatogenesis in pure and hybrid Drosophila pseudoobscura. Zeitschrift für Zellforsch und Mikroskopische Anat. 1934;21:169–223.

    Article  Google Scholar 

  4. Muller HJ. Bearing of the Drosophila work on systematics. In: Huxley JS, editor. The new systematics. Oxford: Clarendon Press; 1940. p. 185–268.

    Google Scholar 

  5. Coughlan JM, Matute DR. The importance of intrinsic postzygotic barriers throughout the speciation process. Philos Trans R Soc B. 1806;2020(375):20190533.

    Google Scholar 

  6. Presgraves DC. The molecular evolutionary basis of species formation. Nat Rev Genet. 2010;11:175–80.

    Article  CAS  PubMed  Google Scholar 

  7. Maheshwari S, Barbash DA. The genetics of hybrid incompatibilities. Annu Rev Genet. 2011;45:331–55.

    Article  CAS  PubMed  Google Scholar 

  8. Johnson NA, Porter AH. Rapid speciation via parallel, directional selection on regulatory genetic pathways. J Theor Biol. 2000;205:527–42.

    Article  CAS  PubMed  Google Scholar 

  9. Blevins T, Wang J, Pflieger D, Pontvianne F, Pikaard CS. Hybrid incompatibility caused by an epiallele. Proc Natl Acad Sci U S A. 2017;114:3702–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Jones FC, Grabherr MG, Chan YF, Russell P, Mauceli E, Johnson J, et al. The genomic basis of adaptive evolution in threespine sticklebacks. Nature. 2012;484:55–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Mack KL, Nachman MW. Gene regulation and speciation. Trends Genet. 2017;33:68–80.

    Article  CAS  PubMed  Google Scholar 

  12. Vrana PB, Fossella JA, Matteson P, Del Rio T, O’Neill MJ, Tilghman SM. Genetic and epigenetic incompatibilities underlie hybrid dysgenesis in peromyscus. Nat Genet. 2000;25:120–4.

    Article  CAS  PubMed  Google Scholar 

  13. Landry CR, Hartl DL, Ranz JM. Genome clashes in hybrids: insights from gene expression. Heredity. 2007;99:483–93.

    Article  CAS  PubMed  Google Scholar 

  14. Ortíz-Barrientos D, Counterman BA, Noor MAF. Gene expression divergence and the origin of hybrid dysfunctions. Genetica. 2007;129:71–81.

    Article  PubMed  Google Scholar 

  15. Bird A. DNA methylation patterns and epigenetic memory. Genes Dev. 2002;16:6–21.

    Article  CAS  PubMed  Google Scholar 

  16. Siegfried Z, Cedar H. DNA methylation: a molecular lock. Curr Biol. 1997;7:R305–7.

    Article  CAS  PubMed  Google Scholar 

  17. Walsh CP, Chaillet JR, Bestor TH. Transcription of IAP endogenous retroviruses is constrained by cytosine methylation. Nat Genet. 1998;20(2):116–7.

    Article  CAS  PubMed  Google Scholar 

  18. Suzuki MM, Bird A. DNA methylation landscapes: provocative insights from epigenomics. Nat Rev Genet. 2008;9:465–76.

    Article  CAS  PubMed  Google Scholar 

  19. Deniz Ö, Frost JM, Branco MR. Regulation of transposable elements by DNA modifications. Nat Rev Genet. 2019;20:417–31.

    Article  CAS  PubMed  Google Scholar 

  20. Yin Y, Morgunova E, Jolma A, Kaasinen E, Sahu B, Khund-Sayeed S, et al. Impact of cytosine methylation on DNA binding specificities of human transcription factors. Science (80- ). 2017;356(6337):eaaj2239.

    Article  Google Scholar 

  21. Tillotson R, Bird A. The molecular basis of MeCP2 function in the brain. J Mol Biol. 2020;432:1602–23.

    Article  CAS  PubMed  Google Scholar 

  22. Dennis K, Fan T, Geiman T, Yan Q, Muegge K. Lsh, a member of the SNF2 family, is required for genome-wide methylation. Genes Dev. 2001;15:2940–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13:484–92.

    Article  CAS  PubMed  Google Scholar 

  24. Jjingo D, Conley AB, Yi SV, Lunyak VV, King JI. On the presence and role of human gene-body DNA methylation. Oncotarget. 2012;3:462–74.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Zilberman D. An evolutionary case for functional gene body methylation in plants and animals. Genome Biol. 2017;18:1–3.

    Article  Google Scholar 

  26. Schmitz RJ, Lewis ZA, Goll MG. DNA methylation: shared and divergent features across eukaryotes. Trends Genet. 2019;35:818–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Boman J, Zhu Y, Höök L, Vila R, Talavera G, Backström N. Environmental stress during larval development induces DNA methylation shifts in the migratory painted lady butterfly (Vanessa cardui). Mol Ecol. 2023;32:3513–23.

    Article  CAS  PubMed  Google Scholar 

  28. Monk M, Boubelik M, Lehnert S. Temporal and regional changes in DNA methylation in the embryonic, extraembryonic and germ cell lineages during mouse embryo development. Development. 1987;99:371–82.

    Article  CAS  PubMed  Google Scholar 

  29. Bird A, Taggart M, Frommer M, Miller OJ, Macleod D. A fraction of the mouse genome that is derived from islands of nonmethylated. CpG-rich DNA Cell. 1985;40:91–9.

    CAS  PubMed  Google Scholar 

  30. Thomson JP, Skene PJ, Selfridge J, Clouaire T, Guy J, Webb S, et al. CpG islands influence chromatin structure via the CpG-binding protein Cfp1. Nature. 2010;464:1082–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Saxonov S, Berg P, Brutlag DL. A genome-wide analysis of CpG dinucleotides in the human genome distinguishes two distinct classes of promoters. Proc Natl Acad Sci. 2006;103:1412–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Weber M, Hellmann I, Stadler MB, Ramos L, Pääbo S, Rebhan M, et al. Distribution, silencing potential and evolutionary impact of promoter DNA methylation in the human genome. Nat Genet. 2007;39:457–66.

    Article  CAS  PubMed  Google Scholar 

  33. Nagae G, Isagawa T, Shiraki N, Fujita T, Yamamoto S, Tsutsumi S, et al. Tissue-specific demethylation in CpG-poor promoters during cellular differentiation. Hum Mol Genet. 2011;20:2710–21.

    Article  CAS  PubMed  Google Scholar 

  34. Villicaña S, Bell JT. Genetic impacts on DNA methylation: research findings and future perspectives. Genome Biol. 2021;22:127.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Cowles CR, Hirschhorn JN, Altshuler D, Lander ES. Detection of regulatory variation in mouse genes. Nat Genet. 2002;32:432–7.

    Article  CAS  PubMed  Google Scholar 

  36. Wittkopp PJ, Haerum BK, Clark AG. Evolutionary changes in cis and trans gene regulation. Nature. 2004;430:85–8.

    Article  CAS  PubMed  Google Scholar 

  37. Floc’hlay S, Wong ES, Zhao B, Viales RR, Thomas-Chollier M, Thieffry D, et al. Cis-acting variation is common across regulatory layers but is often buffered during embryonic development. Genome Res. 2021;31:211–24.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Devens HR, Davidson PL, Byrne M, Wray GA. Hybrid epigenomes reveal extensive local genetic changes to chromatin accessibility contribute to divergence in embryonic gene expression between species. Mol Biol Evol. 2023;40:msad222.

  39. Nadachowska-Brzyska K, Burri R, Smeds L, Ellegren H. PSMC analysis of effective population sizes in molecular ecology and its application to black-and-white Ficedula flycatchers. Mol Ecol. 2016;25:1058–72.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Qvarnström A, Rice AM, Ellegren H. Speciation in Ficedula flycatchers. Philos Trans R Soc B Biol Sci. 2010;365:1841–52.

    Article  Google Scholar 

  41. Svedin N, Wiley C, Veen T, Gustafsson L, Qvarnström A. Natural and sexual selection against hybrid flycatchers. Proc R Soc B Biol Sci. 2008;275:735–44.

    Article  Google Scholar 

  42. Alund M, Immler S, Rice AM, Qvarnstrom A. Low fertility of wild hybrid male flycatchers despite recent divergence. Biol Lett. 2013;9:20130169.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Ålund M, Segami Marzal JC, Zhu Y, Krishna Menon PN, Jones W, Qvarnström A. Tracking hybrid viability across life stages in a natural avian contact zone. Evolution (N Y). 2023.

    Article  Google Scholar 

  44. McFarlane SE, Sirkiä PM, Ålund M, Qvarnström A. Hybrid dysfunction expressed as elevated metabolic rate in male ficedula flycatchers. PLoS One. 2016;11:e0161547.

    Article  PubMed  PubMed Central  Google Scholar 

  45. van der Heijden E, McFarlane SE, Valk T van der, Qvarnström A. Divergent mitochondrial and nuclear OXPHOS genes are candidates for genetic incompatibilities in Ficedula Flycatchers. bioRxiv. 2019:588756.

  46. Li Q, Li N, Hu X, Li J, Du Z, Chen L, et al. Genome-wide mapping of DNA methylation in chicken. PLoS One. 2011;6:e19428.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Mugal CF, Arndt PF, Holm L, Ellegren H. Evolutionary consequences of DNA methylation on the GC content in vertebrate genomes. G3. 2015;5:441–7.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Laine VN, Gossmann TI, Schachtschneider KM, Garroway CJ, Madsen O, Verhoeven KJF, et al. Evolutionary signals of selection on cognition from the great tit genome and methylome. Nat Commun. 2016;7:1–9.

    Article  Google Scholar 

  49. Derks MFL, Schachtschneider KM, Madsen O, Schijlen E, Verhoeven KJF, van Oers K. Gene and transposable element methylation in great tit (Parus major) brain and blood. BMC Genomics. 2016;17:332.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Viitaniemi HM, Verhagen I, Visser ME, Honkela A, Van Oers K, Husby A, et al. Seasonal variation in genome-wide DNA methylation patterns and the onset of seasonal timing of reproduction in great tits. Genome Biol Evol. 2019;11:970–83.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Höglund A, Henriksen R, Fogelholm J, Churcher AM, Guerrero-Bosagna CM, Martinez-Barrio A, et al. The methylation landscape and its role in domestication and gene regulation in the chicken. Nat Ecol Evol. 2020;4:1713–24.

    Article  PubMed  Google Scholar 

  52. Lindner M, Verhagen I, Viitaniemi HM, Laine VN, Visser ME, Husby A, et al. Temporal changes in DNA methylation and RNA expression in a small song bird: within- and between-tissue comparisons. BMC Genomics. 2021;22:36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Sun D, Layman TS, Jeong H, Chatterjee P, Grogan K, Merritt JR, et al. Genome-wide variation in DNA methylation linked to developmental stage and chromosomal suppression of recombination in white-throated sparrows. Mol Ecol. 2021;30:3453–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Al Adhami H, Bardet AF, Dumas M, Cleroux E, Guibert S, Fauque P, et al. A comparative methylome analysis reveals conservation and divergence of DNA methylation patterns and functions in vertebrates. BMC Biol. 2022;20:1–18.

    Article  Google Scholar 

  55. Boman J, Qvarnström A, Mugal CF. Regulatory and evolutionary impact of DNA methylation. 2024.

    Google Scholar 

  56. Mugal CF, Wang M, Backström N, Wheatcroft D, Ålund M, Sémon M, et al. Tissue-specific patterns of regulatory changes underlying gene expression differences among Ficedula flycatchers and their naturally occurring F1 hybrids. Genome Res. 2020;31:1727–39.

    Article  Google Scholar 

  57. Mugal CF, Wang M, Backström N, Wheatcroft D, Ålund M, Sémon M, et al. Gene expression in natural F1 hybrids of Ficedula flycatchers and their parental species. 2020.

    Google Scholar 

  58. Hansen KD, Langmead B, Irizarry RA. BSmooth: from whole genome bisulfite sequencing reads to differentially methylated regions. Genome Biol. 2012;13(10):R83.

    Article  PubMed  PubMed Central  Google Scholar 

  59. Dickel DE, Ypsilanti AR, Pla R, Zhu Y, Barozzi I, Mannion BJ, et al. Ultraconserved enhancers are required for normal development. Cell. 2018;172:491–499.e15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Pan Z, Yao Y, Yin H, Cai Z, Wang Y, Bai L, et al. Pig genome functional annotation enhances the biological interpretation of complex traits and human disease. Nat Commun. 2021;12:1–15.

    Article  Google Scholar 

  61. Borgel J, Tyl M, Schiller K, Pusztai Z, Dooley CM, Deng W, et al. KDM2A integrates DNA and histone modification signals through a CXXC/PHD module and direct interaction with HP1. Nucleic Acids Res. 2017;45:1114–29.

    CAS  PubMed  Google Scholar 

  62. Turberfield AH, Kondo T, Nakayama M, Koseki Y, King HW, Koseki H, et al. KDM2 proteins constrain transcription from CpG island gene promoters independently of their histone demethylase activity. Nucleic Acids Res. 2019;47:9005–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Thioulouse J, Renaud S, Dufour AB, Dray S. Overcoming the spurious groups problem in between-group PCA. Evol Biol. 2021;48:458–71.

    Article  Google Scholar 

  64. Jackson BC, Campos JL, Zeng K. The effects of purifying selection on patterns of genetic differentiation between Drosophila melanogaster populations. Heredity (Edinb). 2015;114:163–74.

    Article  CAS  PubMed  Google Scholar 

  65. Burri R, Nater A, Kawakami T, Mugal CF, Olason PI, Smeds L, et al. Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers. Genome Res. 2015;25:1656–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. McManus CJ, Coolon JD, Duff MO, Eipper-Mains J, Graveley BR, Wittkopp PJ. Regulatory divergence in Drosophila revealed by mRNA-seq. Genome Res. 2010;20:816–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Coolon JD, McManus CJ, Stevenson KR, Graveley BR, Wittkopp PJ. Tempo and mode of regulatory evolution in Drosophila. Genome Res. 2014;24:797–808.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Metzger BPH, Wittkopp PJ, Coolon JD. Evolutionary dynamics of regulatory changes underlying gene expression divergence among Saccharomyces species. Genome Biol Evol. 2017;9:843–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Ellegren H, Smeds L, Burri R, Olason PI, Backström N, Kawakami T, et al. The genomic landscape of species divergence in Ficedula flycatchers. Nature. 2012;491:756–60.

    Article  CAS  PubMed  Google Scholar 

  70. Blake LE, Roux J, Hernando-Herraez I, Banovich NE, Perez RG, Hsiao CJ, et al. A comparison of gene expression and DNA methylation patterns across tissues and species. Genome Res. 2020;30:250–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Skene PJ, Illingworth RS, Webb S, Kerr ARW, James KD, Turner DJ, et al. Neuronal MeCP2 is expressed at near histone-octamer levels and globally alters the chromatin state. Mol Cell. 2010;37:457–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Nan X, Ng HH, Johnson CA, Laherty CD, Turner BM, Eisenman RN, et al. Transcriptional repression by the methyl-CpG-binding protein MeCP2 involves a histone deacetylase complex. Nature. 1998;393:386–9.

    Article  CAS  PubMed  Google Scholar 

  73. Husby A. Wild epigenetics: insights from epigenetic studies on natural populations. Proc R Soc B. 1968;2022(289):20211633.

    Google Scholar 

  74. Lea AJ, Vilgalys TP, Durst PAP, Tung J. Maximizing ecological and evolutionary insight in bisulfite sequencing data sets. Nat Ecol Evol. 2017;1:1074–83.

    Article  PubMed  PubMed Central  Google Scholar 

  75. Angeloni A, Bogdanovic O. Enhancer DNA methylation: implications for gene regulation. Essays Biochem. 2019;63:707–15.

    Article  CAS  PubMed  Google Scholar 

  76. Hernando-Herraez I, Heyn H, Fernandez-Callejo M, Vidal E, Fernandez-Bellon H, Prado-Martinez J, et al. The interplay between DNA methylation and sequence divergence in recent human evolution. Nucleic Acids Res. 2015;43:8204–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Thompson PJ, Macfarlan TS, Lorincz MC. Long terminal repeats: from parasitic elements to building blocks of the transcriptional regulatory repertoire. Mol Cell. 2016;62:766–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Bingham PM, Kidwell MG, Rubin GM. The molecular basis of P-M hybrid dysgenesis: the role of the P element, a P-strain-specific transposon family. Cell. 1982;29:995–1004.

    Article  CAS  PubMed  Google Scholar 

  79. Segami JC, Mugal CF, Cunha C, Bergin C, Schmitz M, Semon M, et al. The genomic basis of hybrid male sterility in Ficedula flycatchers. bioRxiv. 2022:2022.09.19.508503.

  80. Irizarry RA, Ladd-Acosta C, Wen B, Wu Z, Montano C, Onyango P, et al. The human colon cancer methylome shows similar hypo- and hypermethylation at conserved tissue-specific CpG island shores. Nat Genet. 2009;41:178–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Haghani A, Li CZ, Robeck TR, Zhang J, Lu AT, Ablaeva J, et al. DNA methylation networks underlying mammalian traits. Science (80- ). 2023;381:eabq5693.

    Article  CAS  Google Scholar 

  82. Rassoulzadegan M, Magliano M, Cuzin F. Transvection effects involving DNA methylation during meiosis in the mouse. EMBO J. 2002;21:440–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Lim B, Heist T, Levine M, Fukaya T. Visualization of transvection in living Drosophila embryos. Mol Cell. 2018;70:287–296.e6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Landry CR, Wittkopp PJ, Taubes CH, Ranz JM, Clark AG, Hartl DL. Compensatory cis-trans evolution and the dysregulation of gene expression in interspecific hybrids of drosophila. Genetics. 2005;171:1813–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. McGirr JA, Martin CH. Ecological divergence in sympatry causes gene misexpression in hybrids. Mol Ecol. 2020;29:2707–21.

    Article  PubMed  PubMed Central  Google Scholar 

  86. Sánchez-Ramírez S, Weiss J, Thomas C, Cutter A. Widespread misregulation of inter-species hybrid transcriptomes due to sex-specific and sex-chromosome regulatory evolution. PLoS Genet. 2021;17:e1009409.

    Article  PubMed  PubMed Central  Google Scholar 

  87. Chase MA, Ellegren H, Mugal CF. Positive selection plays a major role in shaping signatures of differentiation across the genomic landscape of two independent Ficedula flycatcher species pairs*. Evolution (N Y). 2021;75:2179–96.

    Google Scholar 

  88. Lemos B, Meiklejohn CD, Cáceres M, Hartl DL. Rates of divergence in gene expression profiles of primates, mice, and flies: stabilizing selection and variability among functional categories. Evolution (N Y). 2005;59:126–37.

    CAS  Google Scholar 

  89. Chen J, Swofford R, Johnson J, Cummings BB, Rogel N, Lindblad-Toh K, et al. A quantitative framework for characterizing the evolutionary history of mammalian gene expression. Genome Res. 2019;29:53–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Barreto FS, Pereira RJ, Burton RS. Hybrid dysfunction and physiological compensation in gene expression. Mol Biol Evol. 2015;32:613–22.

    Article  CAS  PubMed  Google Scholar 

  91. Smith TA, Martin MD, Nguyen M, Mendelson TC. Epigenetic divergence as a potential first step in darter speciation. Mol Ecol. 2016;25:1883–94.

    Article  CAS  PubMed  Google Scholar 

  92. Greenspoon PB, Spencer HG, M’Gonigle LK. Epigenetic induction may speed up or slow down speciation with gene flow. Evolution (N Y). 2022;76:1170–82.

    CAS  Google Scholar 

  93. Planidin NP, de Carvalho CF, Feder JL, Gompert Z, Nosil P. Epigenetics and reproductive isolation: a commentary on Westram et al., 2022. J Evol Biol. 2022;35:1188–94.

    Article  PubMed  PubMed Central  Google Scholar 

  94. Frésard L, Leroux S, Servin B, Gourichon D, Dehais P, Cristobal MS, et al. Transcriptome-wide investigation of genomic imprinting in chicken. Nucleic Acids Res. 2014;42:3768–82.

    Article  PubMed  PubMed Central  Google Scholar 

  95. Guerrero-Bosagna C, Morisson M, Liaubet L, Rodenburg TB, de Haas EN, Košťál Ľ, et al. Transgenerational epigenetic inheritance in birds. Environ Epigenetics. 2018;4(2):dvy008.

    Article  Google Scholar 

  96. Gore AV, Tomins KA, Iben J, Ma L, Castranova D, Davis AE, et al. An epigenetic mechanism for cavefish eye degeneration. Nat Ecol Evol. 2018;2:1155–60.

    Article  PubMed  PubMed Central  Google Scholar 

  97. Di Tommaso P, Chatzou M, Floden EW, Barja PP, Palumbo E, Notredame C. Nextflow enables reproducible computational workflows. Nat Biotechnol. 2017;35:316–9.

    Article  PubMed  Google Scholar 

  98. Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010.

    Google Scholar 

  99. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17:10.

    Article  Google Scholar 

  100. Kawakami T, Smeds L, Backström N, Husby A, Qvarnström A, Mugal CF, et al. A high-density linkage map enables a second-generation collared flycatcher genome assembly and reveals the patterns of avian recombination rate variation and chromosomal evolution. Mol Ecol. 2014;23:4035–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011;27:1571–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  102. Okonechnikov K, Conesa A, García-Alcalde F. Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics. 2016;32:292–4.

    Article  CAS  PubMed  Google Scholar 

  103. Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32:3047–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  104. MacManes MD. The Oyster River Protocol: a multi-assembler and kmer approach for de novo transcriptome assembly. PeerJ. 2018;2018:e5428.

    Article  Google Scholar 

  105. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  106. Song L, Florea L. Rcorrector: efficient and accurate error correction for Illumina RNA-seq reads. Gigascience. 2015;4:48.

    Article  PubMed  PubMed Central  Google Scholar 

  107. Nurk S, Bankevich A, Antipov D, Gurevich A, Korobeynikov A, Lapidus A, et al. Assembling genomes and mini-metagenomes from highly chimeric reads. Lect Notes Comput Sci (including Subser Lect Notes Artif Intell Lect Notes Bioinformatics). 2013;7821 LNBI:158–70.

    Google Scholar 

  108. Robertson G, Schein J, Chiu R, Corbett R, Field M, Jackman SD, et al. De novo assembly and analysis of RNA-seq data. Nat Methods. 2010;7(11):909–12.

    Article  CAS  PubMed  Google Scholar 

  109. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  110. Emms DM, Kelly S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 2019;20:1–14.

    Article  Google Scholar 

  111. Smith-Unna R, Boursnell C, Patro R, Hibberd JM, Kelly S. TransRate: reference-free quality assessment of de novo transcriptome assemblies. Genome Res. 2016;26:1134–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  112. Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2.

    Article  PubMed  Google Scholar 

  113. Holt C, Yandell M. MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. BMC Bioinformatics. 2011.

    Article  PubMed  PubMed Central  Google Scholar 

  114. Uebbing S, Künstner A, Mäkinen H, Backström N, Bolivar P, Burri R, et al. Divergence in gene expression within and between two closely related flycatcher species. Mol Ecol. 2016;25:2015–28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  115. Botero-Castro F, Figuet E, Tilak MK, Nabholz B, Galtier N. Avian genomes revisited: hidden genes uncovered and the rates versus traits paradox in birds. Mol Biol Evol. 2017;34:3123–31.

    Article  CAS  PubMed  Google Scholar 

  116. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.

    Article  CAS  PubMed  Google Scholar 

  117. Suh A, Smeds L, Ellegren H. Abundant recent activity of retrovirus-like retrotransposons within and among flycatcher species implies a rich source of structural variation in songbird genomes. Mol Ecol. 2018;27:99–111.

    Article  CAS  PubMed  Google Scholar 

  118. Weissensteiner MH, Bunikis I, Catalán A, Francoijs KJ, Knief U, Heim W, et al. Discovery and population genomics of structural variation in a songbird genus. Nat Commun. 2020;11:3403.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  119. Boman J, Frankl-Vilches C, Santos MDSD, de Oliveira EHC, Gahr M, Suh A. The genome of blue-capped cordon-bleu uncovers hidden diversity of LTR retrotransposons in zebra finch. Genes (Basel). 2019;10:301.

    Article  CAS  PubMed  Google Scholar 

  120. Peona V, Blom MPK, Xu L, Burri R, Sullivan S, Bunikis I, et al. Identifying the causes and consequences of assembly gaps using a multiplatform genome assembly of a bird-of-paradise. Mol Ecol Resour. 2021;21:263–86.

    Article  PubMed  Google Scholar 

  121. Prost S, Armstrong EE, Nylander J, Thomas GWC, Suh A, Petersen B, et al. Comparative analyses identify genomic features potentially involved in the evolution of birds-of-paradise. Gigascience. 2019;8:1–12.

    Article  CAS  Google Scholar 

  122. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    Article  CAS  PubMed  Google Scholar 

  123. Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14(4):417–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  124. Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics. 2012;28:2184–5.

    Article  CAS  PubMed  Google Scholar 

  125. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

    Article  PubMed  PubMed Central  Google Scholar 

  126. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  127. Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res. 2016;4:1521.

    Article  PubMed Central  Google Scholar 

  128. Hackenberg M, Previti C, Luque-Escamilla PL, Carpena P, Martínez-Aroza J, Oliver JL. CpGcluster: a distance-based algorithm for CpG-island detection. BMC Bioinformatics. 2006;7:1–13.

    Article  Google Scholar 

  129. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  130. Siepel A, Bejerano G, Pedersen JS, Hinrichs AS, Hou M, Rosenbloom K, et al. Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. Genome Res. 2005;15:1034–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  131. Craig RJ, Suh A, Wang M, Ellegren H. Natural selection beyond genes: identification and analyses of evolutionarily conserved elements in the genome of the collared flycatcher ( Ficedula albicollis ). Mol Ecol. 2018;27:476–92.

    Article  CAS  PubMed  Google Scholar 

  132. R Core Team. R: a language and environment for statistical computing. 2020.

    Google Scholar 

  133. Ewens WJ. On estimating p values by Monte Carlo methods. Am J Hum Genet. 2003;72:496–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  134. Ge SX, Jung D, Jung D, Yao R. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics. 2020;36:2628–9.

    Article  CAS  PubMed  Google Scholar 

  135. Huminiecki L, Lloyd AT, Wolfe KH. Congruence of tissue expression profiles from gene expression Atlas, SAGEmap and TissueInfo databases. BMC Genomics. 2003;4:1–10.

    Article  Google Scholar 

  136. Fraser HB. Improving estimates of compensatory cis–trans regulatory divergence. Trends Genet. 2019;35:3–5.

    Article  CAS  PubMed  Google Scholar 

  137. Nater A, Burri R, Kawakami T, Smeds L, Ellegren H. Resolving evolutionary relationships in closely related species with whole-genome sequencing data. Syst Biol. 2015;64:1000–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  138. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution (N Y). 1984;38:1358.

    CAS  Google Scholar 

  139. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  140. Krueger F. SNPsplit v0.3.2. 2017.

    Google Scholar 

  141. Hebestreit K, Dugas M, Klein H-U. Detection of significantly differentially methylated regions in targeted bisulfite sequencing data. Bioinformatics. 2013;29:1647–53.

    Article  CAS  PubMed  Google Scholar 

Download references


The authors thank William Jones, Murielle Ålund, S. Eryn McFarlane, David Wheatcroft, Niclas Backström, and Luohao Xu for field and lab work related to this study. The authors are also grateful to Rory J. Craig for contributing to the updated gene annotation and Hans Ellegren, who together with A.Q. and C.F.M. conceived of the study and contributed to the study with the whole-genome bisulfite sequencing data generation.


Open access funding provided by Uppsala University. This study was funded by grants from the Swedish Research Council (2013- 8271 to Hans Ellegren and 2012-3722 to A.Q.) and the Knut and Alice Wallenberg Foundation (2014/0044 to Hans Ellegren). Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala. The facility is part of the National Genomics Infrastructure (NGI) Sweden and Science for Life Laboratory. The SNP&SEQ Platform is also supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation. Computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725.

Author information

Authors and Affiliations



JB: methodology, software, formal analysis, writing—original draft, and visualization. AQ: conceptualization, resources, writing—review and editing, and supervision. CFM: conceptualization, writing—original draft, and supervision.

Corresponding authors

Correspondence to Jesper Boman or Carina F. Mugal.

Ethics declarations

Ethics approval and consent to participate

All sampling procedures were approved by the Swedish Board of Agriculture (Jordbruksverket—DNR 21-11).

Consent for publication

Not applicable.

Competing interests

We declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1. Sample and sequencing information.


Additional file 2: Supplementary Results. Fig. S1. Figure 1A-L in main manuscript colored by sample group (species and hybrids). Fig. S2. Methylation gene profile with gene body split into exons and introns. Fig. S3. Tissue-specific patterns of the association between genetic- and methylation differentiation. Fig. S4. CGI proportion gene profile. Fig. S5. Patterns of genetic and epigenetic change at misexpressed genes. Table S2. Frequency of hypermethylation of tissue-specific DMRs. Table S7. Classification system for the mechanism of DNA methylation divergence. Table S8. Number of fixed difference loci by divergence class. Table S9. Number of differentially expressed genes. Table S10. Determinants in cis of misexpressed genes.

Additional file 3: Table S3. tsDMR GO analysis for brain.

Additional file 4: Table S4. tsDMR GO analysis for kidney.

Additional file 5: Table S5. tsDMR GO analysis for liver.

Additional file 6: Table S6. tsDMR GO analysis for testis.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Boman, J., Qvarnström, A. & Mugal, C.F. Regulatory and evolutionary impact of DNA methylation in two songbird species and their naturally occurring F1 hybrids. BMC Biol 22, 124 (2024).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: