Immune phenotype and variation with age
Within-individual variation in immune phenotype with age
To test whether the leukocyte concentration and profiles varied with age of an individual, we used the leukocyte concentration (log-transformed) as dependent variable in a Linear Mixed Model (LMM) and the relative numbers of lymphocytes, neutrophils, monocytes and eosinophils as dependent variables in four Generalised Linear Mixed Models (GLMMs) with a Poisson distribution (appropriate for the observed distribution of count data). The age was partitioned into “average age” and “delta age” components; where average age represents the among-individual age effect and delta age the within-individual age effect (van de Pol & Wright, 2009). The average age was calculated as the average of all ages at which an individual’s leukocyte concentration or profile were measured. The delta age corresponds to the difference between the individual age at measurement and its average age (i.e. delta age = age - average age). To investigate the potential non-linear within-individual variations of the leukocyte concentration and profile with age, we followed the equation (3) recommended in Fay et al. (2022). This equation proposes to model a within-individual quadratic effect of age as age² - (average age)².
Body mass at capture, capture date, year of capture and the interaction between capture date and year of capture were further included as fixed effects. Because individuals were sampled several times over the years, we included individual’s identities as random intercepts. Our models also included the interaction between “average age” ands sex as well as the interactions between sex and both the linear and quadratic effect of “delta age” to test for sex-differences in the leukocyte variation with age.
Among-individual variation in immune phenotype with age
Models including average and delta age cannot properly investigate the occurrence of selective appearance and disappearance of peculiar individual phenotypes from the population (Fay et al., 2022; van de Pol & Verhulst, 2006). For this purpose, we replaced the “average age” and “delta age” (linear and quadratic) effects in the previously selected within-individual models by the ‘age at access to dominance’ , to assess selective appearance, and the ‘age at last observation’ , to assess selective disappearance (Fay et al., 2022; van de Pol & Verhulst, 2006).
The functions “lmer” and “glmer” in the package “lme4” (Bates et al., 2015) were used to fit the models (Bolker et al., 2009). Final models were selected using a backward elimination procedure. We measured zero-inflation and variance inflation factors (VIFs) in all our models using the R package “performance” (Lüdecke et al., 2020). 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 effective 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.
Immune phenotype and survival probability
We tested whether the risk to die 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 profiles 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 right-censor stop points (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 of death. 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 a capture permit forbidding the monitoring of their families in 2017 and 2018.