Introduction

Phenological shifts, i.e. altered timing of seasonal life cycle activities or events, are often induced by climate change1 and can influence animal reproduction and population viability2,3,4. These shifts might be caused by changing abiotic factors and generally entail modified biotic interactions5. Earlier spring and longer frost-free seasons can advance the flowering of plants and egg-laying in birds6, which might prolong resource availability5 with associated fitness advantages. Such shifts can cause temporal mismatches between trophic levels, such as between plants and pollinators or predators and prey, thereby negatively affecting fitness7,8,9,10.

Birds are useful models to showcase how global climate change is affecting living systems11,12,13. Range shifts in both breeding and wintering areas14,15 as well as changes in body mass or size16,17 are among the well documented ecological effects of climate change on avian biology. The timing of avian migration and patterns of reproduction reveal the most conspicuous effects of climate change. For example, long distance migrants have been shown to migrate earlier to their wintering areas, while medium distance migrants were recorded later than expected for autumn migration18,19. During the last decades, many bird species have shifted their breeding time to earlier egg-laying, which has been associated with increasing ambient temperature4. Overall, there is evidence at latitudinal and regional scales for impacts of climate on avian behaviour with acknowledgement that complex underlying interactions could occur in relation to landscape characteristics, animal distributions, and behavioural activities at different spatial scales20,21. Sedentary non-migrating local populations offer a unique perspective to measure effects of climate on breeding phenology while controlling for possible interaction effects from migration or occurrence in other geographical areas at different times of their lives22.

Phenotypic plasticity has been considered as one possible mechanism driving changes in bird phenology23,24,25. Individuals are predicted to adjust their breeding phenology to suit prevailing or anticipated environmental conditions that enhance their reproductive output and/or survival25. Associations between laying date and fitness are well documented4, and a pedigree might help to disentangle the causal relationship between the two variables26,27,28, as both fitness and laying date may change in response to environmental conditions29. However, species with a constrained breeding season (i.e. predictable environment) seem to have a limited ability beyond existing plasticity to respond to changing environmental conditions29, whereas birds from a temperate zone should be adapted to cope with greater environmental variability13,30. In fact, birds breeding in seasonal environments are challenged to time their reproduction to match future demands of their offspring with the quality of the post-hatch environment31. Optimal breeding time windows may select for the capacity to respond to environmental cues such as ambient temperature, and for reproductive plasticity to respond to inter-annual variation in such cues32.

Mild temperatures and food abundance are the main limiting factors for the reproductive activities of many songbirds breeding in temperate zones33. Studies have shown earlier egg-laying dates associated with warmer ambient temperature in many bird species34,35. However, there is considerable variation in these patterns within and among avian species36,37. For instance, in a long-term study on 168 bird species, Parmesan and Yohe11 showed that 47% had earlier egg-laying date, 8% had later egg-laying date, and 45% showed no significant change with ambient temperature. In the tree swallow (Tachycinta bicolor), egg-laying across North America has been advancing in association with increasing ambient temperature38,39, a pattern also found in flycatchers (Ficedula spp.) across Europe40. Similar results were obtained for common eider (Somateria mollissima) in Canada whereby warmer spring temperatures predicted both earlier mean annual laying dates and earlier suitable conditions for the ducklings’ post-natal growth; this lead to increased breeding success and offspring survival as well as successful recruitment into the breeding population, suggesting that more low-quality females profited from nesting in warmer years41.

In avian biology, phenology matters because earlier clutches generally have more eggs, higher hatching and fledging rates, and earlier offspring are more likely to recruit into the breeding population42,43,44. The fitness advantages of earlier broods can be attributed to the quality of the environment, and/or of the parents26. The ‘date hypothesis’44,45,46 predicts that timing of breeding per se affects all individuals in the same way, through a deterioration of the environment as the season advances (i.e., within-subject hypothesis). The ‘individual quality hypothesis’ reflects quality differences between breeders, irrespective of timing per se (i.e., between-subject hypothesis). In other words: does egg-laying link to high reproductive success or do females that lay early also have higher reproductive success47? Individual quality (e.g. better body condition, or previous breeding experience)48,49 can be linked to a sequential onset of breeding according to the individual quality of the parents50.

In the present longitudinal study, we consider the effects of climate parameters on the breeding phenology and fledging success of a free-flying population of greylag geese (Anser anser) over 29 years (1990–2018). The birds are individually marked, which allows accurate long-term monitoring51. The free-flying flock does not migrate and receives food supplementation twice per day. The geese breed in open nests or breeding facilities on lakes and ponds at various locations in the study area (i.e., Alm valley, Austria). These free-flying geese are exposed to predation and to ambient temperature and snowfall, which has an effect on physiological parameters52 and could affect the timing of egg-laying and the duration of the egg-laying time window53. Both Arctic and generally temperate breeding geese have been shown to adjust the timing of nesting to the availability of food resources in order to successfully raise their goslings54,55. Many goose species are specialized Arctic breeders and their goslings might also feed on insects56; however, goose species are generally considered a keystone herbivore57 and greylag goslings exclusively forage on sprouting grass and herbs58, which have the highest protein content when they start to grow. In the present study, we investigate the relationship between ambient temperature and snow cover on breeding phenology (date of the first egg laid, duration of the egg-laying time window) and fledging success of individually marked geese nesting at the study site. The Alm valley experiences diverse climatic conditions and seasons across the year, as is typical for the central European mountainous region, with rather long and cold winters and short and warm summers59, which we predict will exert a significant influence on the breeding phenology and fledging success of greylag geese. The present study aims to: (1) identify climate patterns (mean winter and annual temperature, annual snow cover) at a local scale in the Alm valley in Upper Austria across 29 years, (2) measure the association between ambient temperature and snow cover and breeding phenology, and (3) measure the relationship between ambient temperature and snow cover on clutch size and fledging success in greylag geese. We expect (i) the average annual temperature will have increased locally, consistently with the global temperature patterns, and snow cover to have decreased. Following the date hypothesis, we expect that (ii) geese will have an earlier onset of egg-laying and a longer egg-laying time window associated with mild winters, which may allow more females in the flock to breed41; that (iii) earlier egg-laying will be associated with larger clutch sizes; and, that (iv) goslings will have higher fledging success under milder conditions. Following the individual quality hypothesis, we further expect that older (i.e., more experienced) females will be the ones to lay earlier and to have higher reproductive success, specifically following milder winters.

Results

Climate across 29 years

Over the course of the 29-years study period, average annual temperature increased (parametric coefficients: estimate 7.31 ± 0.13, P < 0.001; smooth term: F(24,28) = 4.03, P = 0.017; estimated smoothing curves of cubic regression terms in Fig. 1a). The mean ambient temperature increased by 2 °C from an average of 6.9 °C in 1990 to 8.9 °C in 2018. Similar trends were found when considering only winter temperature (Dec–Feb; supplementary Figure S1). Snow cover was highly variable, but did not show a significant trend over time (parametric coefficients: estimate 7.42 ± 0.93, P < 0.001; smooth term: F(25,29) = 2.72, P = 0.062; R2 = 0.217; Fig. 1b).

Figure 1
figure 1

(a) Average annual temperature [°C]; and (b) annual snow depth [cm] over the course of the study period (1990–2018). 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.

Full size image

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.

Figure 2
figure 2

Egg-laying time period (1st quartile, median, 3rd quartile, minimum and maximum) between 1990 to 2018 in relation to average winter temperature (Dec–Feb) in greylag geese in the Alm Valley, Austria.

Full size image
Figure 3
figure 3

Relationship between the first egg laid within the flock and (a) average winter temperature (estimate − 0.43 ± 0.25 SE, − 0.788 LCI − 0.075 UCI); and, (b) annual snow depth (estimate 0.38 ± 0.18 SE, 0.010 LCI 0.742 UCI). (c) Relationship between the length of the egg-laying time window and the number of breeding pairs of the flock (estimate 0.72 ± 0.33 SE, 0.048 LCI 1.384 UCI) in greylag geese in the Alm valley, Austria, between 1990 and 2018. Lines are the predicted relationship based on model results (Table 1) with 95% CIs in shaded grey. Background scatter represents raw data.

Full size image
Table 1 Linear effects model for (I) the first egg laid within the flock; and, (II) the length of the egg-laying 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 greylag geese flock) over 29 years (1990–2018) in the Alm valley, Austria, in relation to weather and the number of pairs attempting to breed within the flock.
Full size table

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, within individual effect: estimate − 0.13 ± 0.04, − 0.217 LCI − 0.054 UCI; between individual effect: estimate − 0.19 ± 0.06, − 0.309 LCI − 0.064 UCI; Table 2). Neither average annual temperatures nor snow cover explained much variance in the data.

Figure 4
figure 4

Relationship between breeding phenology and (a) winter temperature and (b) female age. Higher winter temperatures (Dec–Feb) are associated with an earlier onset of breeding (estimate − 0.33 ± 0.06, − 0.452 LCI − 0.202 UCI) and older females lay earlier than younger females (within individual effect: estimate − 0.13 ± 0.04, − 0.217 LCI − 0.054 UCI; between individual effect: estimate − 0.19 ± 0.06, − 0.309 LCI − 0.064 UCI) in greylag geese in the Alm valley, Austria, between 1990 and 2018. Lines are the predicted relationships based on the output of the LMM (Table 2), with 95% CIs in shaded grey. Background scatter represents raw data (n = 614 breeding records of 155 individual females over 29 years).

Full size image
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).
Full size table

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).

Figure 5
figure 5

Relationship between the (a) average clutch size in the flock and annual variation over time (estimate 0.72 ± 0.25 SE, 0.208 LCI 1.232 UCI); and, (b) proportion of fledged goslings in the flock and average annual temperature (estimate 0.47 ± 0.28 SE, 0.029 LCI 0.901 UCI) in greylag geese in the Alm valley, Austria, between 1990 and 2018. The line is the predicted relationship based on the output of the linear model for average clutch size and beta regression model for proportion fledged goslings (Table 3), with 95% CIs in shaded grey. Dots represent the raw data.

Full size image
Table 3 Linear effects model for (I) average clutch size in the flock; and (II) beta regression model for the proportion of fledged goslings per season over 29 years (1990–2018) in the Alm valley, Austria, in relation to the timing of breeding, the number of pairs attempting to breed, weather and annual variation.
Full size table

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, linear relationship: estimate 2.31 ± 1.15, 0.052 LCI 4.562 UCI, Table 4(II)). There was also a quadratic relationship apparent (estimate − 2.77 ± 1.16, − 5.055 LCI − 0.482 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).

Figure 6
figure 6

Relationship between clutch size and (a) the timing of breeding and (b) annual snow depth. Earlier clutches had more eggs (estimate − 0.26 ± 0.05, − 0.357 LCI − 0.156 UCI) and higher snow cover was associated with larger clutch sizes (estimate 0.12 ± 0.04, 0.037 LCI 0.202 UCI) in greylag geese in the Alm valley, Austria, between 1990 and 2018. Lines are the predicted relationships based on the output of the LMM (Table 4), with 95% CIs in shaded grey. Background scatter represents raw data (n = 559 breeding records with known clutch sizes of 145 individual females over 29 years).

Full size image
Table 4 (I) Linear mixed effects model for clutch size; and (II) generalized linear mixed effects model for the ratio of fledged goslings (binomial distribution with log-link function) over 29 years (1990–2018) in the Alm valley in relation to the timing of breeding, female breeding age, and weather.
Full size table
Figure 7
figure 7

Relationship between proportion of fledged goslings (ratio of fledged eggs) and (a) the timing of breeding (linear relationship: estimate 2.31 ± 1.15, 0.052 LCI 4.562 UCI; quadratic relationship: − 2.77 ± 1.16, − 5.055 LCI − 0.482 UCI) and (b) female age (within individual effect: estimate − 0.28 ± 0.14, − 0.559 LCI − 0.001 UCI; between individual effect: estimate 0.33 ± 0.18, − 0.018 LCI 0.685 UCI) in greylag geese in the Alm valley, Austria, between 1990 and 2018. Lines are the predicted relationship based on the output of the GLMM (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).

Full size image

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 areas60. 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 optimum61—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 offspring31,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 defences31. However, postponing egg-laying might allow for clutches with more eggs because females have more time to consume resources to invest into egg production63,64. Later egg-laying may also offer advantageous foraging conditions during the incubation period26,45,48,65,66,67.

Individual quality has been shown to play a role in the timing of egg laying31. Parents with early onset of breeding tend to be of higher individual quality50, for example through being older and having more breeding experience49,68, better body condition45,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 breeding42,46. While climate change is known to drive range shifts and phenological shifts with population consequences across different species and taxa71,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 time77,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 offspring80,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 quality50, 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 stress82 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 survival83.

Senescence is defined as a decrease in physiological function that leads to age-specific decreases in fitness components84,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 disappearance86. It is possible that the progressive appearance of better-quality individuals with age49 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 warming58. 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 date87. 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 areas88,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 warming91, 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.

Methods

Study area and focal animals

The study area is located at 550 m above sea level in the valley of the river Alm at the northern edge of the Austrian Limestone Alps (47°48′E, 13°56′N). The flock of greylag geese was introduced in the valley by late Konrad Lorenz and co-workers in 197351. The birds are unrestrained and experience natural predation mostly from red foxes (Vulpes vulpes) with losses up to 10% of the flock and 90% of the goslings per year92. The geese are individually marked with coloured leg rings and are habituated to the close presence of humans93. Individual life-history data have been monitored since 1973 and provide reliable information about the age of the breeding birds, social relationships among individuals within the flock (i.e. paired or not) as well as information on reproductive performance (i.e. breeding attempts, clutch sizes and number of fledged offspring).

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).

Table 5 Information about the breeding events (ordinal date of the first and last laid egg) and productivity (percentage of fledged goslings) per year of investigation.
Full size table

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 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 died51. 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.

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 located approx. 3 km to the South from 1990 to 2007 (13°57′20″E, 47°46′03″N) and then approx. 1 km to the North from 2007 to 2018 (13°57′03″E, 47°49′26″). The location was changed for logistic reasons. The device is run by the Office for Hydrographical Services (‘Hyrdographischer Dienst’) of the Government of Upper Austria. The station delivered daily values for air temperature (average in °C) and precipitation (sum in mm) including snow cover (sum in cm). Annual and winter (Dec–Feb) mean values were calculated for temperature and used in further analysis together with a measure of annual snow depth.

Statistical analysis

All statistical analyses were done in R version 4.1.094.

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, betareg99 and car100 packages, model predictions were extracted with effects101 and were visualised with lattice102, ggplots2103 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 package104. Additionally, we used Variance Inflation Factors (VIF) implemented in the performance105 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 collinearity106,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 ΔAICc > 2 to the first ranked model, otherwise we used model averaging with the package MuMIn108 and AICcmodavg 109 (see Grueber et al.110 for details on multimodel inference) across all models with ΔAICc < 2 and visualised a summary of model averaged effect sizes with sjPlots111. 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-linear61—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 ΔAICc < 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 ΔAICc values. Akaike weights (ωi) were calculated to assess the relative likelihood for each model considered112. Thus, ωi reflects the model’s probability given the full model list rather than only those below a given threshold of ΔAICc; added R2 values (pseudo R2 for beta regression models and marginal and conditional R2 for mixed effect models) were extracted with the performance105 package. A table of best candidate models (up to ΔAICc < 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 averaging113. We followed the guidelines provided by Burnham and Anderson112 whereby ΔAICc 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 ΔAICc < 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 variance110. 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 egg-laying 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 betareg99 and narrowed the data115:

$${x}{^{\prime}}=\frac{xx\left(length\left(x\right)-1\right)+0.5}{length(x)}$$

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 systems116,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 systems45,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 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–186).