<figure><img src="image_1.png" /></figure>

HAL
open science

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

Jesper Boman, Anna Qvarnström, Carina Mugal

► To cite this version:

Jesper Boman, Anna Qvarnström, Carina Mugal. Regulatory and evolutionary impact of DNA methylation in two songbird species and their naturally occurring F1 hybrids. BMC Biology, 2024, 22 (1), pp.124. 10.1186/s12915-024-01920-2. hal-04604325

HAL Id: hal-04604325

https://hal.science/hal-04604325v1

Submitted on 17 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# Regulatory and evolutionary impact of DNA methylation in two songbird species and their naturally occurring F₁ hybrids

Jesper Boman¹*, D. Anna Qvarnström² and Carina F. Mugal¹,³*

**Abstract**

**Background** 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 F₁ hybrids.

**Results** 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.

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

**Keywords** DNA methylation, Speciation, Transcriptomics, Epigenomics

**Background**

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–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

<figure><img src="image_4.png" /></figure>

BMC

© The Author(s) 2024. **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 http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

<footnote>*Correspondence:
Jesper Boman
jesper.boman@ebc.uu.se
Carina F. Mugal
carina.mugal@univ-lyon1.fr
¹Department of Ecology and Genetics (IEG), Division of Evolutionary Biology, Uppsala University, Norbyvägen 18D, Uppsala SE-752 36, Sweden
²Department of Ecology and Genetics (IEG), Division of Animal Ecology, Uppsala University, Norbyvägen 18D, Uppsala SE-752 36, Sweden
³CNRS, Laboratory of Biometry and Evolutionary Biology (LBBE), UMR 5558, University of Lyon 1, Villeurbanne, France</footnote>
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–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–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–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 F₁ 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 F₁ hybrids. The pied flycatcher (*Ficedula hypoleuca*) and the collared flycatcher (*F. albicollis*) are two species of songbirds that diverged approximately 1 MYA [39]. Inter-specific mating occurs in sympatry, e.g., on the island of Öland in the Baltic Sea [40]. F₁ 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–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 F₁ 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 differ-
ential gene expression among species as well as in hybrid
misexpression.

Results

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 occur-
ring F₁ hybrids, sampled on the island of Öland in the
Baltic Sea (Fig. 1). In total, we sequenced 6 collared fly-
catchers, 5 pied flycatchers, and 3 F₁ hybrids (41.5 billion
reads in total) [55]. All F₁ 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 (Addi-
tional 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 sam-
ples [56, 57], thus forming matched genomic data sets.

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 annota-
tion (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 pro-
moters, referred to as CGI and *Other* promoters, respec-
tively. 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 expres-
sion 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 sepa-
rately 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 tes-
tis 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 meth-
ylation level in exons (Additional file 2: Fig. S2). Genes
with CGI promoters showed a drop in DNA methyla-
tion 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 cat-
egories, 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.

To further assess the relationship between DNA
methylation and gene expression, we computed a gene

<figure><img src="image_5.png" /><figcaption>**Fig. 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</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 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 (A–E) and Other (F–J) 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</figcaption></figure>

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.

It is possible that functionally important tissue-specific patterns are shared between PIE and COL and potentially also present in viable F₁ 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⁻¹⁶; 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.

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
## Tissue-specific DMRs

A Brain

<figure><img src="image_7.png" /></figure>

B Heart

<figure><img src="image_8.png" /></figure>

C Kidney

<figure><img src="image_9.png" /></figure>

D Liver

<figure><img src="image_10.png" /></figure>

E Testis

<figure><img src="image_11.png" /></figure>

F

<figure><img src="image_12.png" /></figure>

## Between-species DMRs

G

<figure><img src="image_13.png" /></figure>

**Fig. 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; **A–E**). 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 **F**, **G**, "fixed differences" refer to the 100 ± bp vicinity of fixed differences
<table><caption>**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)</caption><thead><tr><th></th><th>Brain</th><th>Heart</th><th>Kidney</th><th>Liver</th><th>Testis</th></tr></thead><tbody><tr><td><strong>Number of spDMRs</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>COL vs. PIE</td><td>11,330</td><td>9182</td><td>10,087</td><td>10,479</td><td>7430</td></tr><tr><td>COL vs. HYB</td><td>8474</td><td>6863</td><td>6051</td><td>4869</td><td>5559</td></tr><tr><td>PIE vs. HYB</td><td>9064</td><td>7545</td><td>7314</td><td>8766</td><td>6190</td></tr><tr><td><strong>Frequency of hypermethylation</strong> (higher in sample group 1, e.g., COL)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>COL vs. PIE</td><td>0.68 ***</td><td>0.67 ***</td><td>0.70 ***</td><td>0.67 ***</td><td>0.63 ***</td></tr><tr><td>COL vs. HYB</td><td>0.62 ***</td><td>0.29 ***</td><td>0.50</td><td>0.57 ***</td><td>0.83 ***</td></tr><tr><td>PIE vs. HYB</td><td>0.42 ***</td><td>0.16 ***</td><td>0.28 ***</td><td>0.28 ***</td><td>0.66 ***</td></tr></tbody></table>

Family-wise (0.1) adjusted p-value levels: * < 0.05 ** < 0.005 *** < 0.0005
$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><caption>**Table 2** Between-groups principal component analysis. $R^2$ is the proportion of variance explained. Values were rounded to two decimal points</caption><thead><tr><th>Annotation</th><th><i>R</i><sup>2</sup><sub>between tissues</sub></th><th><i>R</i><sup>2</sup><sub>between species</sub></th><th><i>R</i><sup>2</sup><sub>between species controlled for tissue</sub></th></tr></thead><tbody><tr><td>Promoters (CGI)</td><td>0.75 *</td><td>0.03</td><td>0.12 *</td></tr><tr><td>Promoters (Other)</td><td>0.65 **</td><td>0.05</td><td>0.14 *</td></tr><tr><td>Genes</td><td>0.88 *</td><td>0.02</td><td>0.13 *</td></tr><tr><td>Introns</td><td>0.55 *</td><td>0.09 *</td><td>0.20 *</td></tr><tr><td>CNEEs</td><td>0.61 *</td><td>0.04</td><td>0.11 *</td></tr><tr><td>CDSs</td><td>0.67 *</td><td>0.04</td><td>0.11 *</td></tr><tr><td>5' UTRs</td><td>0.68 *</td><td>0.04</td><td>0.13 **</td></tr><tr><td>3' UTRs</td><td>0.64 *</td><td>0.05</td><td>0.15 *</td></tr><tr><td>TEs</td><td>0.51 *</td><td>0.11 *</td><td>0.23 *</td></tr><tr><td>CR1s</td><td>0.49 *</td><td>0.12 *</td><td>0.24 *</td></tr><tr><td>LTRs</td><td>0.55 *</td><td>0.10 *</td><td>0.22 *</td></tr><tr><td>Fixed differences</td><td>0.47 *</td><td>0.16 *</td><td>0.30 *</td></tr></tbody></table>

Family-wise (0.1) adjusted $p$-value levels: * < 0.05 ** < 0.005 *** < 0.0005

## 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.

## 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
<figure><img src="image_14.png" /><figcaption>**Fig. 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</figcaption></figure>

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).
<figure><img src="image_15.png" /><figcaption>**Fig. 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</figcaption></figure>

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 allelic-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.

### 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.

### Discussion

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].

<table><caption>**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</caption><thead><tr><th>Tissue</th><th>Promoter type</th><th>P<sub>DMR freq</sub></th><th>GB<sub>DMR freq</sub></th><th>P<sub>non-CpG Fst</sub></th><th>GB<sub>non-CpG Fst</sub></th></tr></thead><tbody><tr><td>Brain</td><td>CGI</td><td>1</td><td>0.456 (-)</td><td>0.041 *</td><td>0 ***</td></tr><tr><td>Heart</td><td>CGI</td><td>0.023 *</td><td>0.039 *</td><td>0.1</td><td>0 ***</td></tr><tr><td>Kidney</td><td>CGI</td><td>0.006 *</td><td>1</td><td>1</td><td>0 ***</td></tr><tr><td>Liver</td><td>CGI</td><td>0.003 **</td><td>0 ***</td><td>0.125</td><td>0 ***</td></tr><tr><td>Testis</td><td>CGI</td><td>0.367</td><td>1 (-)</td><td>0.001 **</td><td>0 ***</td></tr><tr><td>Brain</td><td>Other</td><td>1 (-)</td><td>0.069</td><td>0.013 *</td><td>0 ***</td></tr><tr><td>Heart</td><td>Other</td><td>0.660</td><td>0.2</td><td>1</td><td>0 ***</td></tr><tr><td>Kidney</td><td>Other</td><td>0.171 (-)</td><td>0.023 *</td><td>1 (-)</td><td>0 ***</td></tr><tr><td>Liver</td><td>Other</td><td>0.020 (-)</td><td>1</td><td>1</td><td>0 ***</td></tr><tr><td>Testis</td><td>Other</td><td>1</td><td>1</td><td>0.456</td><td>0 ***</td></tr></tbody></table>

<footnote>(–) average value lower in DE genes, *P* Promoter, *GB* Gene body</footnote>

<footnote>Family-wise (0.1) adjusted *p*-value levels: * < 0.05 ** < 0.005 *** < 0.0005</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 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</figcaption></figure>

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 spD-MRs 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 $F_1$ 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 $F_1$ 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 $F_1$ hybrids of a wide range of species [14, 56, 84–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 F₁ 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 F₁ hybrid flycatchers [44, 45] and somewhat in line with results in a copepod with known mitochondrial dysfunction in F₂ 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–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.

## Conclusions

In this study, we investigated the relationship between genetic differentiation, DNA methylation, and gene expression in two songbird species and their F₁ 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.

## Methods

### 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 F₁ 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 Agri-culture (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 fly-catcher 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 F₁ 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),

$$ \text{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 TransABySS 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 coridon-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 CpG\_CLUSTER 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

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:

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

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

### 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 χ² 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 over-expressed (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 F₁ 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 F₁ 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 F₁ hybrids H. If we assume that the phenotypic value P is pairwise-independent between A, B, and H then,

$$cov(P_H - P_A, P_H - P_B) = var(P_H).$$

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_H - P_A and another for P_H in P_H - P_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 (E. 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<sub>ST</sub> [138] implemented in vcftools v0.1.16 [139].

## Allele-specific methylation estimation

We estimated allele-specific methylation in F₁ 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 F₁ 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 F₁ 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 $F_1$ 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.

**Abbreviations**

<table><tr><td>BCA</td><td>Between-groups principal component analysis</td></tr><tr><td>BDMI</td><td>Bateson-Dobzhansky-Muller incompatibility</td></tr><tr><td>CGI</td><td>CpG island</td></tr><tr><td>CNEE</td><td>Conserved non-exonic element</td></tr><tr><td>COL</td><td>Collared flycatcher</td></tr><tr><td>CpG</td><td>CG dinucleotide</td></tr><tr><td>CR1</td><td>Chicken repeat 1 transposon</td></tr><tr><td>DE</td><td>Differential expression</td></tr><tr><td>DMR</td><td>Differentially methylated region</td></tr><tr><td>F<sub>ST</sub></td><td>Genetic differentiation</td></tr><tr><td>GO</td><td>Gene ontology</td></tr><tr><td>HYB</td><td>F<sub>1</sub> hybrid between female PIE and male COL</td></tr><tr><td>LTR</td><td>Long terminal repeat transposon</td></tr><tr><td>M<sub>diff</sub></td><td>Absolute methylation differentiation</td></tr><tr><td>PIE</td><td>Pied flycatcher</td></tr><tr><td>spDMR</td><td>DMR between species, and between species and hybrids</td></tr><tr><td>TE</td><td>Transposable element</td></tr><tr><td>tsDMR</td><td>Tissue-specific DMR</td></tr><tr><td>TSS</td><td>Transcription start site</td></tr><tr><td>TTS</td><td>Transcription termination site</td></tr><tr><td>UTR</td><td>Untranslated region</td></tr></table>

**Supplementary Information**

The online version contains supplementary material available at https://doi.org/10.1186/s12915-024-01920-2.

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.

**Acknowledgements**

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.

**Authors' contributions**

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.

**Funding**

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.

**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 (https://www.ebi.ac.uk/ena/browser/) 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 (https://www.ncbi.nlm.nih.gov/bioproject/) under accession number PRJNA551584. Scripts are available at the GitHub repository: https://github.com/Jespe/rBoman/DNA_methylation_flycatchers.

**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.

Received: 19 January 2024 Accepted: 15 May 2024

Published online: 29 May 2024

**References**

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.

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.

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.

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

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

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

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

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

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.

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.

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

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.

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

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

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

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

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

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

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

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.

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

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

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

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.

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

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

27. Boman J, Zhu Y, Höök G, 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.

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.

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.

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.

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.

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.

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.

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

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

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

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

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.

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

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

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

43. Ålund M, Segami Marzal JC, Zhu Y, Krishna Menon PN, Jones W, Qvarström A. Tracking hybrid viability across life stages in a natural avian contact zone. Evolution (N Y). 2023. https://doi.org/10.1093/evolut/ QPAD204

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

45. van der Heijden E, McFarlane SE, Valk T van der, Qvarströ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.

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.

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.

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.
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.

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.

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.

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.

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.

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

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.

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.

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

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.

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.

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.

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.

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

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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.

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

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.

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.

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.

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.

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.

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

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

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.

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

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.

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.

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.

89. Chen J, Sofford 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.

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

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.

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.

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.

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.
95. Guerrero-Bosagna C, Morisson M, Liaubet L, Rodenburg TB, de Haas EN, Koštál L, et al. Transgenerational epigenetic inheritance in birds. Environ Epigenetics. 2018;4(2):dvy008.

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.

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

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

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

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.

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

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.

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.

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

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

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

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.

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.

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.

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

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.

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.

113. Holt C, Yandell M. MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. BMC Bioinformatics. 2011. https://doi.org/10.1186/1471-2105-12-491.

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.

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.

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

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.

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.

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.

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 multi-platform genome assembly of a bird-of-paradise. Mol Ecol Resour. 2021;21:263–86.

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.

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.

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.

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

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.

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.

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

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

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

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.

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.

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

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

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.

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.

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

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.

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

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.

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

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

**Publisher's Note**

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