Statistical analyses
To identify and partition sources of variation in autosomal methylation levels, 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 & Wright (2009). Average age was calculated as the average of all ages at which we assessed a bird’s autosomal methylation level, 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 & 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, revealing age-specific selective (dis)appearance of individuals with certain levels of methylation (van de Pol & 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’, genomic position (‘position identity’) and an observation-level random effect (‘observation’). The latter was added to account for overdispersion, which was detected using the R package performance (Lüdeckeet al. 2021), whereas the first two were included to account for pseudoreplication caused by repeated sampling of individuals and genomic positions, respectively. Model evaluation and summary of parameter estimates and statistics were conducted using the R packageparameters (Lüdecke et al. 2020). The marginal effects of the interaction were plotted using the R-package sjplot test (Lüdecke 2018).

Results