Statistical analyses
All statistical analyses were performed with R 3.5.3 (R Core Team, 2014).
Variations with age
To test whether the leukocyte concentration and counts varied with age, we used the leukocyte concentration (log-transformed) as dependent variable and the age as an explanatory variable in a LMM and the counts of lymphocytes, neutrophils, monocytes and eosinophils as dependent variables and the age as an explanatory variable in four GLMMs with a Poisson distribution (appropriate for the observed distribution of count data). Body mass at capture, sex, capture date, year of capture and the interactions between capture date and year of capture and between age and sex, were further included as fixed effects. Because individuals were sampled several times over the years, we included individual’s identity as random intercepts (Table 1A). The functions “lmer” and “glmer” in the package “lme4” (Bates et al., 2015) were used to fit Linear Mixed Models (LMMs) and Generalized Linear Mixed Models (GLMMs) (Bolker et al., 2009). Models including linear, quadratic or no effect of age were considered and compared using Akaike’s Information Criterion (AICc). The age effect models with the lowest AICc (Table S2) were selected. Afterwards, exploratory variables were removed following a backwards elimination procedure. The results obtained with the backwards elimination procedure were confirmed using an AICc selection (Table S3). We measured zero-inflation and variance inflation factors (VIFs) in all our models using the R package “performance” (Lüdecke et al., 2020). The only correlations between fixed effects we observed were between the body mass at capture and the year or the capture date as expected from the huge annual variability in body mass and the increase in body mass due to fat accumulation during the active season. For all models, we checked a posteriori the distribution of the residuals to assess the fit of the models to the observed data. Since we observed moderate overdispersion (all dispersion ratios < 2.58) in some of our models (models for lymphocytes and neutrophils), we estimated all models’ parameters using a Bayesian approach. From the final models, we used the “sim” function from the R-package “arm” to simulate values from the posterior distributions of the model parameters (Gelman & Yu-Sung, 2020). The 95% credible intervals (CI) around the mean were obtained after 5000 simulations. Assessment of statistical support was obtained from the posterior distribution of each parameter. We considered a fixed effect to be important if zero was not included within the 95% CI.
Partitioned age effect
To separate within- from between-individual variation with age, we tested for between- and within-individual age effect, using the same models as above but partitioning the age of each individual into ‘average age’ and ‘delta age’ (following van de Pol & Wright, 2009)) (Table 1B). ‘Average age’ corresponds to the average of all ages at which an individual was sampled, and ‘delta age’ to the difference between its age at sampling and its ‘average age’. The ‘average age’ represents the between-individual age effect, which corrects for the potential selective disappearance of individuals, while ‘delta age’ represents the within-individual age effect (van de Pol & Wright, 2009). Models were selected using both a backwards elimination procedure and an information theoretic approach (see results in Table S4) and statistical support for parameters were estimated as above.
Finally, to test if the between- and within-individual age effects were significantly different, which would indicate selective (dis)appearance, we ran the five selected final models including both age and ‘average age’ as explanatory variables (Table 1C). In these models, ‘age’ represents the within-individual effect and ‘average age’ the difference between within- and between-individual effects (van de Pol & Wright, 2009).
Immune phenotype and survival probability
We tested whether the death probability depended on leukocyte characteristics with mixed-effects Cox right-censored regression models (Nenko et al., 2018; Ripatti & Palmgren, 2000; Therneau et al., 2003). These models included leukocyte concentration or counts as time-dependent covariates and survival as response variable using the “coxme” function in the “coxme” R package (Therneau, 2018). The age at first capture and the sex were also included as fixed effects. Individual identity and year of birth were added as random effects to take into account repeated measurements and cohort effects (Table 2). The data were encoded with a zero as starting point for all individuals and with the years to death, to the end of the study, or to the next capture (for individuals with repeated data) as stop (Therneau, 2018). For the repeated data, the next interval started with the end of the previous interval. A ‘1’ was assigned to the event variable, if the individual died during the interval. We assumed that an individual died if it was neither captured nor observed the following spring (monitored until 2018). A hazard ratio higher than one indicates that the corresponding explanatory variable is associated with an increased risk to die. All individuals were followed until death (n = 27 for leukocyte concentration and n = 43 for leukocyte counts) or still alive in 2018 (n = 4 for leukocyte concentration and n = 6 for leukocyte counts). Three individuals were excluded from this analysis because their fate (alive or dead) was uncertain, due to capture permit forbidding to monitor their families in 2017 and 2018.