Seasonal Variation in Genome-Wide DNA Methylation Patterns and the Onset of Seasonal Timing of Reproduction in Great Tits.

In seasonal environments, timing of reproduction is a trait with important fitness consequences, but we know little about the molecular mechanisms that underlie the variation in this trait. Recently, several studies put forward DNA methylation as a mechanism regulating seasonal timing of reproduction in both plants and animals. To understand the involvement of DNA methylation in seasonal timing of reproduction, it is necessary to examine within-individual temporal changes in DNA methylation, but such studies are very rare. Here, we use a temporal sampling approach to examine changes in DNA methylation throughout the breeding season in female great tits (Parus major) that were artificially selected for early timing of breeding. These females were housed in climate-controlled aviaries and subjected to two contrasting temperature treatments. Reduced representation bisulfite sequencing on red blood cell derived DNA showed genome-wide temporal changes in more than 40,000 out of the 522,643 CpG sites examined. Although most of these changes were relatively small (mean within-individual change of 6%), the sites that showed a temporal and treatment-specific response in DNA methylation are candidate sites of interest for future studies trying to understand the link between DNA methylation patterns and timing of reproduction.

Introduction

In seasonal environments, timing of reproduction has major effects on fitness (Thomas et al. 2001; Zera and Harshman 2001; Lane et al. 2012; Reed et al. 2013; Lu et al. 2016) and is therefore under strong selection, as individuals need to time their reproduction to match the time window when food resources are most abundant (Visser et al. 2006; Plard et al. 2014). As this time window varies between years, plants and animals use environmental cues, such as photoperiod and temperature, to time their reproduction (Shindo et al. 2006; Perfito et al. 2012; Schaper et al. 2012; Stevenson and Prendergast 2013). The within-individual change in phenotype in response to environmental factors, that is, phenotypic plasticity (Pigliucci 2005), allows individuals to adjust to changing environmental conditions and plasticity itself can also be under genetic control (Renn et al. 2008; Espinosa-Soto et al. 2011). However, the molecular mechanisms via which cue perception leads to the onset of reproduction are poorly understood (Caro et al. 2013), although gene regulatory mechanisms are thought to play a major role. Factors that modulate regulatory changes, such as epigenetic mechanisms, are therefore of key interest in understanding the molecular basis of phenological traits (Beaman et al. 2016).

Epigenetic mechanisms can broadly be defined as modifications of the DNA that affect gene expression, without alterations in the DNA sequence itself (Allis and Jenuwein 2016). Of the different types of epigenetic modifications, DNA methylation is the most commonly studied form (Laird 2003) and involves the addition of a methyl group to the fifth carbon of a cytosine residue in the DNA backbone (Klose and Bird 2006). In vertebrates, methylation predominantly occurs when the cytosine is adjacent to guanidine in the DNA (CpG site) although non-CpG methylation is also observed (Shirane et al. 2013; Guo et al. 2014; Derks et al. 2016). The addition of a methyl group at the promoter and/or in the transcription start site (TSS) regions typically impedes binding affinity of transcription factors and reduce transcription rate and therefore gene expression (Jaenisch and Bird 2003). Continuous stretches of CpG dinucleotides in the DNA form so-called CpG islands, which are distributed across different genomic regions (Deaton and Bird 2011). Methylation of CpG sites and CpG islands in other genomic contexts has been shown to have varying effects on gene expression: for example, intragenic methylation is involved in regulation of alternative splicing (Lev Maor et al. 2015), prevention of spurious transcription (Neri et al. 2017), and expression of regulatory ncRNA’s (Deaton and Bird 2011). As mechanisms of active DNA demethylation have become better understood (Wu and Zhang 2010, 2017) a role for DNA methylation as rapid regulator in gene expression due to its environmentally responsive nature is evident (Stevenson and Prendergast 2013; Lynch et al. 2016; Weyrich et al. 2016).

Interestingly, DNA methylation patterns have been linked to timing of reproduction in different taxa in several recent studies. For example, removal of DNA methylation alters the pattern of flowering time in Arabidopsis and Taraxacum officinale (Cortijo et al. 2014; Wilschut et al. 2016). Furthermore, changes in photoperiod and temperature affect DNA methylation patterns of flower buds in Azalea (Meijón et al. 2011). Photoperiod also seems to affect DNA methylation patterns in insects: Nasonia vitripennis exposed to short or long day light cycles showed differences in methylation pattern and diapause response which were abolished by experimental lowering of DNA methylation (Pegoraro et al. 2016). However, although there is increasing evidence for the role of DNA methylation in seasonal timing of reproduction in plants and insects, similar studies in vertebrates are scarce. The only study so far, to our knowledge, to experimentally test the role of DNA methylation in regulating seasonal timing of reproduction investigated the effect of photoperiod on hypothalamic DNA methylation in Siberian hamsters (Stevenson and Prendergast 2013). This pioneering study showed that photoperiodic cues alter the DNA methylation pattern of the deiodinase 3 (DIO3) gene promoter, which is involved in thyroid hormone signaling, and play a key role in reproductive related behaviors and physiology (Dellovade et al. 1996; Barrett et al. 2007).

Although there is accumulating evidence to suggest that temporal changes in DNA methylation can be important in regulating seasonal timing of reproduction, most studies use photoperiod as an environmental cue for explaining methylation changes (Stevenson and Prendergast 2013; Pegoraro et al. 2016). We know much less about the effects of temperature and other environmental cues on fluctuations in DNA methylation patterns as a mechanism for fine tuning seasonal timing of reproduction. This is important to know, because the changes in phenology observed in many plants and animals over the last decades are not due to a shift in photoperiod, but rather due to changes in temperature, or environmental cues correlated to temperature (Reale et al. 2003; Nevo et al. 2012). For example, in the great tit (Parus major), lay dates have advanced with increasing spring temperatures over the last decades in many populations across Europe (Charmantier et al. 2008; Husby et al. 2010; Valtonen et al. 2017). In the long-term study population in Hoge Veluwe national park in the Netherlands mean laying dates have advanced by 0.19 days/year since 1973 (Husby et al. 2010).

The great tit is a particularly well-suited species to examine temporal patterns of DNA methylation in relation to seasonal timing of reproduction, because both causes (Visser et al. 2009) and consequences (Visser et al. 2006) of seasonal timing of reproduction are well understood (Husby et al. 2011; Schaper et al. 2012; te Marvelde et al. 2012) and it has very good genomic resources in the form of both a high-quality annotated genome and methylome (Derks et al. 2016; Laine et al. 2016). In this study, we were interested in temporal variation in DNA methylation to test if, as suggested by other studies, DNA methylation plays a role in determining individual variation in seasonal timing. We therefore 1) test for genome-wide temporal variation in DNA methylation and 2) whether temperature affects temporal changes in DNA methylation. We address this by using female great tits from a line artificially selected for early timing of reproduction that were housed in climate-controlled aviaries and were subjected to contrasting temperatures.

Materials and Methods

Experimental Design

This study included samples (see “Blood sample collection”) from female great tits that are part of a large-scale experiment where females are artificially selected for early or late reproduction using genomic selection. The experimental design for this artificial selection experiment is described in more detail in Verhagen I, Gienapp P, Laine VN, Mateman AC, van Oers K, Visser ME (unpublished data). In short, offspring (F1) from a wild parental generation (P) from the Hoge Veluwe (the Netherlands) with either early (early line) or late (late line) lay dates were collected and genotyped with a 650k single nucleotide polymorphism (SNP) chip to calculate genomic estimated breeding values (GEBVs), the value of an individual in the breeding scheme based on the estimated genomic marker (i.e., SNPs) effects throughout the genome (see Verhagen I, Gienapp P, Laine VN, Mateman AC, van Oers K, Visser ME unpublished data; Gienapp P, Calus MPL, Laine VN, Visser ME, unpublished data, for details). Based on their GEBVs, birds were selected within the early or late reproducing selection line in such a way that birds with the most extreme GEBVs were housed in pairs in outdoor aviaries to generate the F1-generation. Birds from the F1-generation, in their turn, were genotyped and selected for extreme GEBVs within the two selection lines to produce the F2-generation. From these F2 birds, pairs were put in climate-controlled aviaries and subjected to two contrasting temperature treatments from January until July. The two temperature treatments mimicked a cold and warm year in the Netherlands using observed temperatures from the year 2013 and 2014, respectively (see Verhagen I, Laine VN, Mateman AC, Pijl A, de Wit R, van Lith B, Kamphuis W, Viitaniemi HM, Williams TD, Caro SP, Meddle SL, Gienapp P, van Oers K & Visser ME (unpublished data)). Photoperiod changed daily following the natural photoperiod. Aviaries were provided with three nest boxes and from March onward, nesting material was provided.

Blood Sample Collection

A detailed blood sampling scheme is described in Mäkinen et al. (forthcoming) but briefly, birds were divided in two sampling groups representing both temperature treatments and selection lines (see above). Blood was sampled every other week from each individual between January and July. Every pair was sampled within 10 min after capture and then put back in their climate-controlled aviary. During a sampling day, birds (n = 36) were sampled between 8:30 AM and 14:30 PM. Individual birds were randomized by treatment and line to three sampling groups and individuals were sampled around the same time of day every time they were blood sampled. Red blood cells were separated from plasma by centrifuging at 14,000 rpm for 10 min and were stored in Queens buffer at room temperature until being processed. This treatment is sufficient to enrich the sample for red blood cells as the ratio between white to red blood cells in a whole blood sample is normally 1:800 (Vinkler et al. 2010). For the analysis, four time points were chosen based on realized first egg lay dates of the females as described in Mäkinen et al. (forthcoming): 1) the day when day light length >12 h (time point 1), 2) the day when 25% of the females from the warm environment had initiated laying (time point 2), 3) the day when 25% and 50% of the females from the cold and warm environment respectively had initiated laying (time point 3), and 4) the day when 50% of the females from the cold treatment had initiated laying. As our chosen time points did not always coincide with the days of blood sampling scheme, we chose blood samples that were sampled closest to (±7 days) the four chosen time points. This resulted in a total of 63 samples, because we were unable to take one blood sample for time point 4, because one female (warm treatment) was incubating her eggs during the sampling.

Sample Processing for Reduced Representation Bisulfite Sequencing

The detailed description of the sample processing can be found in Mäkinen et al. (forthcoming). In short, DNA was extracted from red blood cell samples using FavorPrepT M 96-well Genomic DNA Kit (Favorgen). Quantity and quality of the extraction was assessed with Nanodrop 2000 (Agilent Biotechnologies) and by gel electrophoresis by running the samples on a 1% agarose gel. Library preparation and sequencing were done at Roy J. Carver Biotechnology Center, University of Illinois at Urbana-Champaign, USA. Samples were individually barcoded and pooled together in four pools of 16 samples, thereby randomizing individuals, sampling days, and treatments for each pool. To obtain enough coverage, each pool was sequenced on two lanes for 100-bp single-end. To avoid effects due to difference in flow cells, the four pools, which were all ran twice, were ran on the same flow cell.

The control and bisulfite-treated libraries were prepared with the Hyper Library Construction kit from Kapa Biosystems, according to the “Reduced Representation Bisulfite Sequencing For Methylation Analysis” protocol (Illumina). Briefly, samples were digested with MspI, which generates CCGG overhangs, prior to bisulfite treatment and library building. The libraries were size selected for fragments of 20–200 bp in length. The libraries were quantified by quantitative PCR and pools were sequenced on a HiSeq2500 using a HiSeq SBS sequencing kit version 4 (Illumina). All pools were sequenced on the same flow cells. Fastq files were generated and demultiplexed with the bcl2fastq v2.17.1.14 Conversion Software (Illumina). Lanes were spiked with PhiX to provide nucleotide diversity and to measure bisulfite conversion success. PhiX reads were subsequently removed from the data set. Adaptors and individual barcodes were also trimmed from the reads. All the reduced representation bisulfite sequencing (RRBS) samples are deposited in Short Read Archive under Biosample accession SAMN07692587 and are described in detail in Mäkinen et al. (forthcoming).

Sequence Alignment and CpG Site Identification

The sequence processing and alignment has been described in detail in Mäkinen et al. (forthcoming) . In short, prior to alignment, the obtained sequences were adaptively trimmed and quality filtered to a minimum of 20 bp and a quality of 25 using Trim Galore! 0.4.2. Genome alignment and methylation site identification were performed with Bismark v0.16.3 using Bowtie2 v2.2.3 in the alignment step. Reads were aligned against the Parus_major1.1 genome (NCBI) using single-end mode and default settings in Bismark (L = 20 and N = 0). CpG sites were then identified from each sample; a site was accepted for downstream analysis when it had a minimum coverage of 10× and when it was identified across all samples. In total, the final data set for downstream analysis consisted of 522,645 CpG sites (supplementary table S2, Supplementary Material online).

We first performed hierarchical clustering on the final data set as well as for a reduced data set with only unrelated individuals using the “cluster” (Maechler et al. 2016) and “factoextra” (Kassambara and Mundt 2017) packages in R. The distance method used was Euclidian and the clustering method Ward.D2. We also calculated mean methylation values over the 522,645 CpG sites for the whole data set, and for each time point, each individual and for both treatments separately.

Methylation Patterns over Time

To examine temporal changes in DNA methylation across the four time points, we used GPrank (Topa and Honkela 2018). GPrank uses Gaussian process modeling to rank and identify temporal patterns and has previously been applied to both temporal transcriptome data (Topa and Honkela 2016) as well as temporal genotype data (Topa et al. 2015). By utilizing beta-binomial Gaussian process regression in a Bayesian framework, this allows for modeling any form of temporal variation in DNA methylation (linear as well as highly nonlinear dynamics) and is therefore well suited and frequently used to analyze data from sequencing-based genetic time series.

Our main interest was to detect changes in DNA methylation across the breeding season and between the two contrasting temperature treatments. In GPrank, the effect of time is modeled using a Gaussian process (Rasmussen and Williams 2006). For each CpG site, we fitted a time-independent model without a temporal change (model 0, i.e., null model) and two different time-dependent models, allowing for a temporal change in DNA methylation (model 1 and model 2; see below). Both the time-independent and time-dependent models included a fixed variance component from the beta-binomial model that captures the uncertainty due to finite sequencing depth (Topa et al. 2015), and an additive white noise, which represents random background noise. In the time-dependent models (model 1 and model 2), methylation levels are allowed to change between sampling time points, which was modeled using a squared exponential covariance in the Gaussian process: it tests correlation between two data points to estimate length scale and amplitude in the data (i.e., how fast the change is in time) to improve point estimate fit with the Gaussian process function. In the time-dependent model 2, we additionally allowed for differences between treatment groups in the pattern of temporal change (compared with model 1) by allowing the treatment groups to have separate time-dependent components and thereby capture possible treatment-specific temporal trajectories. We made comparisons between model 1 and model 0 to test for time dependence and between model 2 and model 0 to test for time dependence with treatment-specific trajectories and between model 2 and model 1 to test for treatment-specific trajectories of the time-dependent sites. Thus, sites were classified as time independent, time dependent, or as time dependent with treatment-specific trajectories (table 1). All three model comparisons were run separately for each CpG site and we selected the best fitting model based on the log-Bayes factor values (log BF), using log BF ≥ 3 which in the log-Bayes factor scale is considered as positive evidence in favor of the alternative model (Raftery and Kass 1995). We also clustered the data set using only time-dependent sites (without or with treatment effect).

Table 1

Model Selection Criteria and Classification of CpG Sites

Model Comparison Criteria	Classification of CpG Site
Model 1–Model 0 < 3 Model 2–Model 0 < 3 Model 2–Model 1 < 3	Time independent
Model 1–Model 0 > 3 Model 2–Model 0 < / > 3 Model 2–Model 1 < 3	Time dependent
Model 1–Model 0 < 3 Model 2–Model 0 > 3 Model 2–Model 1 > 3	Time dependent with treatment

Note.—Left-hand side describes the criteria for the model comparisons, based on log BF, which needed to be fulfilled for a CpG site to be classified as time independent, time dependent, or time dependent with treatment (right-hand side).

Because the GP models cannot take relatedness into account, we examined 5,000 sites across the log BF range of 49–3 and identified as temporally changing based on GPrank model 1 with GLMM using R packages lme4 (Bates et al. 2015) and lme4qtl (Ziyatdinov et al. 2017). We included the IBD matrix (supplementary table 2, Supplementary Material online) as a random effect in lme4qtl to control for potential differences in methylation due to shared sibship. Otherwise, both models used methylation count as the response, time point as a fixed effect and individual repeated measure was added as random effect (for both lme4 and lme4qtl). We compared the estimates between the two models for time point and intercept. As the estimates for both models, with and without controlling fore relatedness were highly congruent (supplementary fig. S1, Supplementary Material online), we conclude that relatedness per se did not influence our analysis when looking at temporal change in DNA methylation patterns. Ideally however, the analyses would be done in a GP model where both relatedness and the repeated measures structure of the data could be considered.

Genomic Locations and Gene Annotation

Genomic locations of identified CpG sites were annotated using the Great tit genome v1.1 (NCBI). We divided the sites into whether they were located in exons, introns, TSS regions (defined as the region spanning 300-bp upstream and 50-bp downstream of the annotated gene start), or promoters (defined as the region 3 kb upstream of gene start excluding the TSS) using the R packages “GenomicFeatures” (Lawrence et al. 2013) and “rtracklayer” (Lawrence et al. 2009). CpG islands were identified from the great tit genome using CpGplot v EMBOSS: 6.6.0.0 (Alan Bleasby, European Bioinformatics Institute, Wellcome Trust Genome Campus, Hinxton, Cambridge CB10 1SD, UK) using default settings: window size 100, minimum island length 200, minimum observed/expected 0.6, minimum percentage 50. Overlap between identified CpG sites and specified regions, was obtained with BEDtools v.2.26.0. (Quinlan and Hall 2010) and the sites were assigned to their specific genomic regions. We used the annotation from genomic locations (promoter, TSS, and within gene body as defined above) for assigning the CpG sites to genes for functional analysis, which was done using Cytoscape plugin ClueGo (Bindea et al. 2009) with chicken (EBI-Quick 31.1.2018) annotation. ClueGo identifies functional groups of genes based on Gene Ontology (GO) terms and enrichment analysis. As the reference set in the functional analysis and enrichment analysis, we used all genes for which at least one CpG site used in the GPrank analysis was associated with. We determined these genes based on the intersections of the CpG sites and the genomic locations described above (i.e., promoter, TSS, intron, and exon of a gene). To account for nonindependence of CpG sites (i.e., genes having many sites and sites supported by different models), genes were categorized to be time dependent or time dependent with treatment-specific trajectories based on what the CpG site with the highest ranking log BF value associated with the gene was classified to be. For a gene to be categorized as time independent, all CpG sites associated with that gene, regardless of the genomic region, had to be classified as time independent.

Results

After filtering, 522,643 sites, covering 3.4% of the known CpG sites in the great tit genome, were available for downstream analyses (supplementary table S1, Supplementary Material online). The genome-wide methylation level across all the sample time points and individuals for these sites was 13.82% (SD = 0.25). Methylation levels were significantly lower in individuals from warm than in cold treatment in the first three time points (t = −4.3186, df = 8,362,200, P value = 1.57e-05; t = −4.154, df = 8,362,100, P value = 3.268e-05, and t = −3.6373, df = 8,362,300, P value = 0.0003, respectively, fig. 1), although the difference was small (delta = −0.0008, delta = −0.0007, and delta = −0.0007 for time points 1, 2, and 3, respectively). For the fourth time point, there was no difference in genome-wide DNA methylation levels (delta = −0.0002, t = −1.3819, df = 7699900, P value = 0.167, fig. 1). In general, we observed a slight increase in mean methylation level over time and this pattern was similar in both the cold and the warm treatments (ANCOVA: F2,5 = 94.09, P value < 0.001, fig. 1). Individual level variation in mean methylation level (supplementary fig. S2, Supplementary Material online) as well as in how methylation changed across time was large (supplementary fig. S3, Supplementary Material online) but hierarchical clustering using all 522,643 CpG sites demonstrated that an individual’s replicated measures were a better predictor of similarity compared with treatment (fig. 2 and supplementary table S2, Supplementary Material online). As a consequence of using selection lines, and thus a limited number of individuals, there are nine full-sibs in the data set (three sister pairs and one sister triplet, supplementary table S1, Supplementary Material online). Removal of these siblings (see Materials and Methods) from the data did not change the replicate driven clustering (supplementary fig. S4A, Supplementary Material online).

. 1.

—Distributions of methylation levels in the analyzed sites in the data. (A) Mean methylation levels across time points per treatment. Error bars represent standard error of the mean. t-Test P values between groups were 1.57e-05, 3.268e-05, 3.268e-05, and 0.167 for time points 1, 2, 3, and 4, respectively, when calculated across all analyzed CpG sites. (B) Histogram of mean methylation level across time points per CpG site.

. 2.

—Cluster dendrogram for the analyzed sites. Sample names (ring number) with indication of the corresponding sampling time point (indicated by _1 to _4 after the ring number). Sample names are colored based on their corresponding temperature treatment with red for warm and blue for cold. The four clusters containing full-sisters are indicated with curly brackets on the right hand side of the figure.

Not surprisingly, most CpG sites were lowly methylated (median 1.13%) although methylation levels spanned the entire range from completely unmethylated to fully methylated (fig. 1). Considering the change in methylation across time points per CpG site (i.e., absolute value of the difference between largest and the smallest methylation level within a site), we observe large variation given the minimum level of methylation within a site (supplementary fig. S3, Supplementary Material online). Mean individual methylation change across all sites was 6.82% (SD = 8.45) and mean methylation change within a site was 18.7% (SD = 18.17). Between treatments, methylation change across time in the CpG sites was slightly lower for warm compared with cold treatment with mean 6.80% (SD = 8.50) versus 6.84% (SD = 8.40) for warm and cold, respectively (delta = 0.04, t = 3.1247, df = 1045100, P value = 0.0018). Thus, in general, there were relatively minor differences between treatments and time points but substantial between-individual variation.

Temporal Variation in DNA Methylation

In total, we found that 47,922 CpG sites out of 522,643 analyzed sites showed support for a time-dependent methylation pattern (i.e., statistical support for time dependence or time dependence in a treatment-specific manner) and these represent 9% of all the analyzed sites (table 2 and supplementary table S3, Supplementary Material online). We repeated the clustering using only these sites and found that the samples formed two main clusters based on treatment (supplementary fig. S4B, Supplementary Material online). The effect of the repeated measurements is still visible as samples of the same individual are clustered together. In a similar manner, the effect of siblings split into both treatment groups is still visible as the siblings from the warm treatment cluster together with their cold treatment siblings in to the cold treatment cluster (supplementary fig. S4B, Supplementary Material online). Thus, for all sites as well as only the time-dependent sites, the main similarity among samples is the repeated sampling of individuals across time.

Table 2

Characteristics of the Analyzed CpG Sites Classified by Their Time Dependency

	Time Independent	Time Dependent	Time-Dependent Treatment-Specific Trajectory
CpG sites assigned to different models	474,723	37,459	10,463
Mean methylation (SD)	0.13 (0.24)	0.17 (0.30)	0.33 (0.33)

Note.—Number of CpG sites classified to each group as well as mean methylation level (and standard deviation) is reported.

Out of the 47,922 sites that showed a temporal pattern, the majority, 37,459 sites, were classified as time dependent without a treatment-specific effect, whereas 10,463 sites only showed a time-dependent effect when the treatments were allowed to have separate trajectories (table 2, supplementary table S3, Supplementary Material online, and fig. 3).

. 3.

—Manhattan plot of statistical support for temporal change in DNA methylation on a per CpG site basis using the Gaussian Process modeling approach, black line at log BF = 3 indicate evidence of good statistical support (Raftery and Kass 1995).

As expected, there was a positive correlation between the number of analyzed CpG sites (corrected for chromosome length) and GC content of the chromosome (r = 0.84, df = 32, P < 0.01, supplementary fig. S5A, Supplementary Material online). This is reflected also in the numbers of time-dependent sites in the different chromosomes. Chromosomes 25LG1 and 25LG2 have the largest number of CpG sites assigned to them in relation to their size (supplementary fig. S5B, Supplementary Material online) as the GC content of these chromosomes is the highest in the whole Great tit genome (Laine et al. 2016). Similarly, GC content is slightly higher also in chromosomes 26, 27, and 28 and LGE22 compared with the other chromosomes (supplementary fig. S4A, Supplementary Material online). This pattern of high GC content is typical for the great tit microchromosomes. In the mitochondria, we observed five sites and these were all time dependent and intragenic (supplementary table S3, Supplementary Material online).

Most of the analyzed CpG sites were intronic (24.8%), whereas for TSS, exon, intergenic, and promoter, the proportions were quite similar (18.3%, 17.3%, 17.4%, and 22.1%, respectively). Methylation levels in the different genomic regions also differed significantly (Wilcox test, see supplementary table S5, Supplementary Material online) and methylation levels were highest in exons and intergenic regions (supplementary fig. S6, Supplementary Material online). The proportions of the genomic regions within time-independent CpG sites are similar to the genomic regions for all analyzed data set with 19.45%, 17.1%, 16.5%, 24.6%, and 22.2% for TSS, exon, intergenic, intron, and promoter, respectively. For the time dependent and time dependent with treatment effect the corresponding values are 7.4%, 18.6%, 25.3%, 27.0%, and 21.6% and 4.3%, 20.6%, 28.7%, 26.8%, and 19.5%, respectively, for TSS, exon, intergenic, intron, and promoter. Methylation levels were higher for CpG sites classified as time dependent with treatment effect (table 2) and this was also seen when splitting the CpG sites based on genomic location and model (supplementary fig. S7, Supplementary Material online). In addition, methylation levels did not differ between the treatments for sites classified as time dependent or time independent, whereas levels differed for sites classified as the time dependent with treatment effect as expected (supplementary fig. S7, Supplementary Material online).

There were differences between genomic regions of the time-dependent CpG sites (47,922 in table 2.) and the time-independent sites. The number of sites annotated to TSS is lower in time-dependent sites (both without or with treatment), whereas more sites are annotated to exon or intergenic regions in the time-dependent class, especially for time-dependent treatment (fig. 4, Pearson’s chi-square test, df = 8, P < 0.001).

. 4.

—Genomic regions of CpG sites in the time-independent and the two time-dependent classes. Colors represent the different genomic locations based on annotation of the Parus major genome v1.1.

CpG sites classified as time dependent (with or without treatment) overlapped 41% of the annotated genes in Parus major v1.1 genome; that is, 53% of the genes which were used in the analysis (i.e., the reference set, supplementary table S4, Supplementary Material online, considering promoter, TSS, and gene body). Genes had on average 15.00 CpG sites assigned to them (SD = 14.12) and typically included CpG sites classified as both time dependent and time dependent with treatment effect (supplementary table S4, Supplementary Material online).

Functional analysis, which is grouping genes based on GO-term enrichment (BH corrected P value <0.05), for time-dependent classified genes identified 106 functional groups mainly involved in regulation and development of different tissue types (supplementary tables S6, Supplementary Material online), whereas for time-dependent treatment classified genes, there were six functional groups identified and these were involved in chromatin modification related functions (GO:0042800 and GO:0048188) as well as transcriptional coactivation (GO:0001223, supplementary table S7, Supplementary Material online). In the time-dependent genes, we identified specific GO-terms not only for developmental functions in different tissues but also for energy balance and reproductive system development (GO:0003006, GO:0048608, GO:0061458, GO:2000243, GO:0008406, GO:0045137, and GO:0046660). Genes underlying these latter GO-terms as well as GO-terms in the six groups from the time dependent with treatment-specific trajectories can be regarded as potential candidates for timing of breeding in the great tit (supplementary tables S6 and S7, Supplementary Material online).

When looking at genes in more detail, we found that DIO2, a candidate gene for timing of reproduction in the great tit (Perfito et al. 2012), showed temporal change in DNA methylation. Two out of the eight CpG sites observed in this gene in our data set were classified as time dependent (supplementary table S4, Supplementary Material online). These sites are located in the first exon of the gene (UTR). The overall methylation level of the eight sites in DIO2 have mean methylation level of 2.4% (SD = 1.5%), and the nonlinear change over time is low in general (supplementary fig. S8, Supplementary Material online). One of the new potential candidate genes identified with the GO analysis, RORA (supplementary table S7, Supplementary Material online, group 1), exhibits a site-specific methylation pattern, which also differs between the treatment groups across time (supplementary fig. S9, Supplementary Material online). It has 30 CpG sites associated to it of which three are time dependent and nine time dependent with treatment. All the sites are located in introns. In contrast to DIO2, methylation levels in RORA are higher for the methylated sites (mean = 23.6%, SD = 27.2%) making the difference between methylated and nonmethylated sites more evident than in DIO2. Both of these genes also indicate linear and nonlinear changes across time.

Discussion

Here, we examined within-individual temporal changes in DNA methylation in great tit females across the breeding season by analyzing repeated methylation measures using a RRBS approach. We were able to sample 522,643 CpG sites with >10× coverage in all samples and 47,922 of these were found to exhibit time-dependent change either without or with treatment-specific trajectories. In general, the methylation levels and temporal changes were slightly higher in the cold compared with the warm treatment. There is a strong temporal component to DNA methylation patterns across the breeding season as 9% of all sites displayed a temporal pattern and thus changes in DNA methylation at some of these sites could potentially be involved in regulating the onset of reproduction.

Temporal Patterns in DNA Methylation

Considering that DNA methylation is viewed as a mechanism for gene regulation, most ecological studies examining DNA methylation have only sampled single time points, leaving us with a poor understanding of the temporal stability of DNA methylation (but see e.g., Rubenstein et al. 2016; Saino et al. 2017). To our knowledge, this is the first animal study to examine genome-wide DNA methylation patterns using a repeated measures approach, which up till now has only been conducted in few clinical studies in humans (Shvetsov et al. 2015; Simpkin et al. 2015; Leenen et al. 2016; Urdinguio et al. 2016). Even in these studies, CpG methylation is typically measured between two time points only, which hinders detection of nonlinear changes in methylation. One exception is (Simpkin et al. 2015) who measured CpG methylation at three time points (at birth and at ages 7 and 17). They found evidence of both linear and nonlinear patterns in methylation levels of CpG sites associated to birth weight and gestational age across the sampling period. We measured CpG methylation at four time points between March and June and observed both linear and nonlinear changes, as seen for example with DIO2 and RORA, within this short time period, thereby highlighting the dynamic nature of DNA methylation patterns.

Previous studies using repeated measures in humans have identified from 0.5% to 30.1% of analyzed sites changing methylation patterns over time (Wang et al. 2012; Martino et al. 2013; Urdinguio et al. 2016). In this study, we found that 9% of the analyzed CpG sites changed their methylation patterns over time, with an average individual methylation change of 2%. For comparison, Martino et al. (2013) reported that in human twins, methylation pattern between birth and 18 months had a mean methylation change of 3.1%. A recent review on the role of DNA methylation in gene regulation in clinical context reported that even low methylation changes at a site in the range of <10% can have a phenotypic effect in nonmalignant diseases (Leenen et al. 2016).

A substantial challenge is to understand how much DNA methylation should change in order to have a phenotypic effect. In our case, we often observe that methylation change occurs only in specific time points and in a single, or a few, CpG sites within a gene. Clinical studies on nonmalignant diseases have implicated that even single sites could be important for creating a disease phenotype (Leenen et al. 2016). In an ecological context, effects of changes in methylation at single CpG site have for example been found for ERα and yolk testosterone in eastern bluebird (Bentz et al. 2016) and for SERT and exploration behavior in urban great tits, although in great tits there was also SNP variation associated with it (Riyahi et al. 2015). In contrast, in another study on great tits and exploratory behavior, only variation in DNA methylation in sites within the TSS of the DRD4 gene was associated with the phenotype, whereas this was not the case for the sites sampled in the gene body (Verhulst et al. 2016).

Genome-wide methylation patterns correlated to transcriptome expression from the great tit show that CpG methylation was negatively correlated with gene expression in gene bodies and at TSS in brain tissue (Derks et al. 2016; Laine et al. 2016). When both promoter and gene body CpG sites were hypomethylated, expression in the brain was higher (Derks et al. 2016). Furthermore, differential CpG methylation between brain and whole blood (including both erythrocytes and leucocytes) in these regions had a functional role which was shown by differential patterns of hypo- and hyper-methylation between the tissues in regard to expression patterns in the brain (Derks et al. 2016). However, although these studies find methylation and expression are correlated, we do not yet have a good understanding of what level of change in methylation, or in which CpG’s in a region, is sufficient to change transcription and thus lead to a phenotypic response. Thus, although it is clear that the environment can leave a mark on an individuals’ DNA methylome (Bentz et al. 2016; Rubenstein et al. 2016; Weyrich et al. 2016; Gokhman et al. 2017; Pértille et al. 2017), we need more studies to establish potential causative effects of such methylation changes on the phenotype of organisms.

Genomic Locations of Time-Dependent Sites

The majority of time-dependent (with or without treatment) CpG sites were annotated to introns and exons rather than to promoter or TSS. Although this result is different to that reported in other RRBS data sets in nonmodel organisms (measured at single time point only) (Pegoraro et al. 2016), this is probably a result of differences in the type of tissue typically examined in different organisms: whole body versus red blood cells. DNA methylation data from red blood cells in chicken also had a pattern similar to this study, that is, with gene region methylation being more abundant than regulatory region methylation although measured only at a single time point (Pértille et al. 2017). Although more CpG sites were annotated to introns than in exons in sites categorized as time dependent (with or without treatment), the measured methylation levels in the genomic regions were higher in gene bodies compared with promoter regions with very low methylation levels around the TSS, as expected based on previous great tit methylomes (Derks et al. 2016; Laine et al. 2016). Although more of the CpG sites that showed time-dependent patterns were situated within introns compared with exons, methylation level in introns was moderate compared with exons. Both of these observations on intronic CpG sites fit with recent knowledge of intronic CpG sites in regard to regulating gene expression as part of the promoter (Suzuki and Bird 2008). As an example of a promoter including also the first intron, pyrosequencing of the DRD4 gene in great tit selection lines for exploratory behavior identified most of the variation in DNA methylation within in the first intron, whereas methylation level increased toward the 3′ end of the gene (Verhulst et al. 2016).

Given that we are mapping genomic locations of time-dependent sites, it is likely that they exhibit a pattern that differs from the usual expectation of identifying differentially methylated sites in TSS or promoter regions. In an investigation of human leukocyte CpG methylation patterns between birth and age at 5, over 6,000 CpG sites were observed changing (Urdinguio et al. 2016): the sites identified as hypomethylated over time were related to immune system functions and were located in introns, whereas hypermethylated sites were related to development and sequence-specific binding and were located in exons. The time-dependent (without or with treatment effect) CpG sites in the red blood cells in our study were also seen in introns and exons, whereas we observed less than expected sites in TSS compared with other genomic regions. Thus, our result are in line with what is observed in other species where time-dependent changes accumulate in a different manner in different genomic regions (Martino et al. 2013; Urdinguio et al. 2016).

Temporal Change in DNA Methylation and Its Potential Role in Timing of Breeding

Our study shows that there is potential in using reversible methylation patterns to identify genes responding to environmental cues similar to transcriptomic studies. In particular, we can compare the genes that displayed temporal changes in methylation with those that additionally showed a temperature effect. For example, although the GO-terms in the time-dependent class were mostly involved in developmental and regulatory effects in various tissues, the time-dependent treatment class was involved in chromatin modification, more specifically H3K4 methylation (GO:0042800 and GO:0048188). This chromatin mark has previously been associated with transcriptional activation (Muramoto et al. 2010; Liu et al. 2017). Another group (GO:0001223; transcription coactivator binding) contained three genes ESR1, POU2F1, and RORA which provides a potential link between temperature sensing and reproduction. In ESR1, the estrogen receptor, methylation at the promoter was recently shown to be a prominent mechanism that mediates maternal environmental effects (Bentz et al. 2016). RORA is a gene with circadian properties (Jetten 2009) and POU2F1 is a transcription factor connected to progesterone and glucocorticoids (Zhao 2013).

We did not see time-dependent methylation changes in many of the typically investigated candidate genes in relation to timing of reproduction (e.g., BMAL, CLOCK, PER2/3, and TSH) in our data set, with the exception of DIO2 (Perfito et al. 2012). There can be several reasons for this, sampling the wrong tissue (hypothalamus is generally preferable for seasonal timing related genes), missing the correct time points or not knowing enough of the methylation dynamics across life time, as discussed below. When looking at DIO2, which exhibited a vastly different methylation pattern from RORA, it is possible that we did not sample the right time window for it, or for many other genes, as it is likely that the genes facilitating the photoperiodic response have established their methylation pattern already prior to our first sampling day at the end of March (Visser et al. 2010; Stevenson and Prendergast 2013). It is also plausible that certain genes have stable methylation patterns throughout a lifetime as a result of tissue-specific imprinting (Lokk et al. 2014). Two recent studies on birds measured DNA methylation from blood in single genes twice and found no difference over time in the methylation levels (Rubenstein et al. 2016; Saino et al. 2017). In the case of superb starlings (Lamprotornis superbus), glucocorticoid hormone receptor CpG methylation was measured when the individuals were chicks and later on as adults to monitor the temporal pattern and with this sampling scheme they were able to show that early life environmental conditions (i.e., rainfall) had a sex-specific impact on life-history traits (Rubenstein et al. 2016). In a study on migration phenology in barn swallows (Hirundo rustica), individuals were measured twice as adults in two consecutive years in the breeding range (i.e., after arrival from wintering grounds). This raises the question on whether the observation of stable methylation in CLOCK poly-q (Saino et al. 2017) is due to not sampling the right time frame or due to the methylation pattern being stable from early life onward. Thus, predictions of the outcome of methylation change for timing of breeding in great tits and the underlying gene regulation of these functional groups is highly tentative until the link between observed methylation changes and transcriptomic expression has been tested and confirmed in additional tissues (in particular the brain).

Limitations and Caveats

As discussed above, we acknowledge that red blood cells are not the optimal tissue to understand gene regulatory changes related to decision making and timing of reproduction in the great tit. However, it is not possible to do repeated measures in more relevant tissues, that is, in the hypothalamus–pituitary–gonad–liver axis, which are more directly associated with functions important for timing of reproduction (Visser et al. 2010). Although tissue specific DNA methylation patterns develop early on (Li et al. 1992; Bird 2002), DNA methylation is also accumulated throughout life at CpG sites (Christensen et al. 2009; Khor et al. 2016; Metzger and Schulte 2017; Pértille et al. 2017) and thus whether environment driven methylation is tissue specific at first, is not well known. Moreover, the correlation in DNA methylation patterns among different tissue can be high (Lokk et al. 2014). Thus, using a general tissue such as red blood cells can be informative also of methylation patterns in other tissues.

Although our sampling design covers a large proportion of the reproductive season it is of course possible that we may have not chosen the right sampling times to detect changes in methylation relevant to timing of reproduction. For example, it may be that the regulatory changes necessary for egg-laying take place over a longer or shorter timespan relative to laying of the first egg (Williams 2012) than what we can detect with our sampling interval. Thus, identifying the reference points for environmental cues and egg-laying would allow for better resolution when designing the sampling scheme.

One strength of our design is the use of individuals from selection lines housed in climate-controlled aviaries. Individual level variation in DNA methylation patterns is high and thus limiting environmentally induced variation in DNA methylation can increase our ability to detect changes relevant to timing of reproduction as all individuals have experienced similar environmental conditions (apart from the treatment effect of temperature of course). In this study, we only examined DNA methylation from the early selection line and from a single generation (F2). Thus, we are not in a position to address potential effects of selection on the DNA methylation patterns between the selection lines or changes across generations that may have occurred during the selection experiment. Thus, to what extent selection has induced changes in the methylation landscape compared with the “ancestral” (wild) great tits is not clear.

Our main focus here was on the temporal stability and effects of the temperature treatment on DNA methylation patterns and only indirectly linked to the observed timing of reproduction. We are currently working on taking another analytic approach to more directly identify differently methylated sites and regions between early and late reproducing individuals within the temporal sampling design we have. Thus, the link to timing of breeding for the sites and genes identified here of course remain tentative.

Lastly, as we only measured DNA methylation in these samples but not RNA expression we cannot draw firm conclusions of the effects of single sites to gene expression. Although earlier studies have examined how DNA methylation in different genomic regions affects gene expression in the great tit (Derks et al. 2016; Laine et al. 2016; Verhulst et al. 2016), they do not have a temporal expression data set to compare with temporal methylation. This is the next important step toward verifying whether the temporal trends in methylation are functional.

Conclusions

Here, we investigated temporal within-individual changes in erythrocyte DNA methylation in great tits of genomic selection line for early timing of breeding which were exposed to either a warm or cold treatment. We find that many CpG sites changes methylation levels over the sampling period, but that the overall change is small. There was also substantial heterogeneity across CpG sites within genes in their response to time or treatment. Thus, our work calls for a greater understanding of the relationship between changes in methylation and expression at individual level and across different tissues. Future studies that examine this relationship across time and on the same individual would be very valuable.

There is currently great interest in ecological epigenetics and a number of studies on primarily DNA methylation patterns in natural populations have been published over the last few years. However, most studies only sample a single time point and a few targeted CpG sites. Our results highlight that there is large between individual variation in DNA methylation and a strong temporal component to DNA methylation patterns. This, combined with the heterogeneity observed also within genes, suggest that great caution is required when interpreting studies that have sampled few sites at single time points.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.