Many small climate change impacts presage rapid population extinction in a common iconic bird

Domain: biology
Published: 2026-03-31
Editorial outcome: minor_revision

Abstract
Population declines in common bird species abound, but it is unclear if these declines presage extinction. The consequences of climate change are particularly elusive, as there are typically both positive and negative effects on the seasonal- and life-cycles. We examine climate sensitivities of a population of a common and loved Australian bird, the superb fairy-wren Malurus cyaneus, with three decades of year-round fine-scale measurement of reproduction, survival and immigration. Population growth has been affected by many aspects of climate, but declines are most strongly linked to low fecundity in dry springs and reduced survival during unusually warm winters and preceding hot summers. Our comprehensive integrated population model is effective at summarizing these effects and predicting recent population volatility, suggesting it will be useful for assessing future population fluctuations and viability. A high risk of population extinction in the next 30-40 years is predicted. Accumulation of many small effects of climate change causes this decline. Hence a common species of Least Concern is in imminent danger, with virtually no time or way to prevent it. Such predicaments can only be detected in very well-studied species, and superb fairy-wrens may be the canary in the coal mine.

Introduction
Extinction risk is highest for rare species, as their small number and restricted distribution make them vulnerable to local environmental change and chance events. However, this does not mean that common species are safe. Historically, some common species have gone extinct very rapidly. Notably, the passenger pigeon (Ectopistes migratorius) was once the most abundant bird species (~3–5 billion individuals) but went extinct only a few decades after the first concerns about their declining numbers were raised1. Concerningly, ample evidence has emerged that many common species are rapidly declining2–4. Whether and how fast the current decline of common species will lead to extinction is unclear. However, this question is urgent, as the extinctions of common species are likely to have disproportionate impacts on ecosystems5,6. Furthermore, rapid declines leave little time for conservation action to prevent extinction and its cascading impacts.

Most historic rapid extinctions have been caused by human-driven overexploitation, deliberate or inadvertent introduction of predators, and habitat loss2. The increased tempo of anthropogenic climate change is exacerbating some of these threats, but also posing new threats7,8. However, whether climate change can also lead to rapid extinction of common species is unknown and more challenging to study. This is because climate (change) affects many physiological and behavioural traits, as well the environment (resources, predators, phenology, habitat). Consequently, climate change can affect a range of vital rates across the seasonal- and life-cycles of organisms in positive and negative ways, making it hard to identify its overall impact. This problem is exacerbated by the fact that most studies only monitor species during part of the year, meaning effects could be overlooked or obscured by opposing effects in non-monitored parts of the seasonal cycle. Finally, the timescale of much research is insufficient to identify the impacts of extreme events, a key component of climate change9. In this study, we show that many small climate change effects can accumulate to drive rapid population extinction in a common iconic bird.

Results
Superb fairy-wrens
We used field data on a common resident songbird of south-eastern Australia, the superb fairy-wren (Malurus cyaneus), to first identify the demographic pathways whereby climate affects population dynamics via vital rates, and then use this to quantify future extinction risk under climate change. Fairy-wrens are an ideal study species because (i) their vital rates can be measured precisely with colour banding and intense survey10,11; (ii) previous studies found clear climate impacts on their life history traits (e.g., breeding phenology11 and chick body mass12) and vital rates (e.g., adult survival rate in the non-breeding season10); and (iii) declines have been observed in some populations of this species13,14. We have studied reproduction and survival of a fairy-wren population in the Australian National Botanic Gardens, Canberra, for three decades, achieving year-round fine-scale censuses that characterized the weekly timing of mortality, all nesting attempts, and all immigration15. The study area is an elevated inland site with strong seasonal variation in temperature, and alternation between drought and more mesic conditions. The population fluctuated throughout the 1993–2022 study period, though it generally declined until a burst of growth in the final two years (Fig. 1A).

Male numbers were consistently greater than female numbers (Fig. 1A). This male excess is due to cooperative breeding, as males are tolerated on their natal territory, and may help the breeding pair through provisioning of offspring and nest defence for up to 6 years16. If males become dominant, they almost always do so on their natal territory or on an immediately adjacent territory. By contrast, juvenile females must disperse in their first year to obtain a breeding vacancy, and often move over large distances17,18. Females achieve this dispersal in two ways. Females fledged before mid-December typically disperse before the end of the breeding season to a subordinate role in a foreign group. Females fledged after that date overwinter in their natal territory. All young females enter a scramble for a breeding vacancy at the start of spring, except for early dispersers that have already inherited the territory on the death of the incumbent. Early dispersers are also advantaged over overwintering natal birds because they can form an alliance with a helper male and fission the territory in which they have settled, and are much more likely to obtain a vacancy near where they have overwintered17. Females that do not obtain a vacancy in the breeding season after they fledge usually die, generating the male bias in adult sex ratio.

These social influences led us to consider three phases of the annual cycle (Fig. 2) during which climate impacts may occur. I. the recruitment phase (mid-November to March), where territories accumulate extra members through reproduction and immigration of young females; II. the non-breeding phase (Apr–Jul), where there is little movement but decreased survival; and III. the scramble phase (August to mid-November), where the number of breeding territories is settled through fission and fusion of existing territories, and the replacement of breeding females that have died.

Integrated population model
Integrated population models (IPM) combine individual- and population-level data into a single statistical model that also describes their interrelationship governed by demography and population dynamics19. We first developed an IPM to estimate all vital rates and annual population growth rates (PGRs). To facilitate understanding, we only present a general non-technical description of the IPM here; a technical description is in “Methods”. Our model considered five different stage classes to account for substantial differences in (climate impact on) vital rate among types of individuals (Fig. 2). We considered two sexes (male and female) each with two age classes (juvenile and adult, the latter being members of groups attending a nest). Juvenile females were split further into two stage classes that differed in their recruitment prospects: natal females that stayed in their natal territory at the end of the recruitment phase, and non-natal juveniles that moved out of natal territory during the recruitment phase.

Transitions between these stage classes was governed by six stage-specific vital rates: apparent survival rate \phi, the probability that an individual survived (and remained in the study population) over a given phase; immigration rate \delta, the number of immigrants in a given phase relative to the number of adult females (or males) at the end of the scramble phase; fecundity \gamma, the number of juveniles produced by each female that survived and stayed in the study area in the recruitment phase; offspring sex ratio \tau, the proportion of females among observed juveniles; female natal philopatry rate \pi, the proportion of juvenile females produced in the study population that stayed in their natal territories in the recruitment phase; and the probability of males being the dominant in their territory \psi (detailed algebraic definitions in “Methods”). Furthermore, vital rates (and climate impacts on them) were specific to each annual phase. Finally, each vital rate was modelled as a binomial or Poisson random variable to account for uncertainty due to small population size.

Demographic pathway of climate impacts on population dynamics
To identify the main demographic pathways via which climate affects population dynamics, we conducted a retrospective path analysis20,21 to quantify how climate variables affect annual PGR through vital rates. This effect of climate on annual PGR via each vital rate is given by the product of how much that vital rate is affected by the climate variable and how much the change of annual PGR is affected by a change in vital rate (Supplementary Fig. 4). Based on detailed previous studies on this population, we explored the effects of five climate variables on all relevant vital rates in the different annual phases (Supplementary Fig. 1): spring rainfall (Aug–Nov)22,23, summer maximum temperature (Dec–Feb)10,11, summer rainfall (Dec–Feb)11, and winter minimum and maximum temperatures (Apr–Jul)10. Vital rates were also density-dependent, but as this did not significantly affect our conclusions; we have confined the description of density-dependence to “Methods”.

We found 11 demographic pathways via which climate affects population dynamics of superb fairy-wrens (Fig. 3). More spring rainfall in a given year was associated with a higher annual PGR through increased fecundity and immigration rate of juvenile females in the recruitment phase (Fig. 3 and Supplementary Figs. 2A, 2B and 4). The breeding season is very long, females can produce up to four successful broods, and the clutch size can also vary. As reported previously, prolonged breeding, increased clutch sizes and promoted fecundity are strongly positively associated with high spring rainfall11,22,23. Hence both the number of natal and early-dispersing young increased after high rainfall. In extreme droughts, recruitment can be severely limited and the population therefore declines as breeding females die and are not replaced. This often leads to fusion of neighbouring territories17.

The effects of temperature on annual PGR are also important. Population growth was limited by above-average winter and summer maximum temperatures, and by below-average winter minimum temperatures, mostly through reduced adult and juvenile survival rates in the non-breeding phase (Fig. 3 and Supplementary Figs. 2–4). Increased maximum temperatures in the winter worsen the severity of an annual episode of reduced survival that occurs just after the winter solstice, when temperatures are coldest10. The survival may be particularly low when warm conditions are followed by cold snaps10. Finally, there appear to be deleterious effects of warm summers on survival in the subsequent year10. All of the effects on adults reported in a previous study10 are recovered in the current study with more years of data (Fig. 3). Juveniles similarly may be affected by temperatures as adults, especially considering that juvenile survival is usually affected by weather more than adult survival24.

Future population extinction risk under climate change
We assessed the fate of the population under continued climate change. To first validate the predictive value of our model, we fitted a model to our data from 1993 to 2015 and determined how well it predicted the population changes from 2016 to 2022. These last years of the study are particularly suitable for validation as they included climate extremes, and unprecedented volatility, including both the sharpest population decline and the two years when population growth was greatest (Fig. 1B). This cut-off was also convenient because 2016 was the hottest year on record in south-eastern Australia, making it a critical test case of our accuracy in projecting to novel extreme conditions. Both the projected sharp decline and subsequent recovery were qualitatively in agreement with the observed population fluctuations from 2016 to 2022 (Fig. 4A). Although the model slightly underestimates the depths of the recent decline and the subsequent rebound (see next section for more discussion), the characterization of both these events support the view that our IPM could project the population size reliably.

Next, we used the model that was validated and parameterized with historical field data to projected annual population size from 2022 to 2100 by using forecasted climate variables under scenarios of low, intermediate and very high greenhouse gas emission (respectively, IPCC Shared Socioeconomic Pathways 1–2.6, 2–4.5 and 5–8.525). The forecasted maximum and minimum temperatures in summer or winter increased in most emission scenarios and the rate of warming amplified from low, intermediate to very high emission scenarios (Supplementary Table S1). Although spring rainfall generally showed no directional change across years in any of the emission scenarios (Supplementary Table S1), its effect on recruitment can still compound the deleterious effects of temperature on survival. In addition, as a control, we also projected population dynamics under a scenario with no further increased emission-driven directional climate change (see “Methods”). To appreciate the impact of climate change on population dynamics relative to other sources of variability and uncertainty, our ensemble projections averaged over all 21 climate models and estimation uncertainty in all model parameters, and included demographic and environmental stochasticity (see “Methods”).

Under the scenario of no further climate change, superb fairy-wren population size was forecasted to fluctuate with a low risk of population extinction (22.9% probability of female population size decline to zero in 2100; Fig. 4B). In contrast, although the projected population trajectories differed somewhat among different greenhouse gas emission scenarios, climate change in general led to dramatic and rapid population decline and consequently all projections went extinct, even in the low emission scenario (Fig. 4C, D). Additionally, the 50% risk of population extinction occurred progressively earlier from low, intermediate to very high emission scenarios. Under the intermediate and very high emission scenarios, the fairy-wren population size was forecasted to go extinct between 2059 and 2062. This implies that we have only 3–4 decades to prevent likely local extinction, as current global emission reductions exceed the low scenario26.

Are there avenues for population rescue?
The alarming projections of extinction raise the question of whether any actions or factors are likely to prevent or delay extinction. Three possibilities warrant discussion.

The first is genetic adaptation under selection. Although our population projection assumes that the effects of climate on vital rates in the future will be the same as in our study period, evolutionary rescue in vital rates may occur and reduce the risk of extinction. For example, after an extreme cold event, green anole lizards (Anolis carolinensis) become more resistant to cold and show gene changes that are relevant to thermoregulation27, which could help them to survive future cold winters better. A recent synthesis of long-term data from wild birds and mammals also suggests that the potential for rapid adaptation was greater than hitherto believed, but the genetic variance for fitness in superb fairy-wrens was at the lower end of the species considered28. It has also recently been argued that changes in the direction of selection associated with climate fluctuations can impede evolutionary rescue29.

The second is the immigration of warm-adapted genotypes. The breeding distribution of superb fairy-wrens spans more than 20 degrees of both latitude and longitude, and most of Australia’s altitudinal range30,31. The population under study experiences a continental climate exacerbated by relatively high altitude, meaning it experiences colder winters and drier conditions than in much of the species range. However, there is likely connectivity to populations adapted to warmer (and possibly drier) environments. The extreme philopatry of male fairy-wrens retards gene flow, but the obligate dispersal of young females has been shown to introduce novel genotypes into our study population32, so rescue by gene flow is possible.

Third, climate change affects all the species in a habitat. We have evidence that during the period over which superb fairy-wrens have declined there has been a far steeper decline in numbers of their most significant nest predator, the pied currawong (Strepera graculina). Currawong territories are large and hence the population is more difficult to enumerate, but has potentially declined from a peak of more than 20 territories to just 5 in 2020. Nest predation, the most significant cause of nesting failure, has declined accordingly, which may help explain the slight underestimation by our IPM of the population surge in the last two years reported here. Superb fairy-wrens may thus be rescued by climate-driven reductions in predators, but this hypothesis awaits further testing.

Accumulation of many small climate change impacts
Despite some potential for natural rescue, the dire model projections cause considerable alarm. Population extinction is likely to happen very fast; within the next 30–40 years in the intermediate and very high emission scenario. This leaves little time for conservation actions, a realization that is striking given that superb fairy-wrens are very common in a range of habitats and have always had a conservation status of Least Concern33. We previously reported that some aspects of their life-cycle are influenced by climate change10,11, but this has not led to concerns about their future viability thus far (also not by us, the researchers). Only by integrating the many impacts of climate change across the seasons and the full life-cycle, we found that although many individual impacts are small or moderate, together they are expected to have a catastrophic cumulative effect, even in the low-emission scenario.

We identified 11 demographic pathways via which a suite of climate variables operate at different times of the year, with different time-lags, and differential impacts on males and females. This knowledge is not only important for making reliable projections but also for identifying the nature and limits of conservation actions, as many demographic pathways make it hard to come up with solutions. Mitigating climate change impacts is challenging, as it cannot be solved locally and can only be buffered or compensated by improving other environmental conditions. However, mitigating a wide diversity of impacts is even more challenging. Given that most conservation success stories involve species in which a single demographic cause of decline was resolved34,35, we should be modest and realistic (i.e., pessimistic) about what conservation can achieve in complex cases like superb fairy-wrens, especially within such a short timeframe.

Population decline and extinction of common species
The deleterious effects of climate change, and particularly of increased winter temperature on our population of superb fairy-wrens is surprising, as a priori an increase in winter temperatures should move these endothermic birds closer to their thermoneutral zone and lessen their energetic burden. Although we lack the relevant data, which would be extremely expensive to obtain, we suspect the problem posed by sudden warm spells followed by cold spells10 may negatively affect the ectothermic arthropod prey exploited by the birds, possibly by provoking emergence and activity at the wrong time of year. In that sense, our data resemble the large reduction in populations of common bird species in both Europe and North America3,4, where the greatest declines are in species that are wholly insectivorous or feed their young on arthropods36–38. Like bird populations, some arthropod declines have been caused by habitat simplifications through agricultural intensification, and the use of novel insecticides39–41. However, there is also evidence that huge arthropod declines can occur in environments that are not subject to either of these influences. These have been attributed directly to global warming42. The decline in fairy-wrens we report on here is not in agricultural habitat, but their habitat is experiencing the deleterious effects of warming.

Our studied population of superb fairy-wrens is predicted to crash under all climate change scenarios. This result echoes the recent reports that climate change will limit population growth and lead to local extinction of wild animals43,44. For example, a recent study on African wild dogs (Lycaon pictus) found that climate warming will lead to population collapse in both intermediate and very high emission scenarios45. Although it is well-established that climate change threatens the population persistence of wild animals8,46, our study uniquely provides demographic evidence that climate change poses an urgent threat to common species. As ecosystems are likely to be disproportionately impacted by the extinctions of common species5,6, more conservation attention to common species is needed in the future viability studies of wild animals.

Our study highlights the value of year-round monitoring and comprehensive analysis of all vital rates throughout the seasonal- and life-cycles. We expect it to be a common feature in many species that climate change affects them in many ways. However, most species have not been comprehensively studied year-round for decades, meaning that such effects also act during parts of the year when researchers are not in the field, and thus will often not be identified. In that sense, superb fairy-wrens may be the canary in the coal mine.

Methods
Ethics
The study period was covered by multiple permits, including the approvals from the Animal Experimentation Ethics Committee of the Australian National University (to A. Cockburn; A2022/20 (approved on 10 August 2022), A2019/23 (approved on 9 August 2019), A2016/22 (approved on 10 August 2016), A2013/27 (approved on 9 August 2013), F.BTZ.42.10 (approved on 10 August 2010), F.BTZ.06.07 (approved on 10 August 2007), F.BTZ.23.04 (approved on 10 August 2004) and F.BTZ.16.01 (approved on 10 August 2001)) and the Australian National Botanic Gardens (to A. Cockburn and H. L. Osmond; 2020/21-1 (approved on 27 June 2020), 2013/4-1 (approved on 9 December 2013), 2008/9 (approved on 1 January 2008), 2006/1651 (approved on 1 January 2006) and 2004/1549 (approved on 1 January 2004)). In addition, A. Cockburn held permits from the Animal Experimentation Ethics Committee of the Australian National University between 1993 and 2001, and both A. Cockburn and H. L. Osmond held approvals from the Australian National Botanic Gardens from 1993 to 2004.

Overall approach
Our two aims are (i) to quantify how climate affects population dynamics via vital rates in a population of superb fairy-wrens (Malurus cyaneus), and (ii) to project population viability under future climate change. To address the first aim, we developed an integrated population model (IPM) to estimate vital rates and annual population growth rate (PGR) throughout the annual cycle. Subsequently, we conducted a retrospective path analysis to measure the effects of climate on annual PGR through each vital rate. To address the second aim, we extended the previous IPM by including the historical effects of climate variables on each vital rate and combined it with forecasted future climate downscaled to the study area, to project population viability under different future climate change scenarios.

Model system
Superb fairy-wrens are a common resident passerine of south-east Australia. A population of fairy-wrens in the Australian National Botanic Gardens, Canberra (35°16′S, 149°6′E), has been studied continuously since 1988. The study area was expanded substantially in 1991/1992, therefore we focus here on the data from 1993 to 2022. The study area encompasses a total of ~60 hectares divided almost equally between plantation of Australian native shrubs and trees, and natural Eucalyptus woodland. The study population is a sample of a much larger population with free movement across its boundaries to unmonitored territories; we included a given territory if >50% of the territory was within the study area. There were between 36−87 fairy-wren territories in a given year (Fig. 1A).

Superb fairy-wrens forms socially monogamous pairs and live in year-round territories. All birds in the population are individually colour-banded, either as nestlings if they were born in the study area (at 5-7 days old), or as juvenile or adult immigrant caught by mist-netting. We aimed to observe every bird at least weekly throughout the year. If a bird was not seen on routine censuses, we made further attempts to find it, resulting in near-complete detection rates10. The breeding pair on each territory can either breed unassisted or be helped by up to 5 male helpers in this cooperatively breeding species15. During the breeding season (Aug–Mar11) we monitored every nesting effort by every female at least every three days. Although the adult sex ratio is male biased, there is negligible bias in the sex ratio of hatchlings47, and adult males and females have similar survival rates48.

Male offspring are extremely philopatric birds occupying year-round territories, and generally live and die on their natal or an immediate neighbouring territory16. Unbanded males are very rare and largely occur on peripheral territories. Both dominant and subordinate males are sexually active in their first breeding season15. In sharp contrast, female offspring always disperse away from their natal territory in the first year of their life, often moving more than two territories (few over 5 km)17,18. Females start breeding in the first breeding season after they fledge and thereafter only rarely move more than one or two territories17. Therefore, virtually all female immigrants are in their first year of life.

Aside from the first year emigration and immigration of females before they gain a breeding vacancy, the loss of an individual from our study population is generally associated with mortality. The death of a breeder male or female is usually associated with rapid replacement by a helper or a disperser, respectively. Occasionally, birds on the periphery are lost from the census population because of short-distance breeding dispersal or fission of breeding territories, but are censused again when the territory to which they have moved fuses with a censused territory, or they return when their mate dies. This affects just 2% of the 2088 individuals considered (45/2088), of which most (40/45) are females.

Climate data
We used data on daily weather from an Australian Bureau of Meteorology automatic recording station at Canberra Airport (~8 km east of the study area and at similar altitude; http://www.bom.gov.au/climate/data). As the study site has a highly seasonal continental climate, and extreme weather is likely to affect the vital rates, we explored the effects of mean daily maximum temperature and mean daily rainfall in the austral summer (from 1 December to the following 28 or 29 February; hereafter ‘summer maximum temperature’ and ‘summer rainfall’), and mean daily minimum temperature in the austral winter (from 1 April to 31 July; hereafter ‘winter minimum temperature’). The effect of mean daily maximum temperature in winter (hereafter ‘winter maximum temperature’) was also included since we recently reported that higher winter maximum temperatures were associated with higher adult mortality in the non-breeding season10. Additionally, we tested the effect of mean daily rainfall in the austral spring (from 1 August to 30 November; hereafter ‘spring rainfall’) on fecundity because more spring rainfall increases the availability of food through benefiting plant growth and increasing abundance of arthropods for breeding superb fairy-wrens16. Although climate at shorter time intervals within a season (e.g., the coldest month of the year) might be expected to have strong effects on vital rates, preliminary analyses using averaged daily temperature and rainfall summarized at monthly intervals resulted in similar or weaker effects compared with when climate conditions were averaged across whole seasons.

To generate future values of these five climate variables (spring rainfall, summer maximum temperature, summer rainfall, and winter minimum and maximum temperature) from 2022 to 2100, we obtained monthly maximum and minimum temperatures, and rainfall, from all available global climate change models (https://ds.nccs.nasa.gov/thredds/catalog/catalog.html) of the Coupled Model Intercomparison Project Phase 6 (CMIP6)25 with a spatial resolution of 0.25° latitude-longitude (~28 km), and chose the model grid point in which our study area is located. For the greenhouse gas emission scenarios of each climate model, we considered three widely used Shared Socioeconomic Pathways (SSPs): (1) an optimistic low scenario with peak emissions in 2020 which was based on SSP1–2.6; (2) an intermediate scenario with peak emissions in 2040 which was based on SSP2–4.5; and (3) a worst-case unrestrained very high scenario which was based on SSP5–8.5. In our projections, we randomly sampled climate time series from one of 21 global climate change models (Supplementary Table S1) that provided monthly maximum and minimum temperatures, and rainfall for all three greenhouse gas emission scenarios. We did not consider intra-model uncertainty in our projections, as it was not available for all 21 models and is generally much smaller than the variability among climate models49,50. The forecasted climate variables varied considerably among climate models for a given greenhouse gas emission scenario, but the tendencies for changes in the projection period are similar. As ensemble approaches tend to have better forecast performance51, we therefore focus on ensemble projections. To facilitate the understanding of the change of annual climate and seasonal climate variables in each climate change scenario, we summarized the change in each climate variable at both global and local (Canberra) scales in Supplementary Table S1.

For each set of the forecasted climate variables, their mean differed with the mean of corresponding historical climate variables in the overlapping years (2015-2021; N=7), potentially due to the difference in scales: historical climate data from one local weather station versus forecasted climate data at a large grid cell level (0.25° × 0.25° latitude-longitude; i.e., ~760 km2). To avoid the potential influence of different values of the mean between historical and forecasted climate data, and between different climate models for each greenhouse gas emission scenario in the forecasted climate data, we calibrated the values of forecasted climate variables from 2022 to 2100. Temperature variables and rainfall variable were each modified differently52. For each temperature variable, we subtracted the difference between forecasted temperature variable and historical temperature variable in the overlapping years. Whereas, for each rainfall variable, we multiplied the quotient of historical and forecasted rainfall values in the overlapping years. In this way, we retained the distribution and temporal change of the forecasted climate data.

Finally, to assess the impact of future anthropogenic climate change in above scenarios, we also generated a climate scenario of no further increased emission-driven directional climate change where the climate will not further change from 2022 to 2100. In the ‘no further climate change’ scenario, we generated values randomly for each forecasted climate variable from 2022 to 2100 using a Gaussian distribution with the mean and standard deviation of the corresponding historical climate variable from 1993 to 2021.

Integrated population model
To estimate vital rates and annual population growth rates (PGRs), we developed an integrated population model (IPM) which incorporated multiple types of data in a joint likelihood to estimate parameters. This allows for synthesizing data sources, efficient use and internal consistency of available information, and comprehensive error propagation throughout the population model19. Our IPM model included four submodels linked to the relevant datasets: (i) a population submodel utilising the population count data to inform historical changes in population size, (ii) a survival submodel using the individual capture-resighting data to estimate the apparent survival of different stage classes, (iii) a recruitment submodel utilising the recruitment data to estimate fecundity, and (iv) an immigration submodel using the immigration data to determine immigration rates. Below, we describe each sub-model and then the overall model.

Population submodel. The population submodel describes how female and male population size change throughout the annual cycle, and links the population count data with the population size variable in the model. Our model consists of a life-cycle structure that includes three phases in the year (I–III, Fig. 2) to accurately account for the fact that climate acts on different vital rates in different parts of the year. Notwithstanding, for future projections we can focus on the dynamics at the end of the scramble phase (15 November) from one year (t) to the next (t+1), which means that our stage and phase-structured model can be summarized by a single equation for the dynamic of adult females or males from year to year Eqs. (1), (2). At this point, all juveniles have become adults or died, and our model simplifies to a two-sex population model of the dynamics of adults only. We assumed that the sex-specific (S, which could be either F(emale) or M(ale)) population size in the first year followed a Poisson distribution such that N S,I,1 ~ Poisson(y S,I,1), in which y S,I,0 is the sex-specific population count of number of adults on 15 November in the first year of study period (i.e., 1993), 87 for the females and 181 for the males. From the second year, we assumed that female and male population size was governed by population size in the previous year, survival, reproduction, and immigration such that

N F,I,t+1 = Poisson(N F,I,t * (\phi F,I,t * \phi F,II,t * \phi F,III,t + \gamma t * \tau t * \pi t * \phi L,II,t * \phi L,III,t + \delta F,I,t * \phi O,II,t * \phi O,III,t + \delta F,III,t)) (1)

and N M,I,t+1 = Poisson(N M,I,t * (\phi M,I,t * \phi M,II,t * \phi M,III,t + \psi t * \gamma t * (1-\tau t) * \phi J,II,t * \phi J,III,t + \delta M,III,t)) (2)

(see Supplementary Table S2 for a description and estimates of vital rates \phi, \gamma, \tau, \pi, \delta and \psi). From this population model, we can estimate annual PGR P S,t as the change in population size N from one year to next for males or females such that

P S,t = N S,I,t+1 / N S,I,t (3)

In this model, we linked the population count data y with the population size N using a Poisson distribution such that

y S,I,t ~ Poisson(N S,I,t) (4)

The probability of males being the dominant in their territory \psi was informed by the counts of dominant males y D,I,t among all males y M,I,t with a binomial model such that

y D,I,t ~ Binomial(y M,I,t, \psi I,t) (5)

Survival submodel. We used the Cormack-Jolly-Seber model53 to estimate apparent survival over each of the three phases from individual capture-resighting data. The data were collated into an m-array format, which followed multinomial distributions. The cell probabilities of the m-array were calculated from stage-, phase-, and year-specific apparent survival probabilities and stage-, and phase-specific detection probabilities. We assumed that detection probabilities are constant across years due to intense survey effort (detection probabilities are very high, except for the dispersing juvenile females; Supplementary Table S3).

Recruitment submodel. We estimated fecundity \gamma t using recruitment data (i.e., counts of adult females and different juveniles) such that

y JMF,II,t ~ Binomial(y JMF,II,t + y F,I,t, \gamma t / (\gamma t + 1)) (6)

in which y JMF,II,t is the sum of the counts of juvenile males and females. In our model, we assume that subordinate males do not contribute to reproduction, although in reality they can obtain substantial paternity54,55. However, the usurped paternity from dominant males does not affect the number of juveniles produced per group/female (i.e., fecundity), and hence we can ignore it in our analyses. Additionally, we estimated the offspring sex ratio \tau using a binomial model such that

y JF,II,t ~ Binomial(y JMF,II,t, \tau t) (7)

in which y JF,II,t is the counts of juvenile females (natal and non-natal juveniles combined). The female natal philopatry rate \pi was also estimated using a binomial model:

y L,II,t ~ Binomial(y JF,II,t, \pi t) (8)

in which y L,II,t is the counts of natal juvenile females.

Immigration submodel. We linked the immigration rates of females and males to immigration data. Female immigration rates were estimated in the recruitment phase \delta F,I,t such that

y GF,II,t ~ Poisson(N F,I,t * \delta F,I,t) (9)

and in the scramble phase \delta F,III,t such that

y GF,I,t+1 ~ Poisson(N F,I,t * \delta F,III,t) (10)

in which y GF,II,t and y GF,I,t+1 are counts of female immigrants at the end of recruitment phase and scramble phase (1 April and 15 November), respectively. Immigration rate of males were only estimated in the scramble phase \delta M,III,t such that

y GM,I,t+1 ~ Poisson(N M,I,t * \delta M,III,t) (11)

in which y GM,I,t+1 is the counts of male immigrants.

Model implementation. We analysed the IPM in a Bayesian framework using Markov Chain Monte Carlo simulations, which were implemented in JAGS56 using the R package jagsUI57,58. We used vague priors Normal(0, 0.01) for all mean parameters, Gamma(0.01, 0.01) for all precision (inverse of variance) parameters, and Uniform(0, 1) for all detection probability parameters. By running three independent chains with different starting values for 10,000 iterations, with a burn-in of 8,000 iterations, we obtained 6,000 posterior samples without thinning. The convergence of model chains was checked using R-hat and Gelman-Rubin diagnostic statistics59. The R-hat statistics for each parameter were well below 1.03, and thus the chains were well mixed.

Demographic pathway of climate impacts on population dynamics
To properly estimate the effects of climate on annual PGR, and the demographic pathway (via which vital rate) that influence annual PGR the most, we applied a retrospective path analysis20 modified from a recent study21. Our retrospective path analysis model could account for uncertainty in the model inputs by fitting every linear model in the path diagram to each sample in the posterior distribution of each vital rate estimated by the IPM21. The model included paths between (a) climate/density variables to each relevant vital rate and (b) these vital rates to annual PGR (Supplementary Fig. 4). We estimated the coefficients of the first kind of paths in the same way as used by conventional path analysis. However, as annual PGR is completely determined by all vital rates, the path coefficients from vital rates to annual PGR were not estimated by the conventional path analysis, but derived mathematically from our population model (see below).

Density dependence. We considered the density-dependence of vital rates in both females and males, but in different ways. For the females, density-dependence is thought to act in the scramble phase and was incorporated by making the transition (survival) rate of females into a breeding position dependent on the number of vacant territories per female competing for such vacancy. Females that do not obtain a breeding position usually die and there are very few female supernumeraries in this population17,18. The number of vacancies is associated with the number of adult females that died in the previous recruitment and non-breeding phase, while the number of females competing for vacancies includes natal and non-natal females that survived the non-breeding phase. Thus, we calculated the breeding vacancy availability index as the number of adult females that died in the recruitment and non-breeding phase, divided by the number of natal and non-natal females that survived the non-breeding phase. We also included the effect of male group size (n g,t; i.e., the average number of males per female at the start of the scramble phase (1 August), which was calculated as (y M,III,t + y J,III,t) / (y F,III,t + y L,III,t + y O,III,t)) on the survival rate of non-natal juvenile females in the scramble phase. More helper males in the territory in the scramble phase enhances the opportunity of obtaining a breeding position for those non-natal juvenile females by defending the territory and providing food to the breeding pair (or group) in the territory. Thus, we included the effect of male group size on the survival rate of non-natal juvenile females in the scramble phase.

For the males, density-dependence is determined by the female dynamics, as the probability of males being the dominant in their territory \psi is driven by the number of territories/groups, which we assume is limited by the number of females. Therefore, we equated the number of dominant males to the number of dominant females, which implies that \psi = y F,I,t / y M,I,t, with \psi <= 1. Finally, we assume that no other vital rates are affected by group size or density-dependent.

Effects of climate/density variables on IPM estimated vital rates. We tested the effects of relevant climate/density variables on each vital rate (Supplementary Fig. 1) by extending our previous IPM through including climate/density variables. To reduce the number of variables considered in the final retrospective path analysis model, we first conducted a variable selection procedure. We fitted a series of models, in which each vital rate (fecundity, sex- and phase-specific survival rates, and female phase-specific immigrations rates, see Supplementary Fig. 1) was fitted as a logistic or log-linear function of a single climate variable. The retrospective path analysis model only included a variable for a given vital rate if the 90% credible interval of the slope coefficient excluded 021. The rationale for using a 90% instead of 95% credible interval for variable selection is that this reduces the chances of overlooking climate effects (false negatives). Although some climate/density variables were significantly inter-correlated (e.g., spring rainfall and breeding vacancy availability index; Supplementary Table S4), none of the equations for a given vital rate included climate/density variables that were correlated significantly at the alpha = 0.05 level (Supplementary Fig. 4; Supplementary Table S4).

We modelled each estimated vital rate in the IPM with a mean \alpha and a process error \epsilon, where \epsilon t followed a Gaussian distribution with mean 0 and standard deviation \sigma. The process error estimates the among-year variance of a given vital rate that is not due to any modelled effects of climate/density variables (i.e., the fixed effects). Stage- and phase-specific apparent survival rate \phi t was modelled using logistic regression with a logit link function such that

logit(\phi t) = \alpha + \beta * X t + \epsilon t (12)

in which \beta denotes the estimated slope and X t represents annual values of a given climate/density variable, or a combination of climate/density variables with \beta or X t in different vectors.

Female phase-specific immigration rate \delta F,t or fecundity \gamma t was modelled using Poisson regression with a log link function such that

log(\delta F,t or \gamma t) = \alpha + \beta * X t + \epsilon t (13)

Path coefficient of vital rates to annual PGR. In our retrospective path analysis, the path coefficient of a given vital rate to annual PGR measures the contribution of the change of that vital rate to the change of annual PGR. This specific path coefficient can be approximated by the sensitivity of the annual PGR of the model to change in the vital rates. To quantify the sensitivity of annual PGR, we re-expressed annual PGR P S,t as (combining Eqs. (1)–(3)):

P F,t = \phi F,I,t * \phi F,II,t * \phi F,III,t + \gamma t * \tau t * \pi t * \phi L,II,t * \phi L,III,t + \delta F,I,t * \phi O,II,t * \phi O,III,t + \delta F,III,t (14)

and P M,t = \psi t * \gamma t * (1-\tau t) * \phi J,II,t * \phi J,III,t + \delta M,III,t (15)

and determined the partial derivative of P S,t with respect to each vital rate. For example, for male superb fairy-wrens, one unit change of immigration rate \delta M,III,t will contribute exactly one unit change of annual PGR P M,t (i.e., \partial P M,t / \partial \delta M,III,t = 1 as can be seen from Eq. (15)). However, one unit change of adult male survival rate in the recruitment phase \phi M,I,t will contribute less, specifically \partial P M,t / \partial \phi M,I,t = \phi M,II,t * \phi M,III,t \approx 0.72 (see Eq. (15) and Supplementary Table S2). This means that, if climate change decreases \phi M,I,t by 0.01 (1% lower survival), this will only lower annual PGR by 0.72%.

Future population extinction risk under climate change
We projected population dynamics under three climate change scenarios and contrasted these with the scenario of no further climate change to determine the impact of future anthropogenic climate change. For each climate change scenario, we specifically accounted for key sources of uncertainty: uncertainty among the 21 climate models, uncertainty in model parameter estimates, demographic and environmental stochasticity (both due to climate effects modelled explicitly and non-climate effects modelled as residual interannual noise). Demographic stochasticity reduces population growth rate at small population size60, which is particularly relevant for studies quantifying risk of extinction.

As the first step, we extended our previous IPM by including the effects of influential climate/density variables on vital rates, as identified in the retrospective path analysis (Supplementary Fig. 4). Thereafter, we combined the climate and density-dependent IPM with forecasted climate variables to calculate the expected vital rates and population size year by year from 2022 to 2100. For each climate change scenario, we projected population dynamics in five steps: (i) we randomly selected one set of the forecasted climate variables from 21 global climate models, (ii) we sampled one set of vital rates, and the coefficients of each climate variable on the corresponding vital rate from 6000 sets of posteriors produced by the IPM; (iii) at each time we sampled from the previously estimated among-year variance of vital rates to incorporate the environmental stochasticity among years that is not due to the climate variable considered, and (iv) we ran the population model including demographic stochasticity and (v) repeated this procedure 100,000 times. For the scenario of no further climate change, the first step was different: we generated forecasted climate variables using Gaussian distributions with the mean and standard deviation of corresponding historical climate variables each time (i.e., 100,000 times).

Our model included density-dependence of vital rates by calculating breeding vacancy availability index year by year based on the projected population sizes in the relevant phases. For the offspring sex ratio \tau, we assume it was constant over time. Demographic stochasticity was incorporated by modelling population size as a random Poisson variable, which also ensures it is an integer and that extinction is clearly defined (N=0). Our projection assumes that climate impacts on vital rates in the projection period will be the same as in our study period (e.g., no non-additive effects of climate on vital rates), and there is no evolutionary rescue resulting in more adaptive vital rate responses to climate change.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Broad comments:\nThis study aims to determine the factors driving demographic fluctuations for an Australian songbird, leveraging an intensive\nlong-term field dataset. As the authors make the case themselves, I think this is a great illustration of the dangers of global\nchange for even abundant species. Studies such as this demonstrate how conservation and management efforts must take a\nbroad view if we are to preserve natural systems. While some of the findings in this study echo previous work from some of\nthe same authors, the IPM and forecasting (to my knowledge) appears to be novel for this dataset. The paper is generally\nwell written with findings mostly communicated clearly.\nHowever, I do think some additional considerations are needed in this study. Regarding the forecasting component of this\nwork, it would be helpful if more detail was provided regarding some of the decisions and the presentation of results. For\ninstance, why were these particular climate models used and why not use an ensemble approach? In forecasting efforts,\nuncertainty can come from many sources, including 1) uncertainty from different climate models (3 of which are considered\nseparately in this study), 2) uncertainty in climate model projections (i.e., different ‘runs’ of the same model since they are\nsensitive to initial conditions), 3) uncertainty in the coefficients estimated in the process model (IPM in this case), 4)\nuncertainty due to process error (variation in population dynamics that is unexplained in the IPM), and 5) uncertainty due to\nprevious states (since abundance in year t is dependent on abundance in year t-1). Uncertainty in forecasts can be\npartitioned into these different bins – see Zylstra and Zipkin 2021 – which I think would be greatly helpful to illustrate which\nof these elements matters most for these forecasts. I was uncertain which of the above were accounted for in these analyses.\nRegarding 1, I might suggest the authors collate data from a set of climate models, run the projections, aggregate the data\nfrom all of these to make a set of ensemble predictions, as has been done previously. This remove the awkwardness of\nincluding one model for the main text and places the others in the supplement. It also represents a more accurate accounting\nof how we think climate conditions will change into the future. I’d also advocate for integrate more climate models here to\ntake a larger sample of the CMIP6 projections. It doesn’t seem like 2 was accounted for in this case but I would advocate for\nthis being taken into account. I believe 3 and 5 are accounted for but was uncertain whether 4 was. The forecasts in all\ncases appear to be quite certain, without the typical cone that you typically see in forecasting studies, which makes me think\nthat some of the above are missing.\nAdditionally, while the forecasting component of this work stands as one of the highlights, I found it somewhat odd that there\nwas no discussion of the ecological forecasting literature – this field has advanced tremendously in recent years, with some\nwork that might be helpful for this study and would help place some of this into context (as some examples see Dietze et al.\n2017, 2024). In the methods, only a short section is dedicated to forecasting. I think some additional information covering the\nabove would be helpful.\nRelated to the environmental variables chosen, it’s clear from the results that these factors have a strong impact on fairy-\nwren demography. While there is some discussion about declines in other species and food availability, it was unclear to me\nwhy these variables impacted these birds so strongly. I think it might be helpful to dedicate some space to discussion the\n(proposed?) mechanistic link between these abiotic conditions and demography. Relatedly, I don’t see this as a major issue\nbut the variable selection process in choosing the environmental variables relevant here appears to be based on\nsignificance, which departs from the ways these are typically done. There are a number different model selection framework\nfor Bayesian model selection (Hooten and Hobbs 2015, Vehtari et al. 2017).\nDetailed comments:\nL 72: Might be helpful to have histogram (or density plots) or the population growth rates of each to support this statement (or\nrun a t-test).\nFig. 2: It might be helpful to describe what is meant by ‘natal’ here, so the reader does not have to necessarily locate that in\nthe text\nFig.4: Somewhat difficult to visualize the uncertainty, particularly for the forecasts. I would suggest not making the forecasts\ntransparent and separating the two regions in some other way (e.g., a line). Also, what exactly do the ribbons represent (i.e.,\nwhat there a threshold like 95% CI used)?\nFig. 4B: While ‘no further climate change’ is explained in the methods, some brief description would be helpful here.\nL 542: It was unclear to me how the ‘model includes environmental stochasticity due to climate’. I think some clarification\nwould be helpful for the reader\nReferences:\nDietze, M. C. 2017. Ecological forecasting. Princeton University Press.\nDietze, M., E. P. White, A. Abeyta, C. Boettiger, N. Bueno Watts, C. C. Carey, R. Chaplin-Kramer, R. E. Emanuel, S. K. M.\nErnest, R. J. Figueiredo, M. D. Gerst, L. R. Johnson, M. A. Kenney, J. S. McLachlan, I. Ch. Paschalidis, J. A. Peters, C. R.\nRollinson, J. Simonis, K. Sullivan-Wiley, R. Q. Thomas, G. M. Wardle, A. M. Willson, and J. Zwart. 2024. Near-term\necological forecasting for climate change action. Nature Climate Change 14:1236–1244.\nHooten, M. B., and N. T. Hobbs. 2015. A guide to Bayesian model selection for ecologists. Ecological Monographs 85:3–28.\nVehtari, A., A. Gelman, and J. Gabry. 2017. Practical Bayesian model evaluation using leave-one-out cross-validation and\nWAIC. Statistics and Computing 27:1413–1432.\nZylstra, E. R., and E. F. Zipkin. 2021. Accounting for sources of uncertainty when forecasting population responses to\nclimate change. Journal of Animal Ecology 90:558–561.", "reviewer_2": "This manuscript is based on fantastic data and sound statistical methods. The aim is to understand how weather impacts seasonal demographic rates and how this translates into temporal variation in population growth rates. This knowledge is\nthen used to predict how the population might evolve in the face of predicted climate change. The results show that relatively\nsmall changes in climate can lead to the rapid extinction of even common and widespread species. The manuscript\naddresses important issues. It not only provides insights into the pathway from weather via demographic rates to population\nchanges, but also uses this knowledge to more accurately predict the effect of future climatic changes.\nThe manuscript is generally well written. I have one important and some less important points to raise.\nMy main focus is on predicting population size. I was surprised by how precise population predictions are over such a long\ntimeframe (Fig. 4). In my experience, such predictions become highly uncertain even after just a few years. I identified the\nfollowing issues: Firstly, the weather variables used to conduct the population predictions are predictions themselves, yet\nthey are not subject to any uncertainty. Furthermore, the predictions are highly specific: under the no climate change\nscenario, for example, the predicted rainfall for the years 2071, 2072 and 2074 is 2.740160897, 2.029908585 and\n2.478329087 respectively. How can weather predictions for such a distant future be so specific? If uncertainties about future\nweather variable values are included, the population prediction would become more uncertain and the conclusion may\nchange. Therefore, I suggest that either uncertainty is included in the prediction or, if this is not possible because measures\nof prediction uncertainty are unavailable, this issue must be discussed and the conclusions must be made in light of this\nsimplification.\nSecondly, having examined the code, it seems to me that demographic stochasticity is not accounted for in the population\npredictions. This contrasts with the IPM analysis of the collected data. Omitting demographic stochasticity again reduces the\nuncertainty in the population predictions. I suggest either accounting for demographic stochasticity in the population\nprediction or explicitly stating and explaining why it is absent.\nMinor comments:\n1) L 26: ‚integrated projection model‘ does not seem to be correct, this ought to be ‘integrated population model’.\n2) L 63: Please provide the name of the city in addition. People from outside Australia are unlikely to know where the\nAustralian National Botanic Garden is located.\n3) L 103: you may consider providing a reference for the IPMs.\n4) L 111: is ‘foreign’ just ‘non-natal’? Or does ‘foreign’ mean a territory outside the study area? Please clarify.\n5) L 316: ‘long distances’. Could you quantify this? Important for the current study is to know whether such dispersal leads to\nemigration from the study area.\n6) L 403: I think the time indices for N and y should be 1, not 0. Otherwise the population model (eq. 1) needs to run from t=0.\nTechnically valid, but looks unusual.\n7) L 430: usually a Poisson distribution is used for modelling fecundity. Therefore, it would be informative to provide a\njustification of why the Binomial was used.\n8) L 444: move parameter delta closer to the ‘immigration rates’.\n9) L 456: priors are an integrative part of a Bayesian analysis. Therefore, I suggest that a table with the priors is provided\nsomewhere.\n10) L 456: This seems an odd reference for this statement. I don't think you need a reference here.\n11) L 473-482: For the modeling of density-dependence the number of occupied and vacant territories needs to be known.\nHow was the number of vacant territories inferred? I assume it is based on the assumption that territories remained stable\nstable\nover the years and so it can be inferred? If so, please also provide evidence that territories are stable over time.\n12) L 499: what exactly is the ‘vital rate’? The posterior mean of the vital rate? Was the estimation uncertainty accounted for\nwhile fitting the regression models?\n13) Extended Data Fig. 5: Here predictions of weather variables are shown. I do not understand, why the predictions can be\nso precise, e.g. in year 2080 it is value a and in year 2081 it is value b. Rather I would expect that the outcome of the\nprediction is a distribution, and that is also what is written in lines 378-380. So I think the dots represent a single realization\nfrom the distribution. Please clarify.\n14) Extended Data Fig. 6: The same comment as above applies here. Are the dots a single realization from predicted\ndistributions?"}, "author_response": "Dear Reviewers,\nWe have revised the manuscript according to your comments. Most importantly, we acknowledge\nthat our original projection models oversimplified matters, as they did not include the uncertainty\ngenerated by climate models and demographic stochasticity. This was the major criticism in your\ngeneral introductory comments and in some of your technical comments as well. We have therefore\ndiscussed the points you made with respect to uncertainty and stochasticity together. We are\nextremely grateful for your comments, as they have enabled us to improve the modelling\nsubstantially, and to increase the clarity of our arguments and results. The population forecasts now\nhave appropriate variances, but our central prediction of rapid population extinction remains. Please\nfind below our detailed responses to each of your comments, with your text indented. Any new cited\ntext is in italic font, and the line numbers in the new version of the manuscript are given.\nMajor comments of Reviewer #1\nThis study aims to determine the factors driving demographic fluctuations for an Australian\nsongbird, leveraging an intensive long-term field dataset. As the authors make the case\nthemselves, I think this is a great illustration of the dangers of global change for even\nabundant species. Studies such as this demonstrate how conservation and management efforts\nmust take a broad view if we are to preserve natural systems. While some of the findings in\nthis study echo previous work from some of the same authors, the IPM and forecasting (to\nmy knowledge) appears to be novel for this dataset. The paper is generally well written with\nfindings mostly communicated clearly.\nHowever, I do think some additional considerations are needed in this study. Regarding the\nforecasting component of this work, it would be helpful if more detail was provided regarding\nsome of the decisions and the presentation of results. For instance, why were these particular\nclimate models used and why not use an ensemble approach? In forecasting efforts,\nuncertainty can come from many sources, including 1) uncertainty from different climate\nmodels (3 of which are considered separately in this study), 2) uncertainty in climate model\nprojections (i.e., different ‘runs’ of the same model since they are sensitive to initial\nconditions), 3) uncertainty in the coefficients estimated in the process model (IPM in this\ncase), 4) uncertainty due to process error (variation in population dynamics that is\nunexplained in the IPM), and 5) uncertainty due to previous states (since abundance in year\nt is dependent on abundance in year t-1). Uncertainty in forecasts can be partitioned into these\ndifferent bins – see Zylstra and Zipkin 2021 – which I think would be greatly helpful to\nillustrate which of these elements matters most for these forecasts. I was uncertain which of\nthe above were accounted for in these analyses. Regarding 1, I might suggest the authors\ncollate data from a set of climate models, run the projections, aggregate the data from all of\nthese to make a set of ensemble predictions, as has been done previously. This remove the\nawkwardness of including one model for the main text and places the others in the\nsupplement. It also represents a more accurate accounting of how we think climate conditions\nwill change into the future. I’d also advocate for integrate more climate models here to take\na larger sample of the CMIP6 projections. It doesn’t seem like 2 was accounted for in this\ncase but I would advocate for this being taken into account. I believe 3 and 5 are accounted\nfor but was uncertain whether 4 was. The forecasts in all cases appear to be quite certain,\nwithout the typical cone that you typically see in forecasting studies, which makes me think\nthat some of the above are missing.\nMajor comments of Reviewer #2\nThis manuscript is based on fantastic data and sound statistical methods. The aim is to\nunderstand how weather impacts seasonal demographic rates and how this translates into\ntemporal variation in population growth rates. This knowledge is then used to predict how\nthe population might evolve in the face of predicted climate change. The results show that\nrelatively small changes in climate can lead to the rapid extinction of even common and\nwidespread species. The manuscript addresses important issues. It not only provides insights\ninto the pathway from weather via demographic rates to population changes, but also uses\nthis knowledge to more accurately predict the effect of future climatic changes.\nThe manuscript is generally well written. I have one important and some less important points\nto raise.\nMy main focus is on predicting population size. I was surprised by how precise population\npredictions are over such a long timeframe (Fig. 4). In my experience, such predictions\nbecome highly uncertain even after just a few years. I identified the following issues: Firstly,\nthe weather variables used to conduct the population predictions are predictions themselves,\nyet they are not subject to any uncertainty. Furthermore, the predictions are highly specific:\nunder the no climate change scenario, for example, the predicted rainfall for the years 2071,\n2072 and 2074 is 2.740160897, 2.029908585 and 2.478329087 respectively. How can\nweather predictions for such a distant future be so specific? If uncertainties about future\nweather variable values are included, the population prediction would become more\nuncertain and the conclusion may change. Therefore, I suggest that either uncertainty is\nincluded in the prediction or, if this is not possible because measures of prediction uncertainty\nare unavailable, this issue must be discussed and the conclusions must be made in light of\nthis simplification.\nSecondly, having examined the code, it seems to me that demographic stochasticity is not\naccounted for in the population predictions. This contrasts with the IPM analysis of the\ncollected data. Omitting demographic stochasticity again reduces the uncertainty in the\npopulation predictions. I suggest either accounting for demographic stochasticity in the\npopulation prediction or explicitly stating and explaining why it is absent.\nResponse to both reviewers\nWe agree that these issues needed to be addressed and we have now adopted our analytical\nframework to address the criticism by considering uncertainty and incorporating stochasticity.\nWe have now modified the population model to incorporate the “uncertainty from different climate\nmodels” (i.e., inter-model variability) by using your suggested ensemble approach. We were able to\nfind 21 CMIP6 climate models (instead of three) that provide monthly maximum and minimum\ntemperatures, and rainfall for all three greenhouse gas emission scenarios.\nInstead of modelling these separately as we had done in our previous version, for each step in the\nanalysis we now randomly selected one of the 21 models using the following approach – we describe\nthis as follows (Lines 577-586): “For each climate change scenario, we projected population\ndynamics in five steps: (i) we randomly selected one set of the forecasted climate variables from 21\nglobal climate models, (ii) we sampled one set of vital rates, and the coefficients of each climate\nvariable on the corresponding vital rate from 6,000 sets of posteriors produced by the IPM; (iii) at\neach time we sampled from the previously estimated among-year variance of vital rates to\nincorporate the environmental stochasticity among years that is not due to the climate variable\nconsidered, and (iv) we ran the population model including demographic stochasticity and (v)\nrepeated this procedure 100,000 times. For the scenario of no further climate change, the first step\nwas different: we generated forecasted climate variables using Gaussian distributions with the mean\nand standard deviation of corresponding historical climate variables each time (i.e., 100,000 times).”\nHowever, it was not possible to include the “uncertainty in climate model projections” (i.e., intra-\nmodel variability) as several CMIP6 global climate models only offer forecasted climate based on\na single initial condition. However, we note that inter-model variability is in general much larger\nthan intra-model variability and thus this would unlikely change the main prediction on the median\nextinction time much (though it may mean we have underestimated uncertainty somewhat). We now\nmentioned specifically (in condensed form) in Lines 381-383: “We did not consider intra-model\nuncertainty in our projections, as it was not available for all 21 models and is generally much\nsmaller than the variability among climate models49,50” that we were not able to include this source\nof uncertainty.\nWe have now incorporated demographic stochasticity also in the projection model by modelling\npopulation size in the next time step N(t+1) as a random Poisson variable, with the average rate\nbeing equal to N(t)*λ(t). This not only adds an important source of stochasticity when population\nsize is low (as is the case in studies of extinction), but also discretizes the population size variable\ninto a zero-positive integer, meaning that extinction is now clearly defined as N=0 (as opposed to\nan extinction threshold of N<1 when N was not an integer in our previous version).\nWe fully agree that population projections should incorporate all key sources of uncertainty, and we\nnow present this in our analysis. A complete decomposition of uncertainties, as suggested by\nReviewer 1, to identify which sources are most influential is certainly an interesting and valuable\ntopic. However, the focus of our manuscript is on the rapid population extinctions arising from the\naccumulation of many small climate impacts. A full uncertainty decomposition goes beyond the\nscope of this work, especially given that it is a complex issue—sources of uncertainty often interact\nin non-additive ways. We therefore feel it is best addressed in a separate manuscript, as including it\nhere would risk distracting from the central message.\nOther comments by Reviewer 1\nAdditionally, while the forecasting component of this work stands as one of the highlights, I\nfound it somewhat odd that there was no discussion of the ecological forecasting literature –\nthis field has advanced tremendously in recent years, with some work that might be helpful\nfor this study and would help place some of this into context (as some examples see Dietze\net al. 2017, 2024). In the methods, only a short section is dedicated to forecasting. I think\nsome additional information covering the above would be helpful.\nAgree, we have now discussed this point in Lines 288-296: “Our studied population of superb fairy-\nwrens is predicted to crash under all climate change scenarios. This result echoes the recent reports\nthat climate change will limit population growth and lead to local extinction of wild animals43,44.\nFor example, a recent study on African wild dogs (Lycaon pictus) found that climate warming will\nlead to population collapse in both intermediate and very high emission scenarios45. Although it is\nwell-established that climate change threatens the population persistence of wild animals8,46, our\nstudy uniquely provides demographic evidence that climate change poses an urgent threat to\ncommon species. As ecosystems are likely to be disproportionately impacted by the extinctions of\ncommon species5,6, more conservation attention to common species is needed in the future viability\nstudies of wild animals.”\nRelated to the environmental variables chosen, it’s clear from the results that these factors\nhave a strong impact on fairy-wren demography. While there is some discussion about\ndeclines in other species and food availability, it was unclear to me why these variables\nimpacted these birds so strongly. I think it might be helpful to dedicate some space to\ndiscussion the (proposed?) mechanistic link between these abiotic conditions and\ndemography. Relatedly, I don’t see this as a major issue but the variable selection process in\nchoosing the environmental variables relevant here appears to be based on significance,\nwhich departs from the ways these are typically done. There are a number different model\nselection framework for Bayesian model selection (Hooten and Hobbs 2015, Vehtari et al.\n2017).\nAgreed, we have now discussed this point in Lines 272-279: “The deleterious effects of climate\nchange, and particularly of increased winter temperature on our population of superb fairy-wrens\nis surprising, as a priori an increase in winter temperatures should move these endothermic birds\ncloser to their thermoneutral zone and lessen their energetic burden. Although we lack the relevant\ndata, which would be extremely expensive to obtain, we suspect the problem posed by sudden warm\nspells followed by cold spells10 may negatively affect the ectothermic arthropod prey exploited by\nthe birds, possibly by provoking emergence and activity at the wrong time of year.”\nFor the variable selection, we would like to highlight that all climate/density variables were\nidentified from our previous studies (papers; see the “Climate data” and “Density dependence”\nsections in the Methods). Here, we used the results from those studies to examine the effects of\nclimate/density variables on vital rates. Additionally, we followed the method used by Woodworth\net al. (2017) to reduce the number of variables in the final retrospective path analysis model. This\nmethod uses a 90% instead of 95% credible interval for variable selection because this reduces the\nchances of overlooking climate/density effects (false negatives). Finally, we validated the predictive\nvalue of our model and the results are pretty good (Fig. 4A). It indicated that the final set of\nclimate/density variables produced by our variable selection could project the population size\nreliably.\nWoodworth, B. K., Wheelwright, N. T., Newman, A. E., Schaub, M. & Norris, D. R. Winter\ntemperatures limit population growth rate of a migratory songbird. Nat Commun 8, 1–9 (2017).\nL 72: Might be helpful to have histogram (or density plots) or the population growth rates of\neach to support this statement (or run a t-test).\nAgreed. we have now deleted this point and rewrote as: “Male numbers were consistently greater\nthan female numbers (Fig. 1A).” (Line 72).\nFig. 2: It might be helpful to describe what is meant by ‘natal’ here, so the reader does not\nhave to necessarily locate that in the text\nAgreed, we have stated this explicitly in Lines 95-97: “L: natal juvenile females (who remain in the\nterritory where they were born for their first winter); O: non-natal juvenile females (who have\nimmigrated to a subordinate role before their first winter)”.\nFig.4: Somewhat difficult to visualize the uncertainty, particularly for the forecasts. I would\nsuggest not making the forecasts transparent and separating the two regions in some other\nway (e.g., a line). Also, what exactly do the ribbons represent (i.e., what there a threshold like\n95% CI used)?\nAgreed, we have now used the 10/90% percentiles of population size (i.e., the population prediction\nintervals) and as result the uncertainty of population projection is clear in the new figure.\nFig. 4B: While ‘no further climate change’ is explained in the methods, some brief description\nwould be helpful here.\nDone, we have now explained this in Lines 182-183: “the scenario of no further climate change (B;\ni.e., climate will not further change due to emissions)”.\nL 542: It was unclear to me how the ‘model includes environmental stochasticity due to\nclimate’. I think some clarification would be helpful for the reader\nDone, we have revised it in Lines 578-582: “(ii) we sampled one set of vital rates, and the\ncoefficients of each climate variable on the corresponding vital rate from 6,000 sets of posteriors\nproduced by the IPM; (iii) at each time we sampled from the previously estimated among-year\nvariance of vital rates to incorporate the environmental stochasticity among years that is not due to\nthe climate variable considered”.\nReferences:\nDietze, M. C. 2017. Ecological forecasting. Princeton University Press.\nDietze, M., E. P. White, A. Abeyta, C. Boettiger, N. Bueno Watts, C. C. Carey, R. Chaplin-\nKramer, R. E. Emanuel, S. K. M. Ernest, R. J. Figueiredo, M. D. Gerst, L. R. Johnson, M. A.\nKenney, J. S. McLachlan, I. Ch. Paschalidis, J. A. Peters, C. R. Rollinson, J. Simonis, K.\nSullivan-Wiley, R. Q. Thomas, G. M. Wardle, A. M. Willson, and J. Zwart. 2024. Near-term\necological forecasting for climate change action. Nature Climate Change 14:1236–1244.\nHooten, M. B., and N. T. Hobbs. 2015. A guide to Bayesian model selection for ecologists.\nEcological Monographs 85:3–28.\nVehtari, A., A. Gelman, and J. Gabry. 2017. Practical Bayesian model evaluation using leave-\none-out cross-validation and WAIC. Statistics and Computing 27:1413–1432.\nZylstra, E. R., and E. F. Zipkin. 2021. Accounting for sources of uncertainty when forecasting\npopulation responses to climate change. Journal of Animal Ecology 90:558–561.\nReviewer #1 (Remarks on code availability):\nI have briefly looked over the code to ensure the major components were present but not in\nany great detail.\nThank you for your suggested papers and efforts.\nOther comments by Reviewer 2\n1) L 26: ‚integrated projection model‘ does not seem to be correct, this ought to be ‘integrated\npopulation model’.\n2) L 63: Please provide the name of the city in addition. People from outside Australia are\nunlikely to know where the Australian National Botanic Garden is located.\n3) L 103: you may consider providing a reference for the IPMs.\nAll corrected or provided.\n4) L 111: is ‘foreign’ just ‘non-natal’? Or does ‘foreign’ mean a territory outside the study\narea? Please clarify.\nDone, we revised it to: “non-natal juveniles that moved out of natal territory during the recruitment\nphase” in Lines 111-112.\n5) L 316: ‘long distances’. Could you quantify this? Important for the current study is to know\nwhether such dispersal leads to emigration from the study area.\nDone, we now clarified it in Lines 337-339: “female offspring always disperse away from their natal\nterritory in the first year of their life, often moving more than two territories (few over 5 km) 17,18.”\nYes, many juvenile females emigrated from our study area since it covers ~60 hectares.\n6) L 403: I think the time indices for N and y should be 1, not 0. Otherwise the population\nmodel (eq. 1) needs to run from t=0. Technically valid, but looks unusual.\nAgreed, we have corrected this. Thank you for pointing out.\n7) L 430: usually a Poisson distribution is used for modelling fecundity. Therefore, it would\nbe informative to provide a justification of why the Binomial was used.\nWe agree that a Poisson distribution is usually used for modeling fecundity. Indeed, we modelled\nthe true number of reproduced individuals with Poisson distributions (equations 1 & 2). We only\nused a binomial distribution for the age ratio data which is likely to contain observation errors\n(equation 6), following Zimmerman et al. 2010. This approach mainly has two advantages. (1) It\nallows the number of adult females 𝑦 to be 0, although it was not encountered in our study. (2) It\nprovides more precise estimates of fecundity. Also note that our survey was very thorough, so there\nis a good reason to believe detection probability is not different between adult and juvenile\nindividuals (except for the dispersing juvenile females). Therefore, adult and juvenile detection\ncancels each other in equation 6.\nZimmerman, G. S., W. A. Link, M. J. Conroy, J. R. Sauer, K. D. Richkus, and G. S. Boomer. 2010.\nEstimating migratory game-bird productivity by integrating age ratio and banding data. Wildlife\nResearch 37:612–622.\n8) L 444: move parameter delta closer to the ‘immigration rates’.\nWe have two immigration rates that were estimated in two different phases. Therefore, for each\nparameter delta, we prefer to put it right after the relevant phase.\n9) L 456: priors are an integrative part of a Bayesian analysis. Therefore, I suggest that a table\nwith the priors is provided somewhere.\nDone, we have now added it in Lines 482-484: “We used vague priors Normal (0, 0.01) for all mean\nparameters, Gamma (0.01, 0.01) for all precision (inverse of variance) parameters, and Uniform\n(0, 1) for all detection probability parameters.”\n10) L 456: This seems an odd reference for this statement. I don't think you need a reference\nhere.\nDone, we have now deleted the reference.\n11) L 473-482: For the modeling of density-dependence the number of occupied and vacant\nterritories needs to be known. How was the number of vacant territories inferred? I assume\nit is based on the assumption that territories remained stable over the years and so it can be\ninferred? If so, please also provide evidence that territories are stable over time.\n“Superb fairy-wrens form socially monogamous pairs and live in year-round territories”. (Line 324)\nIn this paragraph, we have noted that “the number of vacancies is associated with the number of\nadult females that died in the previous recruitment and non-breeding phase” in Lines 506-507. So,\nthe number of vacancies were calculated right before the scramble phase in each year. It seems that\nthere is no need to assume territories remained stable over the years, but if there is no fission or\nfusion, the territories are extremely stable.\n12) L 499: what exactly is the ‘vital rate’? The posterior mean of the vital rate? Was the\nestimation uncertainty accounted for while fitting the regression models?\nDone, we have now rewritten it and to make it clear in Lines 527-529: “We fitted a series of models,\nin which each vital rate (fecundity, sex- and phase-specific survival rates, and female phase-specific\nimmigrations rates, see Extended Data Fig. 1) was fitted as a logistic or log-linear function of a\nsingle climate variable” and Lines 537-540: “We modelled each estimated vital rate in the IPM with\na mean 𝛼 and a process error 𝜀 , where 𝜀 followed a Gaussian distribution with mean 0 and\nstandard deviation 𝜎. The process error estimates the among-year variance of a given vital rate\nthat is not due to any modelled effects of climate/density variables (i.e., the fixed effects).” Therefore,\nthe estimation uncertainty was also included in the IPM (equations 12 & 13).\n13) Extended Data Fig. 5: Here predictions of weather variables are shown. I do not\nunderstand, why the predictions can be so precise, e.g. in year 2080 it is value a and in year\n2081 it is value b. Rather I would expect that the outcome of the prediction is a distribution,\nand that is also what is written in lines 378-380. So I think the dots represent a single\nrealization from the distribution. Please clarify.\nYes, the previous Extended Data Fig. 5 was a single realization from the distribution. To keep\nconsistent with population projection under climate change scenarios, we have now changed this to\n100,000 realizations in Lines 584-586: “we generated forecasted climate variables using Gaussian\ndistributions with the mean and standard deviation of corresponding historical climate variables\neach time (i.e., 100,000 times).” The temporal change of climate variables was no longer provided\nin the revised Supplementary Information since there are many realizations.\n14) Extended Data Fig. 6: The same comment as above applies here. Are the dots a single\nrealization from predicted distributions?\nWe obtained monthly climate data from global climate models of Coupled Model Intercomparison\nProject Phase 6 (CMIP6) and then calculated forecasted climate variables (e.g., spring rainfall from\nAugust to November). As there are 21 global climate change models, the temporal change of\nforecasted climate variables is now summarized across all 21 projections in Supplementary Table\nS1.\nReviewer #2 (Remarks on code availability):\nThe code is perfectly commented and structured. When I run the code, all the analyses run\nwithout any problems.\nThank you."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Broad comments:\nThis is a revision. I also reviewed the previous version of this manuscript. I appreciate the authors' efforts to address\nprevious concerns brought up by myself and the other reviewer.\nThe authors have done a nice job accounting for the additional sources of uncertainty, integrating forecasts from multiple\nclimate models, and clearly outlining what was done in the forecasting portion of the study. These were the largest\nshortcomings in the previous version of the manuscript, in my opinion. I believe the previous comments from both myself\nand the other reviewer are now adequately addressed. I’d like to commend the authors on this nice study.\nI have only a few other comments. Firstly, regarding Table S1 it is useful that the climate model names are listed in the note\nand I think the estimated trend over time is useful. However, 1) I don’t think it’s necessary to include the number of significant\ntrends in each case. There is a multiple comparisons problem (inflated type I error) when making so many comparisons.\nOr is this a hierarchical model with varying slopes, varying intercepts? 2) It’s stated that the mean and 95% CI are given in the\nnote. However, from what model? Was this is hierarchical model and this was the group-level response? Finally, 3) it was\nunclear to me why the global extend was used for changes in annual mean temp but Canberra for all other variables.\nShouldn’t we only be interested in the study area? These were the data that were used for the projections, themselves,\ncorrect?\nL386: I agree this is the best way to approach this. It might be worth adding a citation noting that ensemble approaches tend\nto have better forecast performance."}, "author_response": "Dear Reviewer,\nWe have revised the manuscript according to your comments. Please find below our detailed\nresponses to each of your comments, with your text indented. Any new cited text is in italic font,\nand the line numbers in the new version of the manuscript are given.\nReviewer #1 (Remarks to the Author):\nBroad comments:\nThis is a revision. I also reviewed the previous version of this manuscript. I appreciate the\nauthors' efforts to address previous concerns brought up by myself and the other reviewer.\nThe authors have done a nice job accounting for the additional sources of uncertainty,\nintegrating forecasts from multiple climate models, and clearly outlining what was done in\nthe forecasting portion of the study. These were the largest shortcomings in the previous\nversion of the manuscript, in my opinion. I believe the previous comments from both myself\nand the other reviewer are now adequately addressed. I’d like to commend the authors on this\nnice study.\nI have only a few other comments. Firstly, regarding Table S1 it is useful that the climate\nmodel names are listed in the note and I think the estimated trend over time is useful. However,\n1) I don’t think it’s necessary to include the number of significant trends in each case. There\nis a multiple comparisons problem (inflated type I error) when making so many comparisons.\nOr is this a hierarchical model with varying slopes, varying intercepts? 2) It’s stated that the\nmean and 95% CI are given in the note. However, from what model? Was this is hierarchical\nmodel and this was the group-level response?\nThank you for your positive comments. Agree, we have now removed the number of significant\ntrends from Supplementary Table S1. We also clarified that we used simple linear regressions in\nSupplementary Information, Line 27-28:\n“The means were calculated from the slopes of simple linear regressions for 21 climate models, with\n95% confidence intervals shown in parentheses.”\nFinally, 3) it was unclear to me why the global extend was used for changes in annual mean\ntemp but Canberra for all other variables. Shouldn’t we only be interested in the study area?\nThese were the data that were used for the projections, themselves, correct?\nAgreed, we have now added Δ annual mean temperature for Canberra in Supplementary Table S1.\nL386: I agree this is the best way to approach this. It might be worth adding a citation noting\nthat ensemble approaches tend to have better forecast performance.\nAdded in Line 352-353:\n“As ensemble approaches tend to have better forecast performance51, we therefore focus on\nensemble projections.”\nThe reference is “Tebaldi, C. & Knutti, R. The use of the multi-model ensemble in probabilistic\nclimate projections. Philosophical Transactions of the Royal Society A: Mathematical, Physical and\nEngineering Sciences 365, 2053–2075 (2007).”"}]

Decision letter