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.