The effect of climate change on avian offspring production: A global meta-analysis

Changes in offspring production were normally distributed and varied between −0.171 and 0.196 SDs per year from the long-term averages (Fig. 2). The funnel plot was symmetrical (Kendall’s rank test: tau = 0.002, P = 0.97), and trim-and-fill analysis indicated no directional bias in the sample of studies included in the analyses.

In 114 of 201 (56.7%) populations, offspring production decreased over time; 17.4% of these did so significantly. Conversely, in 87 populations (43.3%), offspring production increased over time, 10.4% significantly. The grand effect size was slightly, but significantly, less than zero, suggesting an overall reduction in offspring production (Table 1). An analysis of spatial autocorrelation of population effect sizes yielded a Moran I of −0.003 (95% bootstrap percentile interval: −0.110 to 0.107), which was close to the value of −0.005 expected under the null hypothesis of no effect, indicating that directional changes in offspring production were not spatially clustered across geographical regions. The analyses of other life history and demographic traits indicated that nest success decreased and egg-laying started progressively earlier. Changes in clutch size and changes in duration of egg-laying periods were nonsignificant (Table 1).

We conducted univariate meta-regressions to assess the effects of life history and ecological characteristics of populations and their environments on changes in annual offspring production (Fig. 3). Changes in offspring production were positively associated with changes in clutch size and nest success. We also found a negative relationship between temporal changes in offspring production and body mass, with relatively larger species performing worse. None of the remaining life history traits produced a significant beta coefficient. The effect of nesting in nest-boxes was marginally positively significant. Ecological factors, including absolute latitude, environment type, protection status of the study area, and a Human Footprint Index (HFI), were not associated with changes in offspring production. Similarly, long-term trends in temperature and precipitation at different periods of the breeding season did not explain variation in changes in offspring production. In particular, changes in offspring production in relation to changing rainfall were consistently negligible.

In the next step, we used multivariate generalized linear mixed models (GLMMs) to model changes in offspring production and its two proxies: changes in clutch size and changes in nest success. In all models, we used the same sets of random effects (population, species, and phylogeny), and fixed predictor variables: log body mass, number of broods in a breeding season, migratory habits, and temporal change in local ambient temperatures. We considered two types of models: one that included all possible interactions between predictors and a simpler model including only two-way interactions of change in ambient temperature with the remaining fixed factors. For rate of change in local ambient temperature, we sequentially substituted measures taken in five phenological periods (enumerated in Fig. 3). Models were compared using Akaike’s information criterion (AIC) (SI Appendix, Table S2). Models with only two-way interactions outperformed the model with all possible interactions. With respect to temperatures, temperature trends in the nestling period produced the best model of changes in offspring production. Temperature trends in the prelaying period yielded the best model of changes in nest success, and temperature trends in laying period yielded the best model of changes in clutch size. Raw data and interactions between predictors in the model of changes in offspring production and model of changes in nest success are visualized in Fig. 4.

The final model of changes in offspring production (left section of Table 2 and Fig. 4B) revealed that increasing ambient temperatures during the chick-rearing period were associated with increased offspring production in small and/or sedentary birds, whereas most migratory birds (except for the smallest ones) experienced declines in offspring production in warming climates. Increasing offspring production tended to correlate with nesting in nest-boxes (Fig. 3). Therefore, we added nest-boxes as a fixed factor to test its effect in combination with other predictors. It was insignificant (beta = −0.003, 95% c.l.: −0.022 to 0.016), and the resulting model (SI Appendix, Table S3) had ΔAIC = 6.2, indicating a poorer fit.

The final model of changes in nest success (middle section of Table 2 and Fig. 4C) suggested that the proportion of fledged nests declined in larger-bodied species. In single-brooded species, nest success was independent of climate changes. In multi-brooded bird species, nest success improved with increasing temperatures.

The final model of changes in clutch size (right section of Table 2) did not include significant interactions between predictor variables. The only important predictor was body mass, with large species producing smaller clutches.

The results of a sensitivity analysis (bar graphs in Fig. 4) showed that body mass alone captured around 40% of the uncertainty in trends of offspring production and nest success. Body mass also explained 60% of the uncertainty in the model of changes in clutch size. The other predictors had minor first-order effects. In particular, first-order sensitivity indices of changes in local temperatures were negligibly small in all final models. However, coupled effects of climate variability and life history traits were influential.

We defined offspring production as mean annual number of fledglings produced by a breeding female in a population monitored throughout the breeding season. Mean clutch size was calculated for all clutches (first, replacements and second clutches) laid by females. We defined nest success as the proportion of successful clutches (producing at least one fledgling) out of all clutches monitored across a breeding season. Second broods were defined as those initiated after successfully fledging young from a first clutch. We classified a species’ population as single-brooded or multi-brooded depending on the presence and frequency of second broods in a population (SI Appendix, Table S1 and SI Text S1). We created three environmental categories, dividing species into land birds (not associated with water), freshwater species (breeding and foraging in inland water bodies), and seabirds. For each study population, we also calculated the temporal trend in the number of breeding females present at a study area. We divided study sites into protected and nonprotected; if >50% of a study area was under some forms of legal protection (national parks, nature reserves, etc.), we classified it as protected, otherwise as nonprotected. Similarly, we divided studies with nests protected against predators and nests nonprotected; if >50% of nests were protected (e.g., concrete nest-boxes, fences around nests), we classified the population as nest-protected. Finally, we divided studies into those using nest-boxes (if >50% of clutches were laid in nest-boxes) and those using natural nest sites (monitoring open-nesting and cavity-nesting species).

The criteria for a dataset to be included in our project were: 1) a minimum study period of 15 y, during which data were collected for at least 10 breeding seasons; 2) average seasonal number of breeding females was greater than 10; 3) in multi-brooded species, females were individually marked; and 4) research had to be conducted for the whole breeding season to determine accurately the fate of successive broods, and annual offspring production of each breeding female.

Data were obtained from two major sources: 1) published information and 2) unpublished data from authors conducting long-term studies on the breeding biology of single species. 1) To locate information published in the current century concerning long-term changes in offspring production of avian populations, we undertook an exhaustive literature search in Google Scholar using the following key words: breeding, nesting, breeding success, offspring production, population productivity, climate change, and laying dates. Because we wanted to use standardized estimates of changes in offspring production across time (see below), we looked for publications that included original data. From published literature, we collected data on 72 bird populations representing 28 species (SI Appendix, SI Text S2). 2) Our literature search enabled us to identify researchers who potentially had long-term data on annual production of young in the bird populations they studied. Between October 2018 and April 2019, we contacted 313 authors of long-term studies asking them whether they would like to collaborate in our meta-analysis project. A total of 101 authors shared their data on 86 species and 129 populations (for details on searching process, see SI Appendix, SI Text S3). From both data sources, we collected a total of 201 datasets on 104 bird species (for problems with species number, see SI Appendix, SI Text S4) covering the histories of 745,962 clutches. The data were gathered in both hemispheres between 1970 and 2019 (Fig. 1), and different populations were sampled over a period ranging from 15 to 49 breeding seasons (median 26; interquartile range 20 to 35) (SI Appendix, Table S1). Literature sources, including descriptions of field methods used in the studies included in the analyses, can be found in SI Appendix, Table S4.

In addition to the data on offspring production, we collected data on laying dates, clutch size, nest success, breeding density, and proportion of females with second broods or replacement clutches. These data were not provided by all researchers, and not found in all literature sources.

Meteorological data for the study areas, including mean monthly temperatures and total monthly precipitation, were provided by authors or, for published data, obtained from the National Oceanic and Atmospheric Administration (http://www.ncdc.noaa.gov/cdo-web/datatools/findstation) or Tutiempo.net (https://en.tutiempo.net/climate/europe.html). We used data from meteorological stations closest to study areas with complete records across study years. For each area, we calculated mean temperatures (^oC) and total precipitation (mm) corresponding with each of the five breeding periods: 1) prelaying (a month preceding commencement of egg laying); 2) egg laying; 3) egg laying and incubation; 4) nestling period; and 5) whole breeding season (from the commencement of egg laying until the end of parental care). The information about timing and duration of these periods was obtained from main investigators or extracted from the literature (47, 48). We subsequently calculated the rate of change in mean temperature and total precipitation over the study period using linear regressions as an index of climate change intensity.

Data on mean female body mass for each species were obtained from Handbook of Avian Body Masses (49). The distribution of body masses was strongly skewed; therefore, we ln-transformed it for our analyses. Migratory status (full migrant vs. sedentary/partial migrant) followed the International Union for Conservation of Nature classification (50).

For each study area, we calculated a Human Footprint Index (HFI), which quantifies the degree to which humans have impacted the landscape (on 1 km² grid cells) by combining measures of human population density, buildings, crop land, pasture land, night-time lights, railways, roads, and navigable waterways (27). Values of the index range from 0 to 50. HFI data were extracted from ref. 51 using QGIS software (version 3.22.11). For each study area, we calculated a mean HFI score by overlaying the HFI raster layer with a site raster layer. Because HFI data are not available for islands, we were able to calculate the index for only 114 continental bird populations.

Analyses were conducted using R version 4.2.2 (52) and Phylometa version 1.3 (53). They were carried out in two steps: 1) we estimated grand effect sizes for a group of variables describing changes in species’ parameters and 2) using metaregression, we looked for factors explaining heterogeneity in effect sizes (hereafter ESs) from the first step.

A list of all variables used in the analysis (including the species’ parameters) is provided in SI Appendix, Table S5. All ESs studied were slopes of linear regressions, in which response variables were regressed on time (years). Slopes are intuitively interpreted as changes per annum. Annual production of offspring per female and clutch size were standardized and expressed as the number of standard deviations from the long-term mean measured across all study years of a particular population. This allowed comparison across species whose mean offspring production and clutch sizes differed.

Because of between-study heterogeneity due to variation in local environmental conditions, life histories of study species, and research methods, we used random-effects models in all meta-analyses. Sensitivity of the model due to publication bias was assessed using funnel plots and a trim-and-fill method (54).

We controlled for a potential bias due to phylogenetic relatedness among species (SI Appendix, Table S1 and SI Text S4) using a set of 5,000 equally plausible trees, downloaded from birdtree.org (55, 56). Grand ESs of changes in species’ parameters were calculated by an iterative procedure with original phylogenetic trees. For meta-regression modeling, we built a 50% majority-rule consensus tree (56) using the R package “ape,” version 5.5 (57).

Grand ESs in Table 1 were calculated using Phylometa software. We reran 5,000 iterations, each time using a different phylogenetic tree, and a subset of ESs in question, which included only a single record for each species. If a species was represented by several local populations, one was selected at random for each iteration. Grand ESs were means across all resamples, with 95% confidence limits represented by the 2.5th and 97.5th percentiles.

Meta-regression models (Fig. 3, Table 2, and Fig. 4) were run using R package “metafor”, version 3.4 (54). Each model included the same set of random effects: species and study identification (to control for nonindependence of multiple studies of the same species) and species identity associated with phylogenetic correlations (to control for nonindependence due to common ancestry).

In analyses where changes in offspring production and nest success were regressed upon different sets of climate variables, we used Akaike’s information criterion to assess the performance of competing models (58) (SI Appendix, Table S2) and to select climate variables for final meta-regression models (Table 2).

For sensitivity analysis of the final models, we used Sobol’s global variance decomposition method (59) implemented in R package “sensobol” version 1.1.3 (60). This procedure decomposes variance into fractions (Sobol’s indices), which isolate effects of predictor variables and their interactions. Pseudo-random samples for the model’s predictions were generated using a Latin hypercube sampling design. Estimates of Sobol’s indices were based on Azzini’s algorithm, using 10,000 bootstrap iterations, which guaranteed convergence under the criterion that the 95% confidence interval around the most sensitive predictor was less than 5% of its sensitivity index. For each predictor variable, we reported first-order and total-order indices. The former measures the fraction of the output variance explained by a respective predictor alone, whereas the latter includes effects of its interactions with other predictors.

L.H. and K.H. designed research; L.H. and M. Czuchra conducted literature search and collected meteorological data; L.H., P. Adamík, P. Albert, W.J.A., D.A., A.V.A., V.B., J. Bańbura, M. Bańbura, E. Barba, R.T.B., P.H.B., E. Belskii, P. B., M. Bolton, E.K.B., J. Bried, L. Brouwer, M. Bukacińska, D.B., L. Bulluk, K.F.C., I.C., M. Charter, A.C, R.C., D.C.D., F.d.L., A.S.D.G., V.C.D., H.D., M.J.D., T.E., L.M.E., Y.O.E., J.A.F., S.I.G., E.Y.G., M.G., M.P.H., J.P.H., P.K., J.K., W.D.K., H.K., M.K.-A., E.K., I.K., M.K., S.K., B.D.K., X.L., M.P.L., A.L., A. Margalida, A. Marzal, A. Millon, A.P.M., V.C.N., J.T.N., A.N., B.O., D.O, M.Ö., R.A.P., H. Pietiäinen, V.P., J. Porkert, J. Potti, H. Pöysä, T.P., J. Prop, P.Q., J.A.R., P.-A.R., R.N.R., A.R., D.R.R, I.E.S., D.A.S., M.S., J.C.S., F.S., T.S., D.V.S., J.S., P.M.T., M.T., J.Tolvanen, J. Török, M.v.d.P., L.V., M.E.V., D.F.W., N.T.W., J.W., K.W., A.G.W., A.W., D.W., M. Z. performed research and contributed the data; Z.J. calculated Human Footprint Indices; K.H. analysed the data; L.H. and K.H. wrote the manuscript, and all authors edited the manuscript. Valuable scientific contributions throughout the writing and editing process were done by D.A., J. Tolvanen, A. Millon, W.D.K., P.B., M.G., M.S., H. Pöysä, M.Ö., M.T., M.v.d.P., D.A.S., and A. Margalida.

The authors declare no competing interest.