## Abstract

Despite growing concerns regarding increasing frequency of extreme climate events and declining population sizes, the influence of environmental stochasticity on the relationship between population carrying capacity and time-to-extinction has received little empirical attention. While time-to-extinction increases exponentially with carrying capacity in constant environments, theoretical models suggest increasing environmental stochasticity causes asymptotic scaling, thus making minimum viable carrying capacity vastly uncertain in variable environments. Using empirical estimates of environmental stochasticity in fish metapopulations, we showed that increasing environmental stochasticity resulting from extreme droughts was insufficient to create asymptotic scaling of time-to-extinction with carrying capacity in local populations as predicted by theory. Local time-to-extinction increased with carrying capacity due to declining sensitivity to demographic stochasticity, and the slope of this relationship declined significantly as environmental stochasticity increased. However, recent 1 in 25 yr extreme droughts were insufficient to extirpate populations with large carrying capacity. Consequently, large populations may be more resilient to environmental stochasticity than previously thought. The lack of carrying capacity-related asymptotes in persistence under extreme climate variability reveals how small populations affected by habitat loss or overharvesting, may be disproportionately threatened by increases in extreme climate events with global warming.

## 1. Introduction

Despite growing concerns about the increasing frequency of extreme climate events caused by climate change, our empirical understanding of this threat to biodiversity is limited [1–3]. Moreover, it is unknown how extreme climate events will interact with other impacts, such as habitat loss or overharvesting, which reduce population carrying capacity and are likely to intensify as human populations grow [4,5]. As a general rule, populations with small carrying capacities are inherently prone to stochastic extinction due to demographic stochasticity (random independent demographic events) [6–8]. The effect of small carrying capacity on extinction risk will likely be amplified by increased environmental stochasticity (temporal variation in birth and death rates due to environmental perturbations) [9,10], for example, due to increased frequency of extreme droughts. Consequently, the influence of environmental stochasticity on the relationship between carrying capacity and time-to-extinction will be an important mechanism driving population responses to combined land-use and climate change.

The theoretical effects of environmental stochasticity on slopes of the relationship between time-to-extinction of local populations and carrying capacity have been examined extensively, but rarely in empirical studies of actual populations [11]. From theoretical studies, we know that time-to-extinction increases exponentially with carrying capacity in constant environments whereas environmental stochasticity causes the slope of this relationship to decline, becoming asymptotic at low carrying capacity [12–14]. However, recent research shows the effects of environmental stochasticity on populations are often overestimated relative to other sources of demographic stochasticity [15]. Consequently, the impacts of extreme climate events on the relationship between carrying capacity and time-to-extinction may be less dramatic than previously thought. Nevertheless, it remains unknown whether time-to-extinction asymptotes at low carrying capacity under high levels of environmental stochasticity in natural populations. Answers to this question have remained elusive because of the near impossibility of quantifying population demographics simultaneously for many populations of varying size, over gradients of environmental stochasticity, and during rare (e.g. less than 1 in 100 yr) extreme events. These challenges must be overcome to confidently anticipate species extinction risk in an increasingly variable world under climate change.

We investigated the influence of environmental stochasticity on the relationship between time-to-extinction and carrying capacity of local populations in metapopulations of an endangered wetland fish, *Neochanna apoda* (brown mudfish), in New Zealand. *N. apoda* inhabit forest pools that dry frequently and for long durations that can cause up to 84% mortality, as was the case in a recent 1 in 25 year drought [16]. Global climate change models predict increases in extreme drought frequency of up to 45% in New Zealand regions inhabited by *N. apoda* [17,18]. Consequently, droughts affecting local *N. apoda* populations are an extreme form of environmental stochasticity (*sensu* [17–19]), making this system an excellent model to investigate the influences of environmental stochasticity on the relationship between time-to-extinction and carrying capacity. We quantified the effects of environmental stochasticity due to extreme droughts using a large mark–recapture study of over 1 100 *N. apoda* from 41 subpopulations varying widely in size (2–305 fish). Fish survival, recruitment, and dispersal were measured across gradients of drought frequency (0–19 yr^{−1}) and magnitude (0–44 days) over 3 years, including the 1 in 25 year New Zealand-wide extreme drought [16,20,21]. We applied these data to a spatially explicit metapopulation matrix model and classical theoretical models to test the hypothesis that elevated environmental stochasticity due to extreme droughts would create population carrying capacity-related asymptotes in time-to-extinction.

## 2. Material and methods

### (a) Study area and system

Local population dynamics were quantified for *N. apoda* metapopulations in a 9 000 ha temperate peat-swamp-rainforest within South Westland, Tai Poutini National Park, New Zealand (figure 1). Freshwater pools formed by root excavations of large fallen trees are randomly distributed every 5–10 m throughout the forest, and most contain *N. apoda*. A highly variable climate, transitioning from intense rainfall (occasionally exceeding 150–300 mm d^{−1}) followed by days of sunshine cause pools to become dry up to 20 times yr^{−1} for up to 44 consecutive days, with drying frequency and magnitude being negatively correlated with average pool depth (figure 2*a*; [21]). Dispersal of *N. apoda* between pools likely occurs during extreme rainfall events which flood the forest matrix and connect the pools. Living *N. apoda* have also been observed buried within the peaty/mossy matrix between pools, possibly stranded en route to open water during a high-water event or dispersing via subsurface pathways [22]. Thus, the fish assemblage is spatially patchy and likely functions as a metapopulation. The rapid rate of drying, short *N. apoda* lifespan, opportunistic reproductive habits of *N*. *apoda*, and often tiny pool sizes, with high population densities (over 30 fish m^{−3}) make them ideal to quantify the effects of extreme drying simultaneously on many populations of varying size over short time frames [21].

A large-scale mark–recapture study of *N. apoda* has been ongoing in the study area since 31 November 2011 [16,20,21]. Briefly, this work measured *N. apoda* distribution, vital rates, dispersal frequency, and abundance in 41 pools, within four 100 × 20 m forest transects distributed approximately every 2.5 km throughout the forest. Pools within transects were (mean ± s.e.) 56 ± 6 m apart and had surface areas ranging between 0.13 and 28 m^{2} (figure 1). At commencement of the current study, fish in all pools had been marked and recaptured every six months using either visual implant elastomer or passive integrated transponder tags until May 2013, a period that included a range of drought intensities including an extreme 1 in 25 yr^{−1} nationwide drought [21]. Six-monthly survival of *N*. *apoda* during this time, as estimated by Cormack–Jolly–Seber survival models, was positively related to average pool depth, particularly during the extreme drought, due to the declining frequency and magnitude of pool drying as pool depth increased (figure 2*a*,*b*; [21]). Thus, all populations were simultaneously impacted by droughts, but the magnitude of this impact varied spatially with pool depth.

This present work expanded the mark–recapture study until 31 October 2014, to additionally quantify how *N*. *apoda* recruitment and dispersal varied during different drought intensities, and to parametrize the projection matrices in the metapopulation model. All protocols were identical to those described in White *et al*. [16,21], with mark–recapture samples being taken every six months from all 41 pools (i.e. two trapping occasions separated by 14 days every six months). These data enabled us to make seven population size estimates for each pool using closed Lincoln–Petersen models (i.e. one estimate every six months) [23]. In general, population abundance increased with pool surface area according to a logarithmic linear regression relationship which explained over 35% of the spatial variability in population abundance (electronic supplementary material, figure S2). The positive effect of pool surface area on abundance likely results from a greater capture rate of terrestrial invertebrates and leaf subsidies per unit volume in pools with large surface areas, as well as increased space.

In addition to local population size, the mark–recapture data were also used to estimate total metapopulation size and recruitment every six months using open Jolly–Seber (POPAN) models [24]. All Jolly–Seber models were fitted using maximum-likelihood-estimation incorporating spatial and temporal covariates of fish capture probability as described by White *et al*. [20], which minimized uncertainty in vital rates due to measurement error. Meanwhile, the minimum 14 days between mark and recapture was determined empirically in [20] to prevent *N. apoda* trap-avoidance behaviour. These procedures enabled us to achieve the high capture probabilities (average 0.38) [20], advised for precise, accurate estimates of vital rates [25]. During each six-month period, the maximum rainfall deficit (difference between evapotranspiration and rainfall in millimetres) was recorded from the nearest National Climate Database rain gauge (#4054: http://cliflo.niwa.co.nz/), as a measure of total system drought strength, with high numbers reflecting dryer conditions than usual. Deficit values occurring at frequencies less than once every 20 years in a particular region are considered extreme in global climate change projections for New Zealand [17,18]. These data enabled us to determine how total metapopulation recruitment varied with system dryness over time.

### (b) Metapopulation matrix model structure

Our spatially explicit matrix model of *N. apoda* metapopulations was constructed using RAMAS Metapop software [26]. RAMAS Metapop has been widely used to model metapopulation dynamics for applied and research purposes and models the effects of both environmental and demographic stochasticity on population dynamics [27–29]. Environmental stochasticity was defined as temporal fluctuations in birth and death rates caused by environmental perturbations, which can include droughts provided populations are not extinguished by a single drought [13,19]. *N. apoda* survival probability during a 1 in 25 year drought was measured as at least 16% [21] which constitutes an extreme form of environmental stochasticity (*sensu* [13,19]). Demographic stochasticity was defined (*sensu* [13,19]) as chance events affecting individual mortality and reproduction.

Environmental stochasticity due to droughts was modelled using the ‘matrix selection’ procedure [30], whereby projection matrices representing either low, moderate, or extreme drought vital rates (survival and recruitment) were randomly selected at each time step to be used in population projection (figure 3*a*). The projection matrices reflected mean population recruitment and survival rates over six-month time periods when maximum drought intensity was either: low (6.4–36 mm rainfall deficit), moderate (36–94 mm rainfall deficit), or extreme (more than 94 mm rainfall deficit), which occur at probabilities of 0.70, 0.28, or 0.02, respectively, according to historic climate data (National Climate Database, NIWA, electronic supplementary material, figure S1*a*). Vital rates in these matrices were the mean probability of survival and *per capita* number of new recruits for each age class over six months for a pool of average depth, and assumed a mixed sex ratio of 1 : 1 as according to past research on *N. apoda* populations [31]. Uncertainty in vital rates was assumed to be caused by non-drought sources of environmental stochasticity, and this uncertainty was characterized using a standard deviation matrix, which allowed vital rates to vary within each drought category. Because these standard deviations could conceivably be confounded with mark–recapture measurement error, despite our efforts to minimize it, this assumption was later scrutinized in sensitivity analyses. Spatial variation in environmental stochasticity was induced by constraining vital rates as a positive function of pool depth according to drought-specific relationships fitted to the mark–recapture data using maximum-likelihood-estimation in White *et al*. [21] (figures 2*b* and 3*b*). Consequently, all populations were exposed to the same type of drought intensity on each time step (i.e. low, moderate, or extreme) but survival and recruitment was lowest in populations inhabiting shallow pools for any given drought intensity (figure 2*b*). The historic drought time sequence since 1965 showed little evidence of serial autocorrelation (electronic supplementary material, figure S1*a*,*b*), suggesting that modelling drought recurrence as serially independent adequately captured the system's dynamics.

Recent research showed the effects of environmental stochasticity can be influenced by sources of demographic heterogeneity, i.e. variation in fitness among individuals within a population due to body size, age, genetic endowment, etc. [15]. To account for this, demographic heterogeneity was introduced by constraining survival within each projection matrix as a function of age according to *N. apoda* length–survival relationships empirically derived during each drought in White *et al*. [21]. Consequently, demographic heterogeneity varied according to drought intensity, with survival increasing rapidly with age during a low drought, while being independent of age during moderate and extreme droughts [21]. Age-independent *N. apoda* survival during moderate and extreme droughts implies equivalent drought tolerance across cohorts [21]. Lengths were converted to ages for the projection matrices using the otolith-based length–age relationship reported for *N. apoda* by Eldon [31]. Within each drought category, recruitment remained uniform across ages due to insufficient information on how age-recruitment dependency changed with drought conditions. Consequently, we later investigated the sensitivity of model results to uniform versus positive age-dependent recruitment.

Demographic stochasticity was modelled by selecting the number of survivors and recruits of the *i*th age class from a binomial or Poisson distribution, respectively, with the mean survival or recruitment rates, and *N _{i}*(

*T*) as parameters (figure 3

*c*). Consequently, even if average survival probability during a drought,

*p*, was 0.5, there was a probability, , that

*all N*(

_{i}*T*) individuals would survive to

*i*+ 1, and a probability that none survive [32]. This probability is exponentially higher when

*N*(

_{i}*T*) is small thereby increasing the probability of population extinction, in accordance with the scaling of demographic stochasticity with population carrying capacity [9,32].

Inter-population dispersal probability (figure 3*d*) decayed with distance according to Kitching's [33] dispersal kernel: *m* = *a* × e(–*D ^{c}*/

*b*), where

*m*was fish movement probability,

*D*was distance (metres), and

*a*,

*b*, and

*c*were constants controlling the intercept and form of the relationship, as estimated from the mark–recapture data. The maximum dispersal distance was 120 m thus allowing fish to disperse between pools within transects, and not between transects. Although

*N. apoda*may disperse distances greater than 120 m over longer time scales, using pools as stepping stones, this would only allow them to move up to 1 920 m over their 8-year lifespan, which is less than the minimum distance between transects.

Local population carrying capacity was represented in the model as a constant upper limit of population size, above which population growth rate became negative due to density dependence. Density dependence (figure 3*e*) was modelled using the ceiling function in RAMAS. This allowed exponential growth, depending on environmental and demographic stochasticity until carrying capacity was reached, and negative growth thereafter [26]. In this manner, our representation of carrying capacity and density dependence was in line with that represented in classical theoretical population models that quantified the relationship between local time-to-extinction and carrying capacity [13]. Fish in populations closer to carrying capacity were also increasingly likely to emigrate due to linear density-dependent dispersal (figure 3*e*). Density-dependent emigration slopes are automatically determined by RAMAS for each population separately to ensure positive emigration probabilities for all abundances more than 1 [26].

Population carrying capacity was estimated as the maximum recorded abundance out of seven Lincoln–Petersen mark–recapture estimates taken per population, which varied spatially with pool surface area (electronic supplementary material, figure S2). To check whether using maximum abundance gave undue influence to outliers in the abundance data, we reran the model after re-parametrizing carrying capacities using the 95th quantile of the relationship between pool surface area and population abundance (electronic supplementary material, figure S2). Using a high abundance value, less than the observed maximum, recognizes that populations occasionally overshoot carrying capacity. We also investigated model sensitivity to temporal variation in carrying capacity by defining standard deviations for all population carrying capacities (i.e. carrying capacity = maximum observed abundance ± s.d.). These standard deviations were estimated from the error in the relationship between population abundance and pool surface area (electronic supplementary material, figure S2).

### (c) Matrix model training, validation, and sensitivity analyses

The metapopulation matrix model was validated by comparing the relationship between time-to-extinction and carrying capacity predicted by the model, to that observed in the mark–recapture data. To ensure an independent validation, we temporally split the mark–recapture data into training data, used to parametrize the model, and validation data used to evaluate model predictions. Matrix vital rates, population carrying capacities, and dispersal probabilities for training were estimated using mark–recapture data collected between 30 November 2011 and 2 June 2013. The model trained by these data was run for 1 000 repetitions of 500 time steps (250 yr at 0.5 yr time step). Owing to the upper time step limit in each model run, time-to-extinction (where abundance < 2) of local populations was calculated as the number of time steps-to-extinction divided by the maximum run length (500 time steps of the model), and is thus expressed as probability of local population persistence. Probability of population persistence in the validation dataset was calculated as the number of sampling occasions fish were encountered in each pool out of the seven total samples taken between 3 June 2013 and 31 October 2014. Using these data, we observed whether there were significant differences in the slope of the relationship of local population persistence and carrying capacity by examining the interaction between carrying capacity and data type (i.e. observed versus predicted datasets) in a quasi-binomial generalized linear model. One population had to be removed from this analysis due to the lack of any fish captured in the last 1.5 years (population 7.14, transect 7). One fish (tag ID: 8 636) was captured in this pool in the first 1.5 years, but disappeared thereafter and was not recaptured in any other pool, indicating it had most likely deceased or emigrated from the study area.

Additional validation procedures involved comparing differences in intercepts and slopes of the relationship between observed and predicted mean population abundance using the independent training and validation datasets. Observed abundance was calculated as the mean of the Lincoln–Petersen mark–recapture abundances measured during the final 1.5 years of mark–recapture, while predicted abundance was calculated as the mean final population abundance estimated using the metapopulation matrix model. Per cent bias (*sensu* [34]), the average tendency for predicted values to be larger or smaller than their observed counterparts as a percentage thereof, was also calculated.

Pending model validation, a final metapopulation matrix model was re-parametrized using the entire mark–recapture dataset. Analyses were conducted on this model to determine how sensitive the slope of the relationship between local population persistence and carrying capacity was to parameter and structural uncertainty. Thus, we decreased all projection matrix vital rates, individually, to their lower 95% confidence limits as quantified using the mark–recapture models. Spatial and temporal uncertainty in carrying capacities were also investigated as described above. The decision to use uniform recruitment-at-age was scrutinized by comparing results to a model incorporating a positive relationship between age and recruitment estimated from Eldon [31] (electronic supplementary material, figure S3). Finally, we examined sensitivity to increases in moderate and extreme drought probability. Sensitivity was quantified using the percentage change in the slope of the relationship between local population persistence and carrying capacity before and after parameter/setting adjustment using quasi-binomial generalized linear models.

### (d) Effects of environmental stochasticity and carrying capacity on time-to-extinction

We tested the hypothesis that elevated environmental stochasticity due to extreme droughts would create carrying capacity-related asymptotes in time-to-extinction using both the final metapopulation matrix model and by applying our estimates of environmental stochasticity to the classical population model described by Lande [13]. These two avenues of modelling were chosen so that changes in the scaling of time-to-extinction with carrying capacity could be observed at both small and large temporal scales (e.g. 250 yr for the matrix model and infinite for the Landes model). Matrix model simulations were run and time-to-extinction was calculated as the probability of local population persistence as described above. Dispersal probability was set to 0 for these simulations so that population persistence was estimated for isolated populations in keeping with classical studies examining the relationship between carrying capacity and time-to-extinction under environmental stochasticity [12–14]. A quasi-binomial generalized linear model was used to analyse the effect of carrying capacity (log_{e0} carrying capacity), environmental stochasticity and their interaction, on local population persistence.

For the quasi-binomial model, environmental stochasticity was measured for each local population as the variance in *per capita* population growth rate (*r*: Δ*N _{t}*

_{+1}/

*N*), resulting from the random draws of different drought strength projection matrices. In a typical matrix model run involving consecutive time steps, population growth rate quantified from projection matrices is confounded with population size for a given time step, which varies over time independently of environmental stochasticity [35]. To avoid such confounding effects,

_{t}*per capita*population growth rates were quantified from a sample of 10 000 independent projection matrices per population, with initial population sizes set at 10 000 and carrying capacity set to 20 000 (electronic supplementary material, figure S4). In this manner, our measure of environmental stochasticity was equivalent to the variance in population growth rate caused by the effects of environmental variability on individual fitness

*sensu*Lande [13] and Lande

*et al.*[19].

Finally, we applied our estimates of environmental stochasticity for each local population, measured above, to a theoretical population model that quantified the form of the relationship between local time-to-extinction and carrying capacity under environmental stochasticity. Several authors have shown that time-to-extinction scales as a power of carrying capacity (*K*) proportional to *K ^{c}*, where

*c*is 2

*r*/

*V*– 1, and where

_{r}*r*and

*V*are the mean and environmental variance in population growth rate, respectively [12–14]. The

_{r}*c*parameter quantifies the form of the relationship between time-to-extinction,

*T*and carrying capacity,

*K*according to the equation (2.1) presented by Lande [13]: 2.1

Thus, if environmental variance in growth rate is sufficiently higher than the mean growth rate (i.e. *c* < 1), time-to-extinction asymptotes at decreasingly lower *K*. Otherwise, time-to-extinction increases exponentially with *K* (i.e. when *c* > 1) [13]. We used our estimates of *r* and *V _{r}* from the matrix model to quantify

*c*for all populations, and observed the relationship between

*K*and time-to-extinction for the populations with the highest and lowest

*c*using equation (2.1).

## 3. Results

### (a) Metapopulation model training, validation, and sensitivity

Using data from the entire duration of the mark–recapture study, the six-monthly *per capita* recruitment estimates from the Jolly–Seber mark–recapture model averaged 0.37 (±0.06 s.d.), 0.13 (±0.02 s.d.), and 0.13 (±0.01 s.d.), during the low, moderate, and extreme drought strengths, respectively. Thus, recruitment during moderate and extreme droughts was consistently 35% of that during low drought strengths; even moderate drought strengths curtailed recruitment, and this was not worsened by the 1 in 25 year extreme drought. To reduce model complexity, we modelled recruitment in both moderate and extreme droughts as a consistent proportion of recruitment during low droughts. Meanwhile, total metapopulation size and local population carrying capacity averaged 596 (±97 s.d.) and 25 (range: 1–305), respectively, as estimated from the Jolly–Seber and Lincoln–Petersen mark–recapture models, respectively. A total of 97 fish movement events, ranging from 3 to 112 m were recorded during the entire duration of the mark–recapture study, and all observations were grouped into 3 m distance bins. These bins were multiplied by 1.62 to adjust for unobserved movements due to the probability of fish being undetected, which averaged 0.62 according to White *et al*. [20], thus yielding a total number of non-dispersing fish of 439 (i.e. 596 – 157). Dispersal probability using this entire dataset was a negative function of distance, best fit by the equation: *y* = 9.68 × e(–*D*^{0.18}/0.10) (electronic supplementary material, figure S5). Recruitment, dispersal, and population abundance statistics changed only slightly when mark–recapture analyses were restricted to the first 1.5 years for model training and validation purposes (electronic supplementary material, table S1).

When parametrized using vital rates and dispersal data measured during the first 1.5 years of the study, the predicted relationship between carrying capacity and persistence of local populations closely resembled that observed during the final 1.5 years of study (electronic supplementary material, figure S6). A quasi-binomial regression on these data indicated no significant differences between predicted and observed datasets, in terms of either slopes (*F*_{1,78} = 0.16, *p* = 0.69) or intercepts (*F*_{1,78} = 12.5, *p* = 0.60). Mean observed abundance was also closely related to predicted abundance according to an intercept and slope of 2.2 and 0.94, neither of which were significantly different from 0 (*t*_{1,39} = –0.92, *p* = 0.36) and 1 (*t*_{1,39} = 1.77, *p* = 0.09), respectively (electronic supplementary material, figure S7). Although the intercept was nearly significantly different, predicted abundance values were only negatively biased by 1.47 fish; about 9.99% on average, compared with observed values. Consequently, the metapopulation matrix model closely approximated the relationship between local population carrying capacity and persistence we independently observed in real populations, and accurately predicted independently observed abundance. Therefore, the final model was parametrized using means and standard deviations reported for the full mark–recapture dataset (electronic supplementary material, S1).

Sensitivity analysis indicated that the slope of the relationship between local population carrying capacity and persistence was sensitive to reductions in low drought period recruitment and survival, moderate drought period survival and removal of the standard deviation matrix which changed by 23, 21, 28, and 16%, respectively (electronic supplementary material, table S2). However, the changes were only significant for moderate drought survival (electronic supplementary material, figure S8). Slopes were sensitive to structural alterations including altered population carrying capacities, temporal variation in carrying capacity, age-dependent recruitment, and increased moderate drought probability (electronic supplementary material, table S2). However, the changes were not significant (electronic supplementary material, figure S9), thus our conclusions regarding the scaling of local population persistence with carrying capacity were robust to structural parameter uncertainty.

### (b) Effects of environmental stochasticity and carrying capacity on population time-to-extinction

The relationship between local population persistence and carrying capacity was significantly reduced by elevated environmental stochasticity as shown by a significant interaction between population carrying capacity and environmental stochasticity (carrying capacity effect: *F*_{1,37} = 2 285.2, *p* < 0.001, *R*^{2} = 0.92, environmental stochasticity effect: *F*_{1,37} = 53.5, *p* < 0.001, *R*^{2} = 0.02, interaction effect: *F*_{1,37} = 22.8, *p* < 0.001, *R*^{2} = 0.01). This interaction was significant regardless of parameter uncertainty (electronic supplementary material, table S2). Thus population persistence increased with carrying capacity, and elevated environmental stochasticity caused the slope of this relationship to decrease (figure 4*a*). However, despite this decrease in slope, even the highest recorded level of environmental stochasticity had surprisingly minor impact on persistence of the largest population (*K* = 305), which, on average went extinct for only 1.4 time steps (out of 500) per model run (figure 4*a*). Thus, environmental stochasticity had a surprisingly marginal impact (figure 4*b*), with vulnerability to droughts being almost entirely dependent on carrying capacity (figure 4*a*).

Mean and environmental variance of *per capita* population growth rate (*r*) in the matrix model ranged between 0.06–0.13 and 0.009–0.035 Δ*N _{t}*

_{+1}/

*N*, resulting in

_{t}*c*ranging between 2.53 and 26.80. Although the slope of the relationship between time-to-extinction and carrying capacity declined predictably as environmental stochasticity increased, this was not sufficient to cause the relationship to asymptote at low carrying capacity (figure 5). This conclusion was robust to model uncertainty as indicated by minimum

*c*> 1 for all model parameters and settings for which sensitivity was investigated (electronic supplementary material, table S2). Thus, time-to-extinction increased indefinitely with carrying capacity for all populations, even those experiencing up to 19 complete drying events per year. Consequently, the highest levels of environmental stochasticity measured in this study were insufficient to curtail time-to-extinction of populations with large carrying capacities.

## 4. Discussion

Using empirical estimates of environmental stochasticity, this study confirmed the long-standing theory that slopes of the relationship between time-to-extinction and carrying capacity decrease as environmental stochasticity increases [11]. However, we further showed that the decrease in slope caused by environmental stochasticity was far lower than that required to substantially impact large populations. Under sufficiently high levels of environment stochasticity, theory predicts that time-to-extinction asymptotes at relatively low carrying capacities, such that even large populations may go extinct in environmentally variable habitats [12–14]. In our case, time-to-extinction increased with population carrying capacity, and the slope of this relationship decreased for populations experiencing the highest environmental stochasticity. However, the decrease in slope was not sufficient to cause substantial declines in time-to-extinction for the largest populations, despite some pools drying multiple extended periods a year. This means that populations with large carrying capacities may be more resilient to extreme climate event frequency than previously thought.

The lack of carrying capacity-related asymptotes in time-to-extinction was surprising given the extreme environmental variability affecting the populations we studied. The population with the largest carrying capacity (305) was also the worst impacted by droughts, with population abundance dropping over 80% during the most extreme drought. Declines in population abundance >80% are usually grounds to increase extinction risk rankings for populations that are normally too large to be considered ‘critically endangered’ under International Union for the Conservation of Nature (IUCN) extinction risk criteria [36]. Such criteria are designed to incorporate the uncertainty in extinction risk resulting from asymptotic scaling of time-to-extinction with carrying capacity when large populations are exposed to extreme environmental variability [13]. However, the mean and variance in population growth rate in this affected population remained relatively stable at 0.06 and 0.04 Δ*N _{t+}*

_{1}/

*N*, respectively, values which were considerably different to those theoretically required to cause such asymptotic scaling relationships [12–14]. The high population stability in the face of extreme environmental variability we observed is very likely a product of life-history and physiological adaptations, such as cutaneous respiration and burrowing [37], which reduces the effects of environmental extremes on birth and death rates in similar emersion-tolerant species [38]. Paradoxically, this implies that extinction risk criteria may be relaxed for populations in extremely variable environments, provided they are well adapted to disturbances.

_{t}A particular strength of our study is the size and detail of the dataset used to drive empirical modelling. Numerous pools which varied in size, but were easily sampled combined with the serendipitous large drought that occurred during our study, made it possible to measure population dynamics for many populations varying widely in size and environmental stochasticity, and during a rare extreme event [21]. Moreover, the high capture probability of *N. apoda* enabled us to isolate the effects of environmental variation on population dynamics, leading to higher precision and accuracy of vital rate estimates, and the effects of drought thereon [20]. Many classic population studies have likewise benefitted from similar ‘island’ topography, enabling precise measures of population dynamics, and over longer time scales than studied here [39,40]. However, unlike most detailed long-term population studies [41,42], our data were not limited to one or a few populations. The species-specific data driving our models may limit the generalizability of the inferences they provide, but the large number of populations varying in carrying capacity and exposure to environmental stochasticity studied here offers rare insights into key theoretical and practical issues, which would not be available from an examination of a single population.

The role of environmental stochasticity in driving population persistence has been heavily debated [6,12,13], particularly recently given the rising concern over the increase in extreme climate event frequency projected to occur with climate change [3,10,43]. However, recent work shows the effects of environmental stochasticity are often overestimated relative to other stochastic drivers of population dynamics [15]. Our research paints a slightly different picture, showing that while the importance of environmental stochasticity may be less than expected, the degree to which this is the case is dependent on carrying capacity. The lack of asymptotic scaling of time-to-extinction with carrying capacity found here under extreme climate variability reveals how large populations may be more resilient to extreme climate variability than previously thought. However, small populations remained highly vulnerable in variable environments. Consequently, the threat of increasing climate extremes likely remains significant, but may be more visible or acute where populations are small or shrinking, such as those constrained by habitat loss and exploitation [4,44].

## Ethics

All procedures were approved by the University of Canterbury Animal Ethics Committee (permit no. 2010/23R).

## Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material or can be found in [45].

## Authors' contributions

R.S.A.W., A.R.M., and P.A.M. conceived and obtained funding for the mark–recapture study. R.S.A.W. and B.A.W. conceived the metapopulation matrix model. R.S.A.W. collected/analysed data, constructed and ran matrix models under guidance from A.R.M., B.A.W., P.A.M., and D.J.B. R.S.A.W. wrote the manuscript with contributions from all authors. A.R.M. was project leader, organizing resources and permitting, and guiding overall research direction.

## Competing interests

We declare we have no competing interests.

## Funding

R.S.A.W. was supported by a UC Scholarship during initial sampling and by the UC Roper Scholarship in Science during the later stages. B.A.W. was supported by ARC Future Fellowship FT100100819.

## Acknowledgements

We thank Alan Lilley, Amy Whitehead, and members of the Freshwater (FERG) and Quantitative Ecology (QAECO) Research Groups at the Universities of Canterbury and Melbourne for advice and assistance during this project. We also thank Mike Winterbourn, George Perry, Nick Bond, and the two anonymous reviewers for their comments on earlier versions of this paper. This project was initially funded by the Brian Mason Scientific and Technical Trust with subsequent contributions by the New Zealand Department of Conservation (DOC), the NIWA Sustainable Water Allocation Programme, Northwest Marine Technology Inc. for elastomer tagging equipment, and the University of Canterbury (UC). We are particularly grateful to the University of Canterbury for use of the Harihari field station facilities and DOC for permission to work in Westland Tai Poutini National Park.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3789994.

- Received April 18, 2017.
- Accepted May 15, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.