Sex-specific changes in autosomal methylation rate in ageing common terns

1. Introduction

Senescence is a within-individual decline in survival probability (actuarial senescence) and/or reproductive performance (reproductive senescence) with age. Although the rate and shape of the decline vary both among and within species, this detrimental process occurs in most species across the tree of life (Shefferson et al., 2017). It is hypothesized to have evolved because unavoidable environmentally-driven mortality reduces the strength of selection against poor performance with age (Fisher, 1930; Medawar, 1952; Williams, 1957; Hamilton, 1966). From a genetic point of view, this so-called ‘selective shadow’ could allow for the accumulation of late-acting deleterious mutations over evolutionary time (mutation accumulation, Medawar, 1952), or for selection favouring alleles with beneficial effects early, but detrimental effects late in life (antagonistic pleiotropy, Williams, 1957). Moreover, early-life investment of limited resources in reproduction over that in perfect somatic maintenance and repair would also be evolutionarily beneficial, such that senescence could be a ‘best-of-a-bad-job’ consequence of accumulated, unrepaired damage (disposable soma, Kirkwood, 1977).

Although the classical theories of mutation accumulation and antagonistic pleiotropy assume a purely sequence-level genetic basis to senescence, the disposable soma theory does not. Given its firm foundation in life history theory, resource acquisition and allocation (sensu van Noordwijk and de Jong, 1986) are important facets to consider. Both are known to show phenotypic plasticity in response to environmental variation (e.g., Erikstad et al., 1998; Descamps et al., 2016), and experimental manipulation of one of these facets without adjustment of the other is known to affect rates of senescence (e.g., Boonekamp et al., 2014). As such, part of the senescence process is also expected to be underpinned by more flexible, regulatory processes (Wilson et al., 2008). In line with this, numerous studies across taxa have shown that senescence is underpinned by both genetic and epigenetic processes (Sen et al., 2016; Melzer et al., 2020).

Epigenetic processes are those that affect the regulation of gene expression (Holliday, 2006). This regulation is complex and encompasses several mechanisms, a major one of which is DNA methylation, the addition of a methyl-group to the fifth carbon site of a cytosine in a CpG (5’-C-phosphate-G-3’) dinucleotide context (Jaenisch and Bird, 2003; Miranda and Jones, 2007; Suzuki and Bird, 2008). Such DNA methylation could lead to transcriptional silencing of genes and repetitive elements, as DNA methylation acts as a ligand for methyl-binding domains, which in turn are targets for further chromatin-modifying complexes such as histone deacetylase complexes or chromatin remodeling complexes (Jaenisch and Bird, 2003; Suzuki and Bird, 2008; Brenet et al., 2011; Miller and Grant, 2013). DNA methylation of gene bodies, however, can lead to increased gene expression (Suzuki and Bird, 2008) or be involved in alternative splicing (Lev Maor et al., 2015). Age-specific DNA methylation has been described for many model species (e.g., Maegawa et al., 2010; Hannum et al., 2013; Tharakan et al., 2020), and the methylation status of specific CpGs has been shown to be a powerful predictor of both chronological and biological age, i.e., to function as an epigenetic clock (Bocklandt et al., 2011). Although age-specific methylation therefore seems the norm, at least in model organisms, the direction of methylation changes with age is harder to predict: both non-directional changes (representing epigenetic drift, Fraga et al., 2005; Tan et al., 2016), and global directional changes (representing global hyper- and hypomethylation, Zampieri et al., 2015; Ciccarone et al., 2018; but see Unnikrishnan et al., 2019) have been reported.

Research on senescence, including its epigenetic underpinning, has so far mostly focused on humans or model organisms kept under controlled laboratory conditions. Extending the taxonomic range and incorporating field studies is crucial to understand the evolutionary ecology of senescence (Monaghan et al., 2008). For this extension, birds are an especially interesting taxon. They have longer life spans relative to their body size than mammals (e.g., Lindstedt and Calder, 1976), and various populations of birds with vastly different life histories have been studied over several decades, such that many basic insights into their senescence patterns can now be obtained (Bouwhuis and Vedder, 2017). Bird genomes are rather compact, show high levels of synteny across species (Zhang et al., 2014), and, similar to those of other vertebrates, are globally methylated (Suzuki and Bird, 2008; Li et al., 2011; Sepers et al., 2019). With respect to age-specific DNA methylation in birds, we are aware of the existence of an epigenetic clock for short-tailed shearwaters (Ardenna tenuirostris; De Paoli Iseppi et al., 2019), and of findings of early-life within-individual increases in global methylation in both great tit (Parus major; Watson et al., 2019) and zebra finch (Taeniopygia guttata; Sheldon et al., 2020) nestlings. However, we are not aware of any study reporting age-specific DNA methylation patterns in avian adulthood, or in avian late-life specifically.

When studying age-specific trait expression, many studies, including most of those producing epigenetic clocks or otherwise studying age-specific differences in DNA methylation patterns, use cross-sectional samples and analysis tools. Patterns revealed by cross-sectional approaches, however, are the result of a combination of within- and among-individual processes. If we aim for understanding the within-individual process of senescence, we need to account for the effect of compositional changes of a population, for example when birds with a certain level of DNA methylation are more likely to die and selectively disappear from the study population (e.g., Vaupel and Yashin, 1985; Forslund and Part, 1995). Mixed-effect models applied to (partly) longitudinal data to specifically test whether within-individual patterns and population-level patterns are the same, or differ, are a powerful analytical tool to do so (van de Pol and Wright, 2009). Across taxa, efforts have increasingly been made to validate and complement findings from cross-sectional analyses of senescence through longitudinal studies (Nussey et al., 2008; Bouwhuis and Vedder, 2017; Gaillard et al., 2017), and this methodological turn is also reflected in studies of human DNA methylation (Bollati et al., 2009; Tan et al., 2016). Longitudinal studies of DNA methylation and ageing in model species and natural populations are, however, rare (but see Lemaître et al., 2021) and needed (Bell et al., 2019).

Here, we report on a longitudinal study of autosomal methylation levels in the common tern (Sterna hirundo). Common terns are long-lived migratory seabirds whose patterns of senescence have been the topic of various studies. Although breeders of both sexes show little sign of reproductive senescence - they breed earlier in the year and fledge more offspring as they grow older (Nisbet et al., 2002; Zhang et al., 2015c; Nisbet et al., 2020) - breeding and survival probabilities are known to decline with age (Zhang et al., 2015b; Vedder et al., 2021b). In addition, there is evidence for sex-specific transgenerational senescence, with daughters of older mothers and sons of older fathers suffering from reduced lifetime reproductive success (Bouwhuis et al., 2015). Studies aiming to identify a molecular basis for these within- and transgenerational effects have so far focused on telomere dynamics, and found that: (i) telomeres shorten with age (Bichet et al., 2020); (ii) telomere length is genetically correlated with lifespan (Vedder et al., 2021a); and (iii) paternal age is negatively correlated with offspring telomere length (Bouwhuis et al., 2018). The explanatory power of telomere length, however, is very low for all of these patterns (e.g., 1.1% of phenotypic variation in lifespan is explained by additive genetic variation in telomere length, Vedder et al., 2021a), such that additional mechanisms are expected. To evaluate whether age-specific changes in global DNA methylation could be such a mechanism, we (i) sequenced and de novo assembled a high-quality chromosome-scale reference genome and (ii) used it to compare within-individual age-specific changes in DNA methylation at shared sites across the genome. Although rates of within-generational senescence do not differ between the sexes (Zhang et al., 2015a), transgenerational effects are known to be sex-specific (Bouwhuis et al., 2015), such that we also considered sex-specificity of any patterns in autosomal methylation status.

2. Materials and methods

2.1. Study population

The data we present were collected as part of a long-term individual-based study of a mono-specific common tern colony located on six artificial concrete islands at the Banter See (53°30’40” N, 08°06’20” E) in Wilhelmshaven, Germany. Fledglings from this colony have been marked with metal rings since 1984, subcutaneously injected with transponders since 1992 (Becker and Wendeln, 1997) and molecularly sexed since 1996 (Becker and Wink, 2003). This marking and sexing, combined with placement of antennae on elevated platforms on the edges of the colony site and around each nest during incubation (shared between both parents), allows for a well-described family structure of all known-sex and known-age philopatric birds (e.g., Moiron et al., 2020).

2.2. Reference genome

To create the chromosome-scale reference genome, blood of one adult female common tern was collected in 100% EtOH and stored at −80°C. 30μg of High Molecular Weight DNA (HMW DNA) was isolated from the whole blood sample using an agarose plug protocol of the Bionano Prep Blood and Cell Culture DNA Isolation Kit (cat no. RE-130-10) modified for avian nucleated erythrocytes. Lysates were embedded into agarose plugs, followed by Proteinase K and RNase A treatments and 1X TE drop dialysis purification. Four sequence datasets were generated following the VGP 1.5 pipeline (Rhie et al., 2021): 67.91x Pacific Biosciences (Pacbio) continuous long reads (CLR); 698.35x Bionano Genomics optical maps; 169.30x 10X Genomics linked-reads; and 79.62x Arima Genomics Hi-C Illumina reads.

To create the Pacbio data, DNA was sheared using a 26G blunt end needle (Pacbio protocol PN 101-181-000 Version 05) to approximately ~40kb fragment length. We used 10μg of this fragmented DNA to generate a large-insert Pacbio library using the Pacific Biosciences Express Template Prep Kit v1.0 (#101-357-000). The library was then size selected (>15 kb) using the BluePippin system (Sage Science). The resulting PacBio Library was sequenced on 10 PacBio 1M v3 smrtcells (#101-531-000) on a Sequel instrument with the sequencing kit 3.0 (#101-427-500) and a 10 h movie with 2 h pre-extension time. Unfragmented HMW DNA was used to generate a linked-read library on the 10X Genomics Chromium (Genome Library Kit and Gel Bead Kit v2 PN-120258, Genome Chip Kit v2 PN-120257, i7 Multiplex Kit PN-120262). We sequenced this 10X library on an Illumina Novaseq S4 150 bp PE lane. uHMW DNA was labeled for Bionano Genomics optical mapping using the Bionano Prep Direct Label and Stain (DLS) Protocol (30206E) and 1 flow cell was run on the Saphyr instrument. Hi-C libraries were generated with the Arima Genomics v1.0 2-enzyme protocol (P/N: A510008), according to the manufacturer’s protocol and sequenced on Illumina HiSeq X.

The resulting four data types were processed using the VGP v1.5 pipeline (Rhie et al., 2021), which includes: assembling Pacbio contigs using FALCON v2018.31.08-03.06; FALCON-Unzip v6.0.0.47841; purging false haplotype duplications with purge_haplotigs v1.0.3+ 1.Nov. 2018; scaffolding with 10X with scaff10x v4.1.0; scaffolding with Bionano Solve DLS v3.2.1; scaffolding with Hi-C data with Salsa HiC v2.2; filling in gaps and polishing for base call accuracy with CLR and Arrow smrtanalysis v6.0.0.47841; and polishing with Illumina short reads with longranger align v2.2.2; and freebayes v1.3.1. The resulting assembly was then manually curated to fix any errors, using gEVAL and Hi-C short read linked-read mapping profiles as described in Howe et al. (2021). BUSCO v4.1.4 with the bird lineage dataset (aves_odb10) was used to assess assembly completeness. The reference genome was submitted to NCBI as GCA_009819605.1, as part of the Vertebrate Genomes Project (VGP).¹ Reference genomes consisting of short-read assemblies (e.g., using Illumina reads) exhibit GC bias, where GC-rich regions such as promoters could be incorrectly assembled or even missing (Kim et al., 2021). Our reference genome, however, has GC-rich promoter regions due to the use of Pacbio long reads (Kim et al., 2021; Rhie et al., 2021).

2.3. Blood sampling and methylation sequencing

Blood of breeding common terns was sampled in the years 2013, 2014 and 2017 using larval stages of the bloodsucking bug Dipetalogaster maximus. Within 9-14 days of clutch completion of each focal bird in each year, bugs were placed into dummy eggs with holes and placed in the nests (Becker et al., 2006; Arnold et al., 2008; Bichet et al., 2019). After 20–30 min of incubation by the focal bird, “bug eggs” were collected and the focal birds’ blood, sucked by the bug, was removed from the bugs’ abdomen using a syringe. Upon collection, whole blood was stored in EDTA buffer (2%) in a fridge (3–7°C) for up to 3 weeks, before the red blood cells were transferred to glycerol buffer (40%) and frozen at −80°C. Each bug was used only once to prevent cross-contamination of blood samples.

In total, we obtained 74 samples from 34 individuals making up 17 breeding pairs: 3 breeding pairs (i.e., 6 individuals) were sampled in all 3 years (i.e., with intervals of 1, 3 and 4 years), 14 breeding pairs (i.e., 28 individuals) in two of the 3 years (with a 1-year interval for 4 breeding pairs, a 3-year interval for 9 breeding pairs and a 4-year interval for 1 breeding pair; see Supplementary Table S1 for details). The age of the sampled females ranged between 4 and 20 years, that of males between 3 and 19 years.

Genomic DNA was extracted from each sample and libraries for RRBS were generated as described in Klughammer et al. (2015), following standard steps such as MspI digestion, end-fill-in, A-tailing and size selection. This included the enrichment of the libraries with Pfu-Turbo Cx Hotstart DNA polymerase to allow the assessment of the bisulfite conversion rate. After clean-up and quality control, libraries were sequenced with an Illumina HiSeq 4000 (50 bp SE).

2.4. Methylation data selection

Across all 74 samples, individuals were sequenced with an average of 51,482,011 (range: 24,942,071–82,760,425) reads (Supplementary Tables S1, S2). After conversion of unmapped bam files to fastq with SamToFastq.jar from picard v1.118, the quality of the sequencing reads was checked using Fastqc v0.11.5 (Andrews, 2010) and Multiqc v1.8 (Ewels et al., 2016). Spiked-in sequences were used to estimate over- and underconversion using RefFreeDMA (Klughammer et al., 2015). Trimgalore v0.3.3 (implemented in RefFreeDMA; Krueger et al., 2021) was used to remove adaptors to guarantee good mapping and to remove low quality bases (≥20) as well as short read fragments (≥16 bp).

We used Bismark v0.22.3 (Krueger and Andrews, 2011) to prepare and index the reference genome for subsequent mapping of the bisulfite-converted reads using Bowtie 2 v2.4.1 (Langmead and Salzberg, 2012), allowing for one mismatch (score_min L, 0, −0.20). Methylation extraction was conducted with Bismark Extractor v0.22.3 (Krueger and Andrews, 2011) with the ‘ignore’ option to remove unmethylated cytosines introduced during the end-repair step. We combined data from the cytosines of both strands at each CpG site with a custom-made python script. The R-package Methylkit v.1.16.1 in R v.4.0 (Akalin et al., 2012; R Core Team, 2021) and its function methRead were used to load and analyse the methylation calls.

Methylation calls were filtered to those with a minimum coverage of 10 reads per CpG position. CpGs with a coverage >99.9th percentile most likely result from PCR bias and were removed. These two steps were performed using the filterByCoverage function in Methylkit. Across all 74 samples, there were on average 5,536,718 CpG positions before, and 661,447 CpG positions after, filtering for >10x and <99.9th percentile coverage (Supplementary Tables S1, S2). Methylation call distributions between samples were normalized using the normalizeCoverage function implemented in Methylkit to reduce the potential bias of systematic oversampling in some samples. We then merged CpG positions of the different samples using the unite function to make sure positions were sequenced in at least 70% of the samples (within and across the sexes) to facilitate a longitudinal analysis approach. Following Meng et al. (2010) and Sziráki et al. (2018), sites that showed little or no variation were removed as well, applying a threshold of a standard deviation <0.15. This resulted in 927,490 observations of 15,700 positions, distributed across all autosomal chromosomes (Supplementary Figure S1), which we used to assess autosomal methylation rates by using the number of positions read as our denominator and the number of cases in which the position was methylated as our dependent variable (see Statistical analyses).

2.5. Statistical analyses

To identify and partition sources of variation in autosomal methylation rates, we used the R package glmmTMB (Brooks et al., 2017) to run a generalized linear mixed model (GLMM) with the BFGS algorithm as an optimizer, using the number of methylated and unmethylated Cs (represented by the sequenced Cs and Ts) for each position as our dependent variable, assuming a binomial error distribution (Lea et al., 2017). As fixed effects we added ‘sex’ (as a two-level class variable using males as the reference category) and age (as a covariate). With respect to the latter, each individual’s age was partitioned into an ‘average age’ and ‘delta age’ component following van de Pol and Wright (2009). Average age was calculated as the average of all ages at which we assessed a bird’s autosomal methylation rate, while delta age was calculated as the difference between the bird’s actual age and its average age (i.e., delta age = age − average age). When adding both age variables as covariates, average age represents the among-individual, and delta age the within-individual age effect (van de Pol and Wright, 2009). If the among- and within-individual age effects were to differ, this would indicate that the effect of age among individuals cannot be explained by changes within individuals, thereby revealing age-specific selective (dis)appearance of individuals with certain levels of methylation (van de Pol and Wright, 2009). We additionally included the interaction between sex and delta age in our model to test for sex-differences in the within-individual age trajectory of methylation. Random effects included were ‘bird identity’ (1–34), genomic ‘position identity’ (1–15,700 scored as the base-pair-number on our scaffold) and an ‘observation’-level random effect (1–927,490 scored at the lowest level of observation, i.e., one value for each genomic position for each individual). The latter was added to (successfully) account for overdispersion detected by the R package performance (Lüdecke et al., 2021), whereas the first two were included to account for pseudoreplication caused by repeated sampling of individuals and genomic positions, respectively. Although we could have theoretically added breeding pair identity as a random effect, doing so would only add information when pairs regularly break up and reform among sampling events. Our common tern pairs all stayed together, such that we would not be able to disentangle the effects from pair and bird identify effects.

Because we found a significant interaction between sex and delta age (see Results), we also split our data set in one for males and one for females and ran sex-specific analyses including ‘average age’ and ‘delta’ age as fixed effects and ‘bird identity’, genomic ‘position identity’ and ‘observation’ as random effects. Moreover, by including positions that were sequenced in at least 70% of the samples, we facilitated a longitudinal analysis approach and a comparison between males and females, but excluded CpGs from the sex chromosome. To also assess age-specificity in methylation rate of the W chromosome, we also ran the sex-specific model for females for the 103 CpG positions we could assess there by including positions that were sequenced in at least 70% of the female samples.

Model evaluation and summary of parameter estimates and statistics were conducted using the R package parameters (Lüdecke et al., 2020). Results were plotted using the R-package sjplot (Lüdecke, 2018).

3. Results

3.1. Reference genome

The reference genome was generated using 68x PacBio sequencing reads, Bionano Genomics optical maps, 10X Genomics linked-reads and Arima Hi-C reads. This allowed us to scaffold the assembly to chromosome-level (Rhie et al., 2021) and we successfully assigned 99.3% of the assembled sequence to 25 identified autosomes, two sex chromosomes and the mitochondrial genome, leaving only 95 scaffolds unlocalised. The total length of the primary haplotype assembly was 1.23 Gbp with a contig N50 of 22.0 Mb and a scaffold N50 of 85.5 Mb. The assembly included 96.0% complete single-copy and 0.4% duplicated orthologs according to the BUSCO analysis. Only 1.3% of the gene models were fragmented, and 2.3% missing (n = 8,338 genes). This represents a high-quality assembly, surpassing the aspired VGP contiguity metrics ~20-fold (Rhie et al., 2021).

3.2. Age-specific global methylation

Across the 74 samples, 29,044,662 RRBS reads (~ 55.6% mapping efficiency) could be uniquely mapped to the reference genome (Supplementary Tables S1, S2). The mean bisulfite conversion rate was 99.2% (Supplementary Tables S1, S2). 15,700 autosomal CpG positions occurred in at least 70% of samples within and across the sexes, and across these positions methylation rates ranged from 0 to 1 in both males and females (Supplementary Figure S2). The average predicted methylation levels per sample ranged from 0.437 to 0.752 (mean: 0.562) in females and from 0.424 to 0.777 (mean: 0.627) in males (Figure 1).

Figure 1. Within-individual changes in predicted average autosomal methylation rate for 17 male (black circles with solid lines) and 17 female (white circles with dotted lines) common terns sampled at 1, 3 and/or 4 year intervals. The x-axis, termed delta age, represents the within-individual difference in age (in years) in relation to the average age at sampling (at negative values the individual was younger than its average age at sampling, at positive values it was older). The predicted methylation rate on the y-axis is extracted from the model presented in Table 1, and based on 15,700 autosomal CpG positions.

When analyzing sources of variation in methylation rates across the 15,700 autosomal CpG positions, we found a significant main effect of sex, as well as a significant interaction between sex and delta age (Table 1). The main effect of sex showed males to have higher methylation rates than females, the interaction term indicated that methylation rates showed different age-specificity in males and females. Subsequent sex-specific analyses showed that methylation rates declined with age within individual females, but did not change with age within individual males (Figure 1; Supplementary Tables S3, S4). In the main model (as well as in the sex-specific models, Supplementary Tables S3, S4), the effect of average age was non-significant and the credible intervals of the average and delta age components strongly overlapped (Table 1), suggesting no selective (dis)appearance of individuals based on their autosomal methylation rates.

Table 1. Results from a generalised linear mixed model with a binomial error distribution testing whether variation in autosomal methylation level is explained by sex (males as a reference category) and the between- (average age) and within-individual (delta age) components of age.

When assessing the age-specificity of methylation rates in the 103 CpG positions of the female W chromosome, neither average nor delta age was significant (Supplementary Table S5).

When comparing the random effects in the main model, most variance was explained by position identity (Table 1), showing that genomic positions vary in their average rate of methylation. Bird identity, on the other hand, explained variation in methylation rates to a much lesser extent (Table 1), such that there is little evidence for consistent variation in methylation among individuals.

4. Discussion

DNA methylation patterns at CpG sites are increasingly used as biomarkers, so-called epigenetic clocks, to predict both chronological and biological age across species and taxa (e.g., Lu et al., 2021). How DNA methylation rate changes within individuals and whether it can explain phenotypic senescence patterns, however, is still largely unknown (Bell et al., 2019). Here, we used blood samples from common terns collected at 1-, 3- and/or 4-year intervals and a longitudinal analysis approach to investigate whether autosomal methylation rates changed with age within individual birds, and whether any change differed between males and females sharing environments and broods. Based on our selection of 15,700 positions, we found female genomes to be generally less methylated than those of males, and to also become even less methylated as these females aged, whereas there was no such age-specific decline of autosomal methylation rate in males. Moreover, we found the estimates for the within- and among-individual components of age to be similar, such that there was no indication for selective (dis)appearance of individuals based on their methylation rate. Finally, we provide evidence for positions in the genome to consistently vary in their methylation rates, whereas evidence for consistent variation among individuals was considerably less.

Our finding of female common terns showing a decrease in methylation rate as they aged fits with findings of global hypomethylation in older compared to younger mammals (Zampieri et al., 2015; Ciccarone et al., 2018; Unnikrishnan et al., 2019). Such methylation loss is thought to partly originate from demethylation of large regions of repetitive sequences, CpG-poor promoters or large hypomethylated blocks of "open sea" regions outside the CpG islands (Bollati et al., 2009; Heyn et al., 2012; Yuan et al., 2015). Whether these changes affect chromatin configuration and thus genome (in)stability, and whether changes in promoter methylation interact with histone modifications and transcription factors to alter expression remains to be investigated (Zampieri et al., 2015; Ciccarone et al., 2018).

Male terns, in contrast, showed no signs of decreased methylation as they grew older. Although many studies developing epigenetic clocks have assumed age-related changes to be similar across the sexes and used mixed-sex datasets to obtain them [e.g., Horvath et al., 2021; Raj et al., 2021), others have found sex-differences in these clocks (e.g., in some human ethnicities (Horvath et al., 2016), baboons (Anderson et al., 2021) or elephants (Prado et al., 2021)]. Moreover, a rare partly longitudinal study in a wild population of roe deer also showed sex-specific epigenetic clock regions, with an accelerated ageing signal in males, which are known to undergo stronger survival senescence in this species (Lemaître et al., 2021). Combined with our findings, this suggests that tests for sex-specificity of methylation should best be the norm.

Interestingly, male and female terns from our study population do not differ in the onset or rate of senescence in survival or breeding probabilities (Zhang et al., 2015b; Vedder et al., 2021b), nor in their average lifespan (7.4 years for both males and females, Bouwhuis et al., 2015), such that sex-specificity in the ageing process is only found in how parental age affects the quality of the offspring that recruit back into the population [with maternal age negatively affecting the reproductive performance of daughters and paternal age negatively affecting survival of sons (Bouwhuis et al., 2015)]. As such, we did not necessarily expect sex differences in the age-specificity of the birds’ autosomal methylation level. The fact that we were able to observe them, raises the question of which site-specific methylation patterns drive the pattern observed on the global level. As mentioned above, global loss is thought to originate from the demethylation of specific regions: repetitive sequences, CpG-poor promoters or large hypomethylated blocks of "open sea" regions outside the CpG islands (Bollati et al., 2009; Heyn et al., 2012; Yuan et al., 2015). CpG island promoters, on the other hand, have been found to show age-specific increases in methylation (Heyn et al., 2012; Day et al., 2013). As such, global demethylation may perhaps be compensated for by such increases in males, but not females. In combination with the fact that we found genomic positions to consistently differ in methylation levels, this stresses the need for moving from the level of global autosomal methylation assessment that we and others (e.g., Watson et al., 2019; Sheldon et al., 2020) have started with, to fully annotating the (common tern) genome and studying patterns across genomic features and focal sites.

Besides finding a sex-specific overall level and within-individual change in autosomal methylation level with age, and significant among-position consistency in methylation level, we found little evidence for consistent among-individual levels of autosomal methylation across years, or for selective (dis)appearance of birds in relation to their methylation level. This suggests that birds differentially change their autosomal methylation from year to year, i.e., show different gaps between epigenetic and chronological age each year, with these changes or gaps perhaps reflecting their individual-specific condition or environment, but not relating to their breeding status (breeding versus non-breeding) in the study population or their local survival. Implementing a random regression analytical framework is data-hungry and not possible with our current dataset, but linking among-year changes in e.g., body mass, other measures of physiology or age-specific reproductive performance seems a promising research avenue, especially when taking the analyses to a site-specific level, such that distinct genotype-methylation-phenotype correlations can be identified.

The strength of our study lies in its longitudinal sampling and analytical approach, which allows us to characterise within-individual changes, rather than infer them from cross-sectional data. At the same time, however, this sampling approach may come with some limitations. Non-destructive, longitudinal sampling in natural populations often relies on using blood as the focal tissue, but how DNA methylation of (in our case) erythrocytes translates to phenotypes mostly remains an open question that needs addressing (Husby, 2020), ideally in experimental study systems. Moreover, because we used RRBS, a high-throughput and low-cost method, to assess methylation, our findings may mostly pertain to methylation of high density CpG regions (Smith et al., 2009; Gu et al., 2011). As such, we may have focally answered the question of what happens at CpG islands during ageing (Beck et al., 2021), something which annotation of the common tern genome will be able to tell. Keeping this in mind, however, our study has provided evidence for a sex-specific overall level and within-individual change in autosomal methylation with age and has shown that CpG positions in the genome vary consistently with respect to their methylation levels, such that future work, changeing the perspective from genome-wide average estimates to the specific genome feature or base pair level, can elucidate whether, where and how much methylation might affect ageing males and females, as such establishing a longitudinal epigenetic clock.

Data availability statement

The dataset and bioinformatic code presented in this study can be found at Zenodo (10.5281/zenodo.7533890 and 10.5281/zenodo.7493357).

Ethics statement

The study was performed under licenses of the city of Wilhelmshaven and the Lower Saxony State Office for Consumer Protection and Food Safety, Germany.

Author contributions

SB collected the blood samples and life history data on the terns. SB and ML designed the study. ML organized sequencing of the samples. BM analyzed the methylation data and performed the statistical analyses with the help of MM. The reference genome was generated as part of the Vertebrate Genomes Project, coordinated by EJ. BH and JM were responsible for sample processing and generated the genome sequence data. CC, GF, and MU-S assembled the genome under coordination of OF. WC, JW, and KH curated the final reference genome. BM and SB wrote the manuscript with input from ML, MM, and EJ. All authors contributed to the manuscript and approved the submitted version.

Funding

This study was supported by the Max Planck Society through a MPRG grant (MFFALIMN0001 to ML), a sequencing grant (PSLIMN6000 to ML and SB), and the “Norddeutscher Wissenschaftspreis 2018” (to BM, SB and ML).

Acknowledgments

We are indebted to Peter H. Becker for setting up the long-term individual-based common tern study at the Banter See and to Götz Wagenknecht and the many scientists, students and field assistants who collected the data across the years. We express our gratitude to Johanna Klughammer, Amelie Nemc and Christoph Bock for their advice on sequencing strategy and for performing the RRBS sequencing. We also thank Michaela Schwarz for help with DNA extraction, Kristian Ullrich for his support in solving bioinformatics hurdles, Melanie J. Heckwolf, Robin Cristofari and our epigenetic journal club for inspiring discussions on DNA methylation and 2 reviewers for constructive comments that helped us improve our paper.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2023.982443/full#supplementary-material

Footnotes

References