Statistical analysis
We used t-tests to compare the two nest box temperature measurements Ta and Tbr of our control and heated nests during the six hours of daily treatment (see above). For each nest we compared growth rates (K ) for nestling mass using nonlinear mixed models (package nlme Pinheiro et al. 2017) following the methodology of Sofaer et al. 2013. This approach uses age and treatment as fixed effects, and nest and nestling identity as random effects to account for the lack of independence among siblings of the same brood and for repeated measures of the same individual over time.
We initially tested the effects of simulated heat waves by fitting linear mixed models using the packages lme4 and lmerTest (Bates et al. 2014; Kuznetsova et al. 2017) with telomere length as the dependent variable, body mass, growth rate, as covariates, and brood size, hatching order, treatment, sex and age plus their interactions as categorical fixed factors. Due to the strong correlation between offspring mass and age during post-natal growth we created a variable called mass index that encapsulated the standardized differences in mass according to the age class of each subject, with a high mass index meaning that individuals are relatively heavy for their age and vice versa (see more details in S3.1). We included random slopes for individuals to account for repeated measures of subjects belonging to the same brood at different times. From this initial model we dropped the non-significant variables and fitted our best model that included only treatment, age, mass index and the interaction between treatment and age, as well as individual identity as random factor (see Section S3.2). To ensure that we were not over-fitting the data by including the interaction between treatment and age, we also performed a model comparison test using Akaike Information Criterion (AIC) (Akaike 1987) and Bayesian Information Criterion (BIC) (Schwarz 1978). This test indicated that our model was preferred by both approaches (AICweights = 99%; BICweights = 80%) to a model with an additive effect of treatment and age and to another with no interaction at all (see Tables S3.2.3). We conducted post-hoc pairwise comparisons between groups using the emmeans package (Lenth et al. 2019), which applies Tukey adjustment and Kenward-Roger degree of freedom adjustments. We also built a General Linear Model with time spent brooding as a dependent variable and treatment as a fixed factor to test the effects of heat waves on time spent brooding by the parents. Since this latter model indicated an effect of treatment on brooding effort, we attempted a test of the effects that variation in brooding rate has on telomere length and attrition by including in our initial model brooding and its interaction with age and treatment as a covariate. This part of the analysis however needs to be interpreted with extreme caution (see sections S6). Indeed, AIC suggested the best model to be the one including the three-way interaction while BIC gave opposite indications; also results were highly sensitive to influential observations (S6.3.3). All models assumptions and performances (homogeneity of variance, normality of residuals, normality of random effects, collinearity, and influential observations) were checked (see Supplemental Materials). All analyses were performed with R version 4.1.2 for Mac (R Core Team 2021).