From individual to population level: Temperature and snow cover modulate fledging success through breeding phenology in greylag geese (Anser anser)

Local weather conditions may be used as environmental cues by animals to optimize their breeding behaviour, and could be affected by climate change. We measured associations between climate, breeding phenology, and reproductive output in greylag geese (Anser anser) across 29 years (1990–2018). The birds are individually marked, which allows accurate long-term monitoring of life-history parameters for all pairs within the flock. We had three aims: (1) identify climate patterns at a local scale in Upper Austria, (2) measure the association between climate and greylag goose breeding phenology, and (3) measure the relationship between climate and both clutch size and fledging success. Ambient temperature increased 2 °C across the 29-years study period, and higher winter temperature was associated with earlier onset of egg-laying. Using the hatch-fledge ratio, average annual temperature was the strongest predictor for the proportion of fledged goslings per season. There is evidence for an optimum time window for egg-laying (the earliest and latest eggs laid had the lowest fledging success). These findings broaden our understanding of environmental effects and population-level shifts which could be associated with increased ambient temperature and can thus inform future research about the ecological consequences of climate changes and reproductive output in avian systems.

Onset of breeding and length of egg-laying time window. On average, geese started laying eggs on the 76th day of the year ± 11.0 days (ordinal date ± SD). The earliest egg-laying was observed in 2002 at the 51st day of the year, and the latest start was on the 116th day of the year in 2005 (Fig. 2). The first egg laid in the flock per season was associated with winter temperatures (Fig. 3a) and with annual snow cover ( Fig. 3b) with earlier onset of breeding after milder winters (higher winter temperatures and less snow (Table 1(I), supplementary Figure S2). The length of the egg-laying time window differed across study years (min. = 18 days in 1994; max. = 62 in 2002, average days = 38.38 ± 10.84 SD) and was positively associated with the number of females that attempted to breed within the flock (Table 1(II), Fig. 3c). However, there was no relationship with ambient temperature or snow cover. Breeding phenology on a pair level (i.e., the date of the first egg laid) was predicted by average winter temperatures and age, whereby egg-laying was earlier when winter temperatures were warmer (Fig. 4a, estimate − 0.33 ± 0.06, − 0.452 LCI − 0.202 UCI) and when the breeding female was older (Fig. 4b,   . Cubic regression spline smoothers with 95% confidence intervals were added to aid visual interpretation. The smoothers explain 35.8% and 28.1%, respectively, of the deviance. Productivity. In our population, a pair incubated up to 14 eggs per breeding season and had a maximum of 9 fledglings (mean ± SD: 5.7 eggs ± 2.1; 1.2 fledged ± 1.8). Fewer eggs were produced earlier in the study period and more eggs were produced later in the study period, thus, annual variation over time was the strongest predictor for average clutch sizes (Fig. 5a, estimate 0.72 ± 0.25 SE, 0.208 LCI 1.232 UCI, Table 3(I), supplementary Figure S3). Within the whole flock (12-37 successfully breeding pairs per year), no fledglings were produced in 1997, 1999 and 2002 while a maximum number of 43 fledged goslings reached in 2015 as well as in 2018. Average annual temperature was the strongest predictor for the proportion of fledged goslings per season (Fig. 5b, estimate 0.47 ± 0.28 SE, 0.029 LCI 0.901 UCI, Table 3(II)). Similar patterns were seen following milder winters with less snow cover, while the timing of breeding could not predict fledging success (neither the first egg in the flock nor the length of the egg-laying time window, supplementary Figure S4). On a pair level, earlier clutches had more eggs (Fig. 6a, estimate − 0.26 ± 0.05, − 0.357 LCI − 0.156 UCI, Table 4(I)) and there was some indication that larger clutches were produced following a winter with more snow cover (Fig. 6b, estimate 0.12 ± 0.04, 0.037 LCI 0.202 UCI). Clutch sizes did not differ with female age. The proportion of fledged goslings (ratio between the number of fledged goslings and the number of hatched eggs) was predicted by the timing of breeding, in a way that more hatchlings fledged earlier in the season and fewer hatchlings fledged with a later onset of breeding (Fig. 7a,  www.nature.com/scientificreports/ UCI), whereby the highest proportion of goslings fledged approx. between ordinal day 60 and 80. Furthermore, female age predicted the proportion of fledged goslings (Fig. 7b), whereby the within subject effect showed that females raised a higher proportion of goslings when they were younger than when they were older (estimate − 0.28 ± 0.14, − 0.559 LCI − 0.001 UCI), while the between subject effects indicated that older females, overall, had higher hatch/fledge ratios (i.e., raise a higher proportion of their produced offspring successfully until fledgling; estimate 0.33 ± 0.18, − 0.018 LCI 0.685 UCI).   www.nature.com/scientificreports/

Discussion
The average annual temperature increased locally by 2 °C across the 29-years study period, which is at the lower end of the values known for warming of the alpine areas 60 . Still, our longitudinal data show an advanced onset of the breeding season (both at individual level and also as the date of the first laid egg within the flock) following warmer winters. Additionally, earlier breeding predicted larger clutch size and higher gosling survival (i.e., hatch/fledge ratio). Phenological responses to ambient temperature can have fitness consequences. For example, earlier hatched young may have higher fitness if environmental quality deteriorates across the season, thus making timing per se the main driver of fitness loss over time ('date hypothesis' 44,46 ). However, performance-related responses to increasing temperatures tend to be non-linear and often follow quadratic relationships with a pronounced optimum 61 -a pattern we found in this study, as gosling survival was lower for females that lay too early or too late in the season. Phenological mismatches between the species and the environment can have negative fitness consequences. Advancing laying date according to climate warming may incur a trade-off between the most advantageous conditions for the parent versus the most suitable conditions for the offspring 31,62 . Early egg-laying dates could favour a competitive advantage among offspring feeding on the first green shoots and grasses low in tannins or other chemical plant defences 31 . However, postponing egg-laying might allow for clutches with more Table 2. Linear mixed effects model for the timing of breeding (ordinal date of egg laying) in relation to female breeding age and weather predictors of greylag geese in the Alm valley, Austria, (n = 619 breeding records of 155 females over 29 years). Best model (AIC c = 1411.3, ωi = 0.467; conditional R 2 = 0.703, marginal R 2 = 0.156) with random intercept, random slope (1 + age|female ID) and their correlation (r = 0.43) to evaluate plasticity. The most parsimonious model featured the mean centred value of female age (within subject effect), the mean female age (between subject effect), and winter temperature whereby older female laid eggs earlier and warmer winters were associated with earlier egg laying. Confidence intervals of parameter estimates not including zero (i.e., considered significant) are displayed in bold. Note that all quantitative input variables were scaled and centred. Null model ranked last [Akaike's information criterion corrected for small sample size (AICc) = 1435.3, ΔAICc = 23.99]; model weight based on complete list of candidate models (n = 40).   (Table 3), with 95% CIs in shaded grey. Dots represent the raw data. www.nature.com/scientificreports/ eggs because females have more time to consume resources to invest into egg production 63,64 . Later egg-laying may also offer advantageous foraging conditions during the incubation period 26,45,48,[65][66][67] . Individual quality has been shown to play a role in the timing of egg laying 31 . Parents with early onset of breeding tend to be of higher individual quality 50 , for example through being older and having more breeding experience 49,68 , better body condition 45,69,70 or occupation of better feeding areas; all of which could yield fitness benefits ('individual quality hypothesis' 45,[48][49][50]68,69 ). Such hypotheses are not necessarily mutually exclusive, but parents may have to face a trade-off considering breeding benefits (which might be related to the date hypothesis) as well as fitness costs (which might be related to the quality hypothesis) associated with the timing of breeding 42,46 . While climate change is known to drive range shifts and phenological shifts with population consequences across different species and taxa [71][72][73][74][75][76] , there is scant evidence of the potential costs of early versus late egg laying in relation to changes in ambient temperature. These potential effects of being too early may become more pronounced and measurable as climate change impacts increase.
Both, age and breeding experience have been shown to influence breeding productivity, with many studies showing age-related fitness parameters increasing over time [77][78][79] . This increase is likely associated with physiological maturation and increased experience with improved parental care and greater foraging efficiency both to reach breeding condition and raise offspring 80,81 . The potential role of the timing of breeding on age-related productivity is less understood and may or may not contribute to the mechanisms explaining this pattern. If timing of breeding is dependent on individual quality 50 , earlier laying might be predicted for individuals of higher quality (i.e., increased breeding experience with age), or for those better able to cope with physiological stress 82 Thus, if individual quality increases with age, one might predict that older birds breed earlier, which is also seen in our study. In addition to the effects of climate, the complex sociality and the long-term monogamy of the greylag goose needs to be considered, as the behavioural fine tuning of female and male partners are essential for optimizing offspring survival 83 .  www.nature.com/scientificreports/ Senescence is defined as a decrease in physiological function that leads to age-specific decreases in fitness components 84,85 . Our data hint at reproductive senescence because reproductive output as measured by offspring survival decreased with age measured at the individual level in females. At the between individual level, however, we found the opposite relationship and older females in the population had higher fledging success relative to younger females. This suggests differences in female quality, perhaps mediated by experience and/or a competitive advantage for access to resources, as a predictor of fledging success. Future research should aim to disentangle possible processes that could underpin these individual versus group-level patterns including selective disappearance 86 . It is possible that the progressive appearance of better-quality individuals with age 49 explains the findings.
Concluding, our longitudinal study addresses organism-environment interactions that influence fitness and species response to climate change in a model system, the greylag goose. Greylag geese can be considered among the more flexible goose species as they breed across a broad geographical area, from the Mediterranean nearly into the Arctic; therefore, they are suspected to be capable of coping with diverse conditions and could profit from warming 58 . Still, there are potential study limitations that need to be addressed. For instance, in addition to natural food availability, also supplemental food is predicted to advance gonadal development and lay date 87 . The study animals have been supplemented with food twice per day year-round since 1973, and therefore we are confident that the effects on gonad development from food supplementation have remained unchanged for almost 50 years. Furthermore, a comparison with other populations at similar and/or different latitudes might help disentangling the role of the environment. Furthermore, resident and migrant birds seem to respond differently to climate changes in their breeding areas [88][89][90] whereby arrival date at the breeding grounds may constrain a flexible laying date in migrants while such a constraint may be absent in resident birds. Therefore, plasticity seems to be an important mechanism to adjust for climate warming 91 , even though phenotypic and genetic mechanisms might both be involved.
Our study adds two important perspectives to this growing field of climate change research: (1) egg laying was earlier in years with higher ambient temperature and older females had earlier egg-laying date; and (2) an earlier onset of egg laying was related to both an increased clutch size and proportion of fledged goslings. Thus, a local increase in ambient temperature was associated with a local increase in productivity through phenological shifts. As this study uses long-term data from individually marked geese, it is possible to disentangle individual-level and flock-level patterns of response in relation to climate, which sharpen focus on casual mechanisms associated with fitness beyond phenology.  www.nature.com/scientificreports/ Collection of breeding data. Data were collected over 29 years, from 1990 to 2018, during which period the flock consisted of on average 149.14 individuals (SD = 15.64; Table 5). The flock was at its minimum in 1998 (147 birds) whereas the maximum number of individuals was registered in 2016 (180 geese). As the number of individuals in the flock may undergo seasonal variation, we set Jan. 15th as the point for the calculation of the total number of individuals in the flock per year. This was inferred according to the information collected in the frame of our long-term monitoring, which includes a regular check (two to three times per week) of the individuals present at the feeding site. Similarly, we set Feb. 15th for the estimation of the number of sexual mature females that could potentially have a nest in the forthcoming breeding season ( Table 5). All data (egg-laying, clutch size, hatching success, fledging success) were recorded 2-3 times per week during the reproductive season of the geese (approx. mid Feb. to end of July). During the mating period, the breeding huts and the traditional nesting locations were checked every two days, eggs weighed, measured and numbered according to their laying sequence. Later on, during the rearing phase, individual families were monitored every two days and the number of the accompanying offspring (i.e., goslings) was documented. The start of egg-laying was calculated as follows: (a) at pair level, the ordinal date of the first laid egg was used; (b) at flock level, the ordinal date of the first egg laid by the first breeding pair was used. The first egg was laid between the Model-averaged coefficients from a set of 3 models for clutch size; and, from a set of 4 models ratio of fledged goslings with ΔAICc < 2.0 presented as estimated values ± (unconditional) SE, lower and upper 95% CIs; confidence intervals of parameter estimates not including zero (i.e., considered significant) are displayed in bold. Note that all quantitative input variables were scaled and centred; (I) Null model ranked 59th [Akaike's information criterion corrected for small sample size (AICc) = 1528.6, ΔAICc = 18.07]; Note that in the clutch size model selection, we did not include 'year' as random intercept (throughout) because the variable caused singularity and explained 0% of the variance, which means that the fitted linear mixed model is very close to being a linear model (as indicated by comparing the lmer() model output with the lm(), results are very similar). In such cases the random term should be removed 119 (Table 4), with 95% CIs in shaded grey. Background scatter represents raw data (n = 380 breeding records with known fledging success of 96 individual females over 29 years). www.nature.com/scientificreports/ 51th and the 75th day of the year, whereas the last egg was laid between the 84th and the 116th day. The close of the egg-laying time window was determined by the ordinal date of the first egg laid by the last breeding pair. Only breeding pairs with known identity were used (individually marked males and females). In total, 614 laying dates of 300 individuals (148 males, 155 females) were used, with an average of 21.17 ± 5.67 (SD) breeding pairs per year. The possibility for unpaired females to successfully breed is low. Specifically, over the 29 years of data collection, a total of 22 nests were maintained by 20 solo females without a partner. Of these 22 nests, two produced fledged goslings (in sum four, one gosling in 1992 and three goslings in 1993). Importantly, clutches from solo females are so-called 'collection clutches' that include dumped eggs by other females, which adds a level of uncertainty about which individual began egg-laying and how many females contributed how many eggs; therefore, they were excluded from all analyses. We also excluded replacement clutches from the analysis (a total of 26 known replacement clutches over 29 years of which five produced a total of nine fledged young: specifically, three fledglings from two different females in 2012, three fledglings from one female in 2013 and three fledglings from two different females in 2017). The difference between the number of males and females included in the study might be related to mate switches, which can occur when one partner has died 51 . In our data approx. 10% of females and 11% of males switched partner. The rarity of mate changes has implications for the random factor structure in the statistical analyses.

Methods
In addition to the timing of egg-laying, we recorded (c) clutch size, and (d) 'gosling survival' (i.e., the ratio between the number of fledged goslings and the number of hatched eggs). We defined a gosling as 'fledged' when it can fly autonomously. The variable gosling survival was used as a proxy for breeding investment versus breeding performance and is informative about the environmental conditions and the quality of parental care during the raising period. Information about the percentage of fledged goslings per year, are reported in Table 5.
Climate data. Climate data were obtained from a meteorological station ' Almsee-Fischerau' which was

Climate change.
To visualise climate trends over the course of the study period, we fitted generalized additive models (R-package 'mgcv') 95,96 with annual snow cover [cm], average annual temperature [°C] and winter temperature respectively as response variables (each following a Gaussian distribution fitted with an identity link function), and year as sole predictor variable (1989-2018) using cubic regression splines (the function itself decides the ideal number of knots required, which resulted in a cubic regression spline that fitted our data best). The GAM models serve sheer illustration purposes of the weather trends in our study area, why we plotted the raw data overlaid with the spline smoother with 95% confidence intervals to aid visual interpretation (see Zuur et al. 97 for a similar approach).
Model selection procedure. All breeding phenology and productivity models (detailed below) were fitted with the base, lme4 98 , betareg 99 and car 100 packages, model predictions were extracted with effects 101 and were visualised with lattice 102 , ggplots2 103 and base plots. Residual distributions of the models were inspected visually to assess model fit and potential deviations from the assumptions of normality and homogeneity of residuals by evaluating the model criticism plots produced by the 'plot' function in the base package and the 'mcp.fnc' in the LMERConvenienceFunctions package 104 . Additionally, we used Variance Inflation Factors (VIF) implemented in the performance 105 package, and considered values of VIF < 5 as low collinearity and an indication that predictors can be fitted in the same model without having problems of collinearity 106,107 . We used an information theoretic approach using Akaike Information Criteria (AIC) to derive the best, most parsimonious model if the next ranked model had a difference of ΔAIC c > 2 to the first ranked model, otherwise we used model averaging with the package MuMIn 108 and AICcmodavg 109 (see Grueber et al. 110 for details on multimodel inference) across all models with ΔAIC c < 2 and visualised a summary of model averaged effect sizes with sjPlots 111 . All quantitative variables were scaled (standardized to mean = 0 and standard deviation = 1) to bring the variables to comparable dimensions and to facilitate the correct interpretation of effect sizes, but were back transformed for visualisation purposes. Our data set contained no missing values, ensuring accurate model comparisons throughout the selection and, if applicable, averaging process. For weather predictors, we initially explored the linear and quadratic relationships-as temperature effects are usually non-linear 61 -in the model selections. In order to do so, we used dependency chain arguments to ensure that quadratic terms were not fitted without considering the linear relationship, but these quadratic effects never featured into any top model sets with ΔAIC c < 2.0, why we removed the quadratic predictors to simplify the process and to reduce the model candidate list of potential predictor combinations. The model set was then ranked using ΔAIC c values. Akaike weights (ω i ) were calculated to assess the relative likelihood for each model considered 112 . Thus, ω i reflects the model's probability given the full model list rather than only those below a given threshold of ΔAIC c ; added R 2 values (pseudo R 2 for beta regression models and marginal and conditional R 2 for mixed effect models) were extracted with the performance 105 package. A table of best candidate models (up to ΔAIC c < 2.0) is presented in the results section while the complete candidate lists are available in supplementary Tables S1-S7. As mentioned above, models below this threshold were extracted and consequently used for model averaging 113 . We followed the guidelines www.nature.com/scientificreports/ provided by Burnham and Anderson 112 whereby ΔAIC c within 0-2 indicate models with substantial levels of empirical support. Model averaging is furthermore a useful method to ameliorate the effect of potentially uninformative parameters (i.e., when a variable with poor explanatory power is added to an otherwise good model and the result is a model with ΔAIC c < 2) 114 . We report the direction of parameter estimates and their magnitudes (effect sizes), and unconditional SEs and CIs (95% confidence limit) from model averaged coefficients. We report unconditional SE because this incorporates model selection uncertainty, as opposed to standard SE which only considers sampling variance 110 . We used confidence intervals to assess the magnitude of the effect and concluded that the estimate is different from zero (i.e., there is a 'significant' effect) when the confidence interval excludes zero.) Breeding phenology and productivity on flock level. To model phenology at flock level as a function of the covariates, two different response variables were analysed with a set of Linear Models (LMs) with identity link function. We used (I) the ordinal date of the first egg laid by the first pair; and (II) the length of the egglaying time window (measured as the timespan between the dates of the first egg laid by the first pair and the first egg laid by the last pair within the flock each year) as response variables. Fixed predictor variables in the saturated models were: annual snow depth (continuous), average annual temperature (continuous), average winter temperature (Dec-Feb; continuous), number of breeding pairs in the flock (discrete) and year (discrete). Their combinations resulted in a candidate list of n = 32 models (complete list in supplementary Tables S1 and S2).
To model breeding productivity at flock level as a function of the covariates, two different response variables were analysed. We used (I) average clutch sizes (i.e., a measure for breeding investment) in a set of LMs with identity link function; and, (II) proportion of fledged young (i.e., a measure for outcome in relation to investment) within the flock per year in a series of beta regression models (i.e., beta distribution in generalised linear model (GLMs)) best suited to fit proportion values within the binomial model group. Note that beta distributions are defined for all values between 0 and 1, but not 0 and 1 themselves, but we had zero goslings fledging in the years 1997, 1999 and 2002. To solve this issue, we followed the recommended lemon squeezer transformation from the vignette of the R package betareg 99 and narrowed the data 115 : Fixed predictor variables in the saturated models were: annual snow depth (continuous), average annual temperature (continuous), average winter temperature (Dec-Feb; continuous), number of breeding pairs in the flock (discrete), year (discrete), ordinal day of the first egg laid in a given year (discrete) and the length of the egg-laying time window (discrete). Their combinations resulted in a candidate list of n = 128 models (complete list in supplementary Tables S3 and S4).
Breeding phenology and productivity on pair level. To model phenology on pair level as a function of the covariates, the timing of breeding (ordinal day of egg-laying, hereafter 'egg-laying day') was analysed with Linear Mixed Models (LMMs) with identity link function. Fixed predictor variables in the saturated model were: annual snow depth (continuous), average annual temperature (continuous), average winter temperature (Dec-Feb; continuous), and age of the breeding female (discrete). We fitted age (in years) to test for a potential earlier onset of breeding with increasing female age (our proxy for 'breeding experience') as known from other avian systems 116,117 . To test if specifically, older (i.e., more experienced) females lay earlier following milder winters, we also included the interaction term between temperature and age ('average annual temperature × female age). Importantly, we were aiming to tease apart the within and between individual variation of age, why we included the mean centred value of age (within subject effect) as well as the mean female age (between subject effect) 47,86,118 . The predictor combinations resulted in a candidate list of n = 40 models (complete list in supplementary  Table S5). Throughout, 'year' and 'female ID' were included as random terms to account for non-independence stemming from data coming from the same breeding females over several years and multiple measures from different females within the same year. We fitted the random intercept, random slope and the correlation between them in the models (structure in lme4: (1 + age|female ID)) to evaluate plasticity. We did not include partner ID because, as mentioned above, each subject is in 90% of the cases tested with the same partner which would create a redundant random effect that also resulted in severe convergence issues in the models. Furthermore, we used female ID rather than pair ID because of our interest in the age effect in shaping phenology. Again, because of the long-term pair bonds in greylag geese, including partner age as well would have created problems with collinearity in fixed effects and unidentifiable random slopes.
To model clutch sizes on pair level, a set of LMMs was fitted with identity link function. Fixed predictor variables in the saturated model were: annual snow depth (continuous), average annual temperature (continuous), average winter temperature (Dec-Feb; continuous), age of the breeding female (discrete; within and between individual effect), and egg laying day (ordinal date; considering the linear and quadratic relationship to account for known season effects and influences of breeding age on productivity in most avian systems 45,46 ). The predictor combinations resulted in a candidate list of n = 96 models (complete list in supplementary Table S6). We included the random intercept for 'year' , and the random intercept, random slope and the correlation between 'age' and 'female ID' in all models (1 + age|female ID).
To model the proportion of fledged goslings (i.e., egg-fledged ratio; using the column bind 'cbind' function with the number of eggs as binomial denominator) on pair level, a set of Generalized Linear Mixed Models (GLMMs) with binomial distribution and log-link function was fitted. Fixed predictor variables in the saturated model were: annual snow depth (continuous), average annual temperature (continuous), average winter www.nature.com/scientificreports/ temperature (Dec-Feb; continuous), age of the breeding female (discrete; within and between individual effect), and egg laying day (ordinal date; linear and quadratic relationship). The predictor combinations resulted in a candidate list of n = 96 models (complete list in supplementary Table S7). We included the random intercept for 'year' , and the random intercept, random slope and the correlation between 'age' and 'female ID' in the models (1 + age|female ID). Sample sizes differed between phenology (n = 619), clutch size (n = 559) and fledged goslings (n = 380) because not all detected breeding attempts (i.e., clutch initiation recorded) could be checked to confirm final clutch sizes, and not all families could be closely followed to determine the number of goslings.
Ethical statement. This study complies with all current Austrian laws and regulations concerning the work with wildlife. Observing the animals and controlling their nests were performed under Animal Experiment Licence Number 66.006/0026-WF/V/3b/2014 by the Austrian Federal Ministry for Science and Research (EU Standard, equivalent to the Animal Ethics Board). We confirm that the owner of the land, the Duke of Cumberland, gave permission to conduct the study on this site. All data were collected non-invasively. Birds were habituated to the presence of humans. The authors adhere to the 'Guidelines for the use of animals in research' as published in Animal Behaviour (1991,41,(183)(184)(185)(186).