More frequent extreme climate events stabilize reindeer population dynamics

Extreme climate events often cause population crashes but are difficult to account for in population-dynamic studies. Especially in long-lived animals, density dependence and demography may induce lagged impacts of perturbations on population growth. In Arctic ungulates, extreme rain-on-snow and ice-locked pastures have led to severe population crashes, indicating that increasingly frequent rain-on-snow events could destabilize populations. Here, using empirically parameterized, stochastic population models for High-Arctic wild reindeer, we show that more frequent rain-on-snow events actually reduce extinction risk and stabilize population dynamics due to interactions with age structure and density dependence. Extreme rain-on-snow events mainly suppress vital rates of vulnerable ages at high population densities, resulting in a crash and a new population state with resilient ages and reduced population sensitivity to subsequent icy winters. Thus, observed responses to single extreme events are poor predictors of population dynamics and persistence because internal density-dependent feedbacks act as a buffer against more frequent events.

xtreme climate events can induce severe population crashes and destabilized population dynamics across animal taxa and biomes [1][2][3][4] . For instance, El Niño years severely affect the dynamics of birds 5 , and droughts have led to extinctions in butterflies 6 . As accumulated evidence now suggests that global warming comes with an increase in the frequency of extreme climate events [7][8][9] , their ecological impacts are also receiving more attention 4 , yet mostly in plants. Given the inherent rareness of extreme events and the anecdotal approaches to study them in animals, the scientific focus has almost exclusively been on single events and their short-term effects 10 . However, especially in longlived species where intrinsic properties regulate population dynamics, the impact of an environmental perturbation could vary with how often it occurs 11 . Thus, ignoring longer-term impacts, as well as anticipated changes in event frequencies under global warming 9 , may ultimately lead to biased predictions of population persistence 4 .
As some types of (previously extreme) events become progressively less rare-typically in climate change hotspots, such as the Arctic 12 -novel opportunities for mechanistic insights and a predictive understanding of their ecological impacts are now emerging. In the Arctic, extreme warm spells and rain-on-snow (ROS) events in winter may cause impenetrable snow-packs and even encapsulate the entire vegetation in thick ground-ice 13 . Such environmental perturbations on the tundra are no longer that rare [12][13][14][15] , and population crashes and destabilized dynamics linked to icing events have been reported for a range of herbivore species 2 , including muskoxen Ovibos moschatus 14,16 , caribou, and wild reindeer Rangifer tarandus [17][18][19] . One intuitive, yet perhaps naive, extrapolation from observations of single population crashes following extreme ROS events would be to expect more variable dynamics and greater extinction risk with continued warming. However, theoretical studies indicate that this may not necessarily be the case 11 . The effects of environmental stochasticity on population growth can be buffered by density-dependent feedbacks [20][21][22][23] . Perturbations may have little or no impact at low density when resource competition is weak 24 . In addition, population responses to environmental stochasticity can depend on demographic structure, with some age classes being less sensitive to environmental fluctuations than others 20,25 . A change in population structure towards more resilient age classes after a population crash, i.e. a new population state 26 , may therefore promote positive population growth rates and reduce the probability of new crashes in subsequent years. Thus, it has been suggested that high frequencies of bad years may lead to less variable population growth and, hence, stabilized rather than destabilized population dynamics 11 . However, empirical support for this prediction is still lacking.
Here, based on demographic population modelling of empirical time-series data 27 , we evaluate how changes in the frequency of rainy and icy winters affect wild reindeer R. t. platyrhynchus population dynamics in Svalbard, a climate change hotspot in the High Arctic 13 . Because of the rapidly warming winter climate and the strong ROS signals in both reindeer demographic performance 28 and abundance 29 , this northernmost ungulate represents an excellent case study for exploring the effects of more frequent extreme events. We show that the impact of an extreme ROS and icing event on reindeer survival, fecundity, and population growth rate is strongly ageand density-dependent. A population crash causes relaxation of density dependence, more resilient age structure, and, thereby, a long-lasting reduction in the population sensitivity to subsequent extreme events. Thus, because effects of environmental stochasticity are modified by internal density-dependent feedback, frequent extreme events dampen the population dynamics and even reduce the extinction risk.

Results and discussion
Exploring the climate-density interaction. As a preliminary analysis, we first explored the impact of ROS and population density on annual reindeer population growth rates over the study period 1994-2014 ( Fig. 1), obtained from the posterior means of an integrated population model (IPM) combining mark-recapture and count data 27,30 (Fig. 2). Because the trend for Arctic greening 31 due to gradually warmer and longer summers 32 is likely to influence the carrying capacity of the reindeer population 28 , we accounted for variation in winter length and a linear change in carrying capacity. As expected, we found a strong negative effect of ROS on annual population growth rate (Fig. 1c Fig. 1 Density-dependent effects of rain-on-snow (ROS) and icing shape reindeer population fluctuations. a Occasionally, icing causes die-offs. This calf died from starvation in the 2001/2002 winter, which was characterized by high reindeer density and high ROS amounts. b Annual fluctuations in ROS (blue line) and total female population size N (black line) in the study period. c Density-dependent effect of ROS on the population growth rates (see Supplementary Table 1). For illustration, population densities (N) are classified into low (yellow symbols), medium (orange), and high (red) and Supplementary Table 1). This effect was diminished at low densities, when food competition is weak or negligible (cf. ref. 33 ) even under icy conditions. Accounting for age-specific effects of climate and density. Second, based on these preliminary findings, we explored the underlying demographic mechanisms by modelling annual agespecific survival and fecundity rates, obtained for 9090 posterior samples from the IPM 27 , as a function of weather and population size, allowing the effect of ROS to depend not only on density but also on age class 20 (Supplementary Figs. 1 and 2, Supplementary  Table 2). As expected in long-lived ungulates 22,34 , vital rates of young and old age classes were the most variable and the most strongly influenced by ROS events (Supplementary Table 2). Third, using these functions for survival and fecundity rates, together with past age structures ( Supplementary Fig. 3) and weather conditions, we reconstructed fluctuations in vital rates and population sizes. Predicted population growth rates (Supplementary Fig. 4) were strongly correlated with observed growth rates (Pearson's correlation r = 0.89), suggesting high predictive power of the population model ( Supplementary Fig. 5).
Population dynamics and persistence under climate change. We evaluated the demographic consequences of a set of ROS scenarios for the reindeer population, using stochastic simulations of the population model (Fig. 3). Global circulation models suggest that mid-winter warm spells and heavy ROS events will become more frequent in the future Arctic 12,14,15 , including Svalbard. Therefore, we varied the frequency of extreme ROS winters (see Supplementary Fig. 6) from very low (virtually never) to low, medium (as observed in the past, i.e. 1962-2014), high, and very high (the likely future scenario). A very high frequency of extreme ROS winters (rightmost panels in Fig. 3) reduces the mean population size by only 11% (Fig. 4a, Table 1) compared to the scenario describing observed historical conditions (i.e. medium ROS winter frequency, mid panels in Fig. 3) and by 25% (Fig. 4a, Table 1) compared to the very low frequency scenario (leftmost panels in Fig. 3). However, very frequent ROS winters also result in a strong reduction in the temporal variability in population sizes and growth rates (Figs. 3 and 4, Table 1). Such a change in population variability occurs because the impact of a bad year interacts with density and age structure (Fig. 5), which, in turn, are likely to be influenced by the time since the previous bad year. This also has implications for the frequency and magnitude of population crashes under different ROS scenarios (Fig. 3). Accordingly, the probability of going extinct (population size N = 0) during a period of 100 years is about 15,000 times higher for the medium ROS scenario, i.e. the observed historical climate, than for the very high ROS scenario anticipated under continued global warming (Table 1). Likewise, the probability of going quasi-extinct (here arbitrarily defined as N < 100) is 10-fold.
Step 1: Preliminary analysis of density and climate effects on population growth rates Linear regression of log(N t+1 /N t ) against density and climate covariates (Fig. 1c) Data Time series of:

Summer temperature
Winter length

ROS
Step 2: Survival and fecundity as functions of density and climate Models with density (N posthunt,t ) and climate covariates fitted to each posterior time series of S i,t and F i,t Step 3: Model checking Step by step reconstruction of past S i,t , F i,t and N prehunt,i,t using model from step 2 Step 4a: ROS scenarios Five ROS scenarios simulated using inverse transformation method Step 4b: Population projection with different ROS scenarios Boxes with rounded corners show the analytical steps used in this paper. Square boxes contain data or estimates that are fed into the models in steps 1-4. S i,t and F i,t refer to survival and fecundity estimates for age class i at time t, N t is the population size at time t, either age-structured before the hunting season (N prehunt,i,t ) or right after the end of the hunting season (N posthunt,t ), and ROS is a measure of rain-on-snow, as defined in the Methods. Note that in this figure, we use the term density in place of population size when considering density effects of population size on population growth (i.e., density dependence). This is to aid the reader in following how N t was used in different steps of our analysis (as a cause of density dependence, or as a response variable), and does not indicate a change in the way N t was estimated  such differences in long-term population dynamics and persistence under contrasting ROS scenarios, we ran deterministic population simulations where we increased the time elapsed between two extreme ROS winters (Fig. 6). We modelled the extreme ROS winter as the maximum ROS amount ever recorded at the local weather station (i.e. 1996, Fig. 1b). Provided a rather high initial density characterized by a low proportion of primeaged animals (i.e. 3-8-year olds) in year t = 0 (i.e. representative of the population state before a known population crash in 1996), an extreme ROS winter leads to a dramatic decline in population size in year t = 1 (Fig. 6a). In particular, many calves and old animals die ( Supplementary Fig. 1) and very few new calves are born ( Supplementary Fig. 2) in crash years 27,28 , also resulting in a marked increase in the proportion of prime-aged animals ( Supplementary Fig. 3). With no subsequent ROS winter (Fig. 6a, leftmost panel), the low density (i.e. relaxation of density dependence) and new age structure allow the population to recover towards an asymptotic population size (Fig. 6a) and stable age structure ( Supplementary Fig. 7, leftmost panel). If a second ROS winter occurs immediately after the first one (Fig. 6a, second leftmost panel), the model predicts no additional decrease in population size before recovery commences. This occurs because the previous year's crash generated a new population state, with low density as well as large proportion of prime-aged animals, showing little sensitivity to ROS. This promotes population recovery (Fig. 5). If the second ROS winter is delayed until t = 2, the previous year's recovery in population size induces a slight population decline due to the harsh feeding conditions. Because of these lagged effects of a population crash, it takes about 7 years before the impact of an icy winter returns to the initial level at t = 0 (Fig. 6b). Accordingly, if there are several extreme ROS winters in a row, population size converges towards a reduced density of 1000-1500 individuals due to relaxation of density dependence (Fig. 7). These interactions between climate effects, density, and age structure explain why our study population did  Low density, high % prime-aged Low density, low % prime-aged High density, high % prime-aged High density, low % prime-aged -0.8 Expected growth rate not crash in recent icy winters (Fig. 1b). Such lagged responses to environmental stochasticity are predicted by theory 11,22 , yet poorly documented (but see 20 ). Thus, to our knowledge, this study provides the first evidence from the wild that fluctuations in population size of long-lived species are dampened when extreme climate events become more frequent 11 .
In contrast to our study area (a hotspot for winter climate change 13,14 ), heavy ROS is still an uncommon phenomenon in many parts of the Arctic 14 , where the anticipated near-future winter warming 12,15 may first imply moving from very low to low frequency of extreme ROS winters (Fig. 3). When ROS is still that rare, ungulate populations are more likely to be in the kind of state (in terms of density, as well as demographic structure) that may lead to a crash when an extreme event occurs. Recent observations of occasional population crashes across the circumpolar Arctic 14,17,18 support this, with potentially large socioeconomic and ecosystem implications 35 . However, one important message from our simulations, supported by observed population dynamics (Fig. 1b), is that the impact of an extreme ROS winter on the population growth rate is far less dramatic if it occurs soon after a previous perturbation (i.e. moving towards high and very high ROS frequencies, Fig. 3). This is a likely characteristic of the near-future winter climate in many coastal Arctic regions 14,15 , including Svalbard.
Extreme environmental perturbations can have unexpected ecological impacts 36 , and therefore represent one of the major challenges in future predictions of ecosystem change 1,37 . Because extreme climate events are by definition rare, studies based on empirical data from the wild are also typically of anecdotal nature 10 . Our case study clearly demonstrates that populationdynamic inference based solely on single events and their shortterm impacts-ignoring potential long-term impacts, as well as the consequences of multiple events-may lead to erroneous conclusions. In particular, our results emphasize how internal density-dependent feedback processes can modify the effects of environmental stochasticity and, hence, buffer populations of long-lived species if extreme events become the norm due to global warming.

Methods
Study area and species. Our study population of wild Svalbard reindeer is located in the Reindalen, Semmeldalen, and Colesdalen valley system in central Spitsbergen (78°N, 15°E), Svalbard, Norway. The area is characterized by U-shaped coastal valleys with High Arctic tundra vegetation of low stature, dominated by mosses, graminoids, dwarf shrubs, and forbs. Because of semi-isolation by the sea, steep mountains, and glaciers, as well as a stationary behaviour, there is very little exchange of animals with other nearby reindeer populations 38 . The reindeer occur alone or in small groups and are not subject to significant levels of predation (polar bear attacks are very rare 39 ), inter-specific competition, or insect harassment. Each fall, there is a very low level of harvest 40 , and some reindeer have been culled for scientific purpose 41 . Both hunted and culled animals are reported and their age is determined 42 .
ROS has previously been reported as the main climatic driver of vital rates and population growth in Svalbard reindeer 2,19,29,40 . The proximate mechanism behind this is that ROS causes icing, which restricts food availability, in turn affecting body mass loss during winter 28 . During the study period, summers have become warmer causing increased plant productivity 2,32,43,44 , representing the most plausible explanation for positive trends in autumn body mass and population size (Fig. 1b) in our study population 28 .
Data. The reindeer data used in this study originate from a mark-recapture study, which protocol applies with and is approved by the Norwegian Animal Research Authorities and the Governor of Svalbard. The reindeer are captured as calves (and recaptured in later years) in a net between two snow mobiles during April each year, and a post-breeding resighting survey is performed in early August 40 . A posterior sample of 9090 estimates of annual survival (Supplementary Fig. 1) and fecundity ( Supplementary Fig. 2), as well as population size (N), for six female age classes during 1994-2014, were obtained from an IPM 27,30,45 (Fig. 2). The posterior means of each annual demographic parameter are hereafter referred to as observed values or observations. The six age classes considered are 0 (calves), 1, 2, 3-8, 9-11, and ≥12 year olds, corresponding to age classes 1-6 in Supplementary Figs. 1 and 2, Supplementary Table 2. The model combines individual mark-recapture data with population counts within a Bayesian state-space modelling framework, accounting for observation error and demographic stochasticity. The small and controlled level of harvest and scientific culling was also accounted for. Detailed modelling and data methodology are described in ref. 27 , with further model updates in ref. 30 .
We obtained daily historical weather data for Longyearbyen airport, ca 25 km from our study population, from the Norwegian Meteorological Institute (freely available at http://eklima.met.no). Based on daily mean temperatures and precipitation during 1962-2014, we first calculated annual winter rain amount by summing the amount of precipitation recorded at temperatures ≥1°C during November-April 2,40 . We added 1 to this value and log-transformed it to get the variable ROS, an icing/winter harshness proxy 29 (Fig. 1b). Second, we estimated the length of the winter based on a 10-day running mean of daily temperatures. We considered the onset of winter as the first day in autumn when the running mean was <0°C and then stayed <0°C for a minimum of 10 consecutive days. Likewise, end of winter was estimated as the first day (after the winter period November-April) when the running mean was >0°C and then stayed >0°C for a minimum of 10 consecutive days.
Population growth rate as a function of density and climate. Based on the observed total population sizes (N t , where t is time) from the IPM 27,30 , we first explored climate-density interactions by fitting a linear regression of annual population growth rate (log(N t+1 /N t ) against log(N t , linearly detrended) and climate covariates, namely ROS and length of the winter (Fig. 2). We found a negative effect of length of the winter (t to t + 1) and a negative interaction effect between the amount of ROS (t to t + 1) and population size (Supplementary Table 1), suggesting a density-dependent ROS effect and supporting the assumption that increased plant abundance due to warmer summers 32,44 has a gradual effect on the carrying capacity K (thereby, the detrended N). To visually illustrate the ROS-density interaction effect (Fig. 1c), we divided population sizes into high, medium, and low N (detrended). The respective linear regression model replacing the covariate N with these population size classes (as factor) provided qualitatively similar results as the model described above (Supplementary Table 1).
Survival and fecundity as functions of density and climate. We obtained 9090 posterior samples from the IPM 27,30 , each of them consisting of annual age-classspecific demographic rates (i.e. survival and fecundity) and population sizes from 1994 to 2014. For each of these posterior samples, we investigated the effects of weather and population density on survival and fecundity, generating in total 9090 population models. First, the effects of weather and population density on survival S of each age class i at year t were estimated separately for each posterior sample (Fig. 2, Supplementary Fig. 1), using linear mixed-effects models (function lme in R package nlme). As survival might be negatively affected by an increase in population size N posthunt (i.e. number of calves + adult females at the end of the hunting season), we tested for an effect of (scaled) population size in year t on survival. We fitted a density-dependent ROS effect 20 (cf. Fig. 1c) using the form ROS ′ t ¼ ROS t e k N posthunt;t (see ref. 21 for a similar approach). This form of the interaction effect ensures that the effect of ROS is strictly negative or positive (depending on the fitted coefficient of ROS') for all values of N posthunt (Supplementary Figs. 8 and  9). In contrast, the simpler and more common specification of an interaction effect, such as k ROS t N posthunt;t , would implicitly lead to a switch in sign of the effect of ROS at some level of N, and the stronger the interaction effect, the more likely it is that the change of sign is well within data range. Note that k was estimated using an optimization function aiming at minimizing Akaike's Information Criteria. Moreover, mean survival is likely to differ among age classes, and the effect of climate may also differ among age classes (e.g. ref. 20 ), thus we included an interaction between age and the density-dependent ROS effect. Finally, year was added as a fixed (numeric) effect to account for trends. Year was also added as a random (intercept) effect to account for non-independence among age classes, leading to a survival model with the following form: where a i is the age class i, N posthunt,t is population size just at the end of the hunting season, and ε S is the residuals of the regression. The posterior distributions of the parameters of Eq. (1) shown in Supplementary Table 2 represent a combination of uncertainty due to finite time series, uncertainty stemming from stochasticity, and uncertainty in the IPM. Likewise, we ran a model of similar structure for fecundity F of each age class i at year t, as females produce at most one calf per year ( Supplementary Fig. 2). Note that there is one age class less than in the survival analysis, as calves do not get pregnant: Again, the posterior distributions of the parameters of Eq. 2 shown in Supplementary Table 2 represent a combination of uncertainties (see above for survival). Six and five age classes were considered for survival and fecundity, respectively 27,30 . We assumed equal survival and fecundity among ages 3-8 years, among ages 9-11 years, and among ages >11 years.
Reconstructing past vital rates and population growth rates. Based on the above models of vital rates, for each posterior sample we estimated age-specific survival and fecundity rates for past conditions of population density, winter harshness (i.e. ROS amount), and winter length (Fig. 2, Supplementary Table 2). Importantly, to account for sources of environmental stochasticity due to processes other than covariates included in the model, we estimated a covariance matrix Σ of the different vital rates (fecundity, survival) for all age classes based on the random year effects and the residuals ε S i;t and ε F i;t of Eqs. 1 and 2. From this covariance matrix Σ, we generated 100 new residuals from a multivariate normal distribution.
These rates then allowed the population size at time t + 1 to be estimated from the population size of each age at time t. The population size just before the hunting season at time t + 1, N sim,t+1 , corresponds to the sum of females of different ages: Each of these terms can be estimated based on the vital rates. The number of calves produced, N sim calves,t+1 (first term on right side of Eq. 3), consists of calves produced by females of each age (except yearlings). N sim calves,t+1 was modelled using a binomial process to allow for demographic stochasticity (i.e. chance events that affect individuals independently): where F sim i,t is the estimated fecundity rate for females of age i at time t, resulting in calves at time t + 1. For instance, F sim 2,t is the probability of a 2-year-old female at time t having a calf at heel at time t + 1 (depending on climate and population size at time t, see Eq. 2). The number of female calves was then drawn from a binomial distribution with a probability 0.5 (i.e. we assumed a balanced sex ratio). N sim yearlings,t+1 (second term on right side of Eq. 3) corresponds to the number of female calves that have survived from time t to time t + 1 and was also modelled using a binomial process to include demographic stochasticity. Moreover, a few female calves are removed from the population by hunting (H) or scientific culling (C). Thus the number of female yearlings is modelled as follows: N sim yearling;tþ1 $ Bin N sim calves;t À H calves;t À C calves;t ; S sim 1;t ð5Þ where S sim i,t is the estimated survival probability of females of age i at time t. Similarly, N sim 3,t+1 , …, and N sim 14,t+1 (all other terms in Eq. 3) correspond to the population size in the previous age that have survived from time t to time t + 1. Thus, for age j ϵ [2,13]: N sim jþ1; tþ1 $ Bin N sim j;t À H j;t À C j;t ; S sim j;t ð6Þ Note that summer mortality in calves (as well as for other age classes) is considered to be close to zero 46 . Thus N sim,t+1 can be estimated from observed values of ROS, length of the winter, and population size, estimating S and F from the models presented above.
Using the models described above, we estimated age-specific survival rates (from 1994 to 2013), fecundity rates (from 1995 to 2014), and the population size (from 1995 to 2014) through a step-by-step approach. Each year was estimated based on observations of age-specific population sizes (provided by the IPM; Supplementary Fig. 3) and observed ROS and winter length the previous year. Thus, in Eqs. 5 and 6, N sim calves,t was replaced with N obs calves,t and N sim j,t was replaced with N obs j,t . The estimated population size N sim,t+1 may then be calculated step by step (Eq. 3) and compared with the observed population size.
The annual age-specific survival rates ( Supplementary Fig. 1) and fecundity rates ( Supplementary Fig. 2) estimated within this framework were closely correlated to the ones observed (i.e. obtained from the IPM). Pearson's correlations ranged from r = 0.72-0.89 for survival rates and r = 0.65-0.84 for fecundity rates. Accordingly, our model was able to reconstruct annual fluctuations in total population size well ( Supplementary Fig. 4), with a strong correlation (r = 0.89) between estimated and observed population growth rates ( Supplementary Fig. 5).
Reindeer population projections with different ROS scenarios. Finally, we analysed how increasing frequency of rainy winters affects the population growth rate (Fig. 2). To simulate realizations of ROS ( Supplementary Fig. 6), we used the observed ROS data from 1962 to 2014 and an inverse transformation method. Given a continuous uniform variable U in [0,1] and a cumulative distribution function D, the random variable X = D −1 (U) has distribution D. We can change the distribution (i.e. different scenarios of ROS) by using a non-uniform distribution of U. Thus, because the cumulative ROS distribution is bounded between 0 and 1, we used a beta distribution f ðUÞ / U α 1 À1 ð1 À UÞ α 2 À1 where α 1 and α 2 are shape parameters to simulate five different ROS scenarios. We used the following parameters to obtain different ROS scenarios (Supplementary Fig. 6): very high frequency of extreme ROS winters (α 1 = 4, α 2 = 1), high frequency (α 1 = 2, α 2 = 1), low frequency (α 1 = 1, α 2 = 2), and very low frequency (α 1 = 1, α 2 = 4). In addition, the medium frequency scenario simulating the historical state (1962-2014) was obtained by setting α 1 = 1 and α 2 = 1 leading to a uniform distribution as a special case of the beta distribution. Note that 0 < U < 1 for all scenarios. Accordingly, our simulations ( Fig. 3, Supplementary Fig. 6) are conservative as the simulated realizations will fall within the range of observed (historical) values of ROS.
Using the estimated mean age-specific population sizes in 2014 and each of our 9090 age-structured population models (i.e. based on 9090 posterior samples, see Modelling step 2), we simulated for each of the five ROS scenarios (i.e. distributions) 100 stochastic population size trajectories for 100-year projections. This resulted in 4,545,000 population size time-series (454,500,000 population sizes) in total. In other words, instead of using a step-by-step approach (as described above), we fed our age-structured population model with the observed age-specific numbers in 2014 and the ROS trajectories. Note that, for each simulated trajectory, winter length was drawn from a normal distribution with the mean and the standard deviation observed between 1994 and 2014, i.e. our reindeer study period.
We explored how increasing amount of ROS affects the population growth rate, under different age structures and population densities. To do so, we considered two initial age structures (low prime-aged [3-8-year olds] proportion = 28%, high prime-aged proportion = 51%). Note that the two age structures reflect the observed demography before and after the 1996 population crash (Fig. 1b,  Supplementary Fig. 3), respectively. For each age structure, we considered two initial population densities (low N = 1200, high N = 1700). We simulated an increasing amount of ROS (from 0 to 4.2, on a log-scale) from the deterministic version of our population model, keeping winter length constant at its mean observed value . From these population projections, the expected population growth rate for each ROS value was determined for each combination of density/age structure (see Fig. 5).
Finally, we examined how the time elapsed since the previous extreme ROS winter affects the population growth rate during a second extreme ROS winter. An extreme ROS winter reflects the conditions in 1996 (i.e. 1995/96), with the highest ROS value recorded (log ROS = 4.2). We simulated different ROS sequences where ROS was kept at its mean recorded value (1994-2014) except for a first extreme ROS winter always occurring at t = 0, and a second extreme ROS winter occurring at t = 1, then at t = 2, and so on until t = 7. We also included a ROS sequence with no second extreme ROS winter and a sequence with many consecutive extreme ROS winters. Keeping winter length constant at its mean observed value, we then used the deterministic version of our population model to calculate time-series of population sizes (Figs. 6a and 7) and age structures ( Supplementary Fig. 7) following a rather high initial population density (N = 1700) and a low initial proportion of prime-aged individuals (28%), reflecting the population state prior to the 1996 crash. We fed our population model with the eight different sequences of ROS and estimated the expected population growth rate for the second extreme ROS winter (Fig. 6b). All analyses were performed with the statistical software R 47 .
Reporting Summary. Further information on experimental design is available in the Nature Research Reporting Summary linked to this article.