Statistical analysis
We conducted all analyses in R 4.1.2 (R Core Team, 2021). To allow the use of RTL as a predictor, we first corrected RTL measurements for technical effects, including storage time, technician identity, and qPCR plate effects (Chik et al., 2023; Sibma, 2021). We ran a linear mixed model using the package lme4 1.1-28 (Bates et al., 2015), with RTL as the response variable, z-transformed such that effect sizes are comparable between studies (Verhulst, 2020). We fitted the following predictors: duration from when the sample was stored as a blood sample until DNA extraction (‘BloodAge’, in years), duration from when the DNA sample was stored until RTL measurement (‘DNAAge’, in years), the squared terms for both storage durations to account for non-linear effects, technician identity as a two-level fixed factor (A/B), and plate identity as a random variable. We fitted this model assuming a Gaussian error distribution. The residual RTL values (‘corrected z-RTL’) were then extracted for the survival models.
We tested for the correlation between post-fledging telomere length and short-term survival in two ways. First, we ran a generalized linear mixed model (GLMM) using lme4 . We fitted annual survival (whether an individual survived one more year after sampling as an adult as a binomial response, 0/1) with a logit link function, and corrected z-RTL as a continuous fixed predictor. As survival changes with age non-linearly, we also fitted age at sampling (in years) and its squared term as continuous fixed predictors, along with sex. Finally, we added bird ID and year of capture as random variables to correct for non-independence in observations from the same bird, and from yearly stochasticity. We checked the variance inflation factor (VIF) of fixed predictors, and concluded that there was minimal collinearity as all VIFs were < 5 (Montgomery et al., 2013). Second, we fitted a time-dependent Cox proportional hazard model (n = 1,211). In brief, at the time of each death event, the model compares covariate values of individuals who died, to all other individuals who were still alive and therefore at risk of dying, to estimate the risk score associated with the covariate value. To run this model, we coded time-to-event (death) for each individual in days, in a step-wise manner, with the first step being the time elapsed between the hatch date and first RTL measurement, and subsequent steps being the time elapsed between two consecutive RTL measurements, and the last step being the time elapsed between the last RTL measurement and the date the bird was last seen, on which it was assumed dead. We excluded birds without an exact hatch date. We then ran the Cox model using the package coxme 2.2-18.1 (Therneau, 2022), right-censoring birds that were still alive at the time of the analysis, and using the same fixed effect and random effect structures as the binomial GLMM.
We then tested whether adult telomere dynamics are associated with lifespan, using a bivariate model, which allows estimation of the covariance among the two response variables, and allows fixed effects to be fitted to only one of the two responses. We did so in a Bayesian framework, fitted with MCMCglmm 2.32 (Hadfield, 2010), following a similar approach by Heidinger et al. (2021). In this model, we only included the 1,214 birds that were dead at the time of analysis. We fitted z-RTL, and individual lifespan as response variables, assuming respectively a Gaussian and a Poisson distribution. For z-RTL only, we fitted age at sampling in years, centred around the population mean, so that the random individual intercept in z-RTL can be interpreted relative to the population mean. We also fitted BloodAge, DNAAge, their squared terms, and technician identity as fixed variables, to correct for the age and technician effects on RTL. For the random effect structure, we fitted a random slope function of RTL by age at the individual level, along with the year of capture and plate ID as random variables to RTL. We did not fit social parent identities as they were found to explain minimal variation in RTL (Chik et al., 2023). Because each individual had one lifespan value, a random effect of individual on lifespan could be estimated as a part of the residual variation in the model. The random-residual effects structure therefore allows the estimation of the among-individual variance and covariance among RTL, rate of RTL change, and lifespan:
\begin{equation} \begin{bmatrix}\sigma_{\text{RTL}}^{2}&\sigma_{RTL,\ RTL:Age}&\sigma_{RTL,\ \ Lifespan}\\ \sigma_{RTL,\ RTL:Age}&\sigma_{RTL:Age}^{2}&\sigma_{RTL:Age,\ Lifespan}\\ \sigma_{RTL,\ \ Lifespan}&\sigma_{RTL:Age,\ Lifespan}&\sigma_{\text{Lifespan}}^{2}\\ \end{bmatrix}_{\text{ID}}\nonumber \\ \end{equation}
To examine the correlations between telomere dynamics and reproductive success, we fitted two bivariate mixed models, with LRS and ARS respectively. The LRS model had the same framework as the lifespan model, and included only the 1,214 individuals that were dead at the time of analysis. For the ARS model, we used the whole dataset. We paired ARS with the z-RTL measurement taken in the same year for each bird. We retained the same fixed effects structure on RTL, and modelled variance and covariance between z-RTL and ARS explained by bird ID and capture year. We also estimated the residual covariance between z-RTL and ARS in this model to examine the within-individual covariation between the two variables.
For all three bivariate models, we used default priors for fixed effects, inverse-Wishart priors for random effects, and adjusted the number of iterations, burn-in, and thinning intervals for each model, such that convergence is reached based on the following criteria: visual inspection of posterior trace plots showed no distinguishable trend, autocorrelation <0.1, and the effective sample size >1000.