How density dependence, genetic erosion and the extinction vortex impact evolutionary rescue

(a) Model description

We explored evolutionary rescue and extinction using stochastic simulations. This approach allowed us to incorporate several stochastic demographic and genetic processes that are important for rescue or extinction [2,26]. We derive a discrete-time individual-based model with diploid genetic inheritance that has general features of the semelparous life-history strategy common to many sexually reproducing organisms. We use non-overlapping generations for comparison with other models of adaptation and evolutionary rescue for semelparous organisms (e.g. [1,9,23]). Our genetic submodel is a finite quantitative locus model with random and independent segregation of alleles, similar to that of Boulding & Hay [23]. Our demographic submodel is based on an expanded stochastic Ricker model [26,27], which models density dependence.

Each individual has an intrinsic fitness determined by genetic and environmental components. The genetic component is determined by each of m biallelic diploid loci. Each allele-copy contributes a positive (1/(2√m)) or negative (−1/(2√m)) amount to the genotype of individual i, g_i, which is the sum of allele-copies across all loci. From these allele frequencies, we estimate the additive genetic variance in the population, σ_a², following [28], as shown in electronic supplementary material, SA). We do not model mutations, since populations are small and timespans are short in our simulation experiment (see below), such that at realistic mutation rates mutation would be rare in our simulations. An individual's phenotype, z_i, is the genotype, g_i, plus a deviation due to environmental variation drawn from a normal distribution with mean 0 and variance ${\sigma }_{\mathrm{e}}^2$. Intrinsic fitness, W_i, the expected fitness of an individual in the absence of density effects, is determined by Gaussian selection [10] such that $W_i=W_{\max }\hspace{0.17em}\exp (-{(z_i-\theta )}^2/2w^2)$, where θ is the optimal phenotype, w² is the variance (w is the width) of the fitness landscape and is inversely proportional to selection strength, and W_max is the maximum intrinsic fitness (i.e. intrinsic fitness of an individual with the optimal phenotype). A change in θ models a change in the environment. The expected number of offspring of an individual, R_i, is determined according to a Ricker function that reduces fitness due to density dependence such that $R_i=W_i{\mathrm{e}}^{-\alpha N}$, where α is the strength of density dependence and N is the population size [27]. Every generation, each female mates with one male chosen uniformly at random from the population. Each female produces a Poisson-distributed number of offspring with mean 2R_i, with parental alleles segregated independently and at random to offspring. Sex is assigned to each individual upon birth with equal probability. Generations are not overlapping, i.e. only one generation is present in each time step. As such, if all females have zero offspring, or if all individuals in the population are of the same sex, then extinction is guaranteed in the next time step.

(b) Experimental design, simulations and analysis

While extinction is straightforward to identify, the classical criterion of evolutionary rescue (population growth rate above replacement) is more difficult to operationalize for individual populations because rescued populations may transiently decline in size or growth rate owing to demographic or genetic stochasticity and revert to an ‘unrescued' state. We operationalize rescue using two separate criteria: (i) expected number of offspring per individual (${\overline{R}}_t$, as defined below) greater than 1 for three consecutive generations (hereafter ‘fitness-based rescue') and (ii) populations exceeding their initial size for three consecutive generations (hereafter ‘size-based rescue').

Parameter values are given in electronic supplementary material, table SB1, justified in electronic supplementary material, SC, and are realistic for natural populations [29,30]. We simulated populations with m = 25 loci and an environmental change that moved the phenotypic optimum from its prior state of 0 to θ = 2.8. Because θ > 0, the environment in our simulations favours a shift toward positive alleles. Populations were founded with positive and negative alleles assigned to each individual's loci with equal probability such that across simulation trials the expected population mean genotype was 0, although owing to sampling there was variation in initial population mean genotype across trials. Thus, to reach the new phenotypic optimum θ = 2.8, a typical population needed to change from 25 positive allele-copies (out of 50) per individual to 39 allele-copies per individual (see electronic supplementary material, SD). All simulations were run with maximum intrinsic fitness W_max = 2, non-genetic phenotypic variance ${\sigma }_{\mathrm{e}}^2=0.5$ and w² = 3.5 as the variance of the fitness landscape. With these parameter values, mean intrinsic fitness in the typical population immediately following the environmental change was approximately 0.65, such that populations were sufficiently challenged to the point that both rescue and extinction were possible on rapid timescales.

To evaluate the effects of negative density dependence on evolutionary rescue and its potential interactions with initial population size and genetic diversity, we conducted a multi-factorial simulation experiment (2 × 2 × 2, i.e. eight treatments). Populations were simulated with the following parameter values: negative density dependence absent (α = 0) or present (α = 0.0035, carrying capacity approx. 165, see electronic supplementary material, SC), hereafter referred to as ‘density-independent populations' or ‘density-dependent populations'; initial size N₀ = 100 or N₀ = 20, hereafter referred to as ‘large' or ‘small' populations in reference to their initial state rather than their eventual size; initial genetic diversity high (${\sigma }_{\mathrm{a}}^2\approx 0.5$, heritability approx. 0.5) or low (${\sigma }_{\mathrm{a}}^2\approx 0.25$, heritability approx. 0.33 [28]), hereafter ‘high diversity' or ‘low diversity' populations, again in reference to their initial state rather than their eventual level of diversity. These initial heritabilities match the range of heritabilities in life-history traits in wild populations [31]. Note that large populations do not necessarily harbour more genetic diversity than small populations independently of the genetic diversity treatment. To simulate populations with initially low diversity, at the beginning of each simulation six loci were set to fixation for the positive (adaptive) allele and six were set to fixation for the negative (maladaptive) allele, leaving only 13 loci (approx. half of the genome) available for selection to act upon while leaving the expected mean population genotype unchanged compared with the high-diversity treatments.

To assess adaptation and extinction on short timescales, we ran 4000 simulations lasting 15 generations for each treatment. All simulations were performed in R, version 4.1.2 [32]. In each simulation, for each generation we recorded population size (N_t), mean genotype (${\overline{g}}_t$) and phenotype (${\overline{z}}_t$), mean intrinsic fitness (${\overline{W}}_t$), additive genetic variance (${({\sigma }_{\mathrm{a}}^2)}_t$) and the proportion of loci at fixation within the population for either allele. Throughout our analysis, we will use an overbar to represent a mean across individuals within one population in one time step. We define adaptation as an increase in ${\overline{W}}_t$. Because fitness is optimized at a positive genotypic value (i.e. θ > 0) and fitness is a monotonic function of maladaptive genetic load, adaptation can be equivalently defined as an increase in ${\overline{g}}_t$ or a reduction in genetic load, $(\theta -{\overline{g}}_t)$.

As a first evaluation of how density dependence, population size, and genetic diversity influence rescue and extinction over 15 generations, we visualized mean population size over time for each treatment. This visualization shows whether the 'U'-shaped population trajectory characteristic of rescue [1] occurs, and if so, how it differs by treatment. All mean population sizes were estimated including extinct populations as size 0, thus incorporating the possibility of extinction into expected population size in each time step. To assess the effects of genetic erosion on population size and extinction, we compared mean population size from each simulated treatment with numerically estimated population size predicted by Gomulkiewicz & Holt's model [1] (in which phenotypic variance and rate of adaptation are held constant), modified to include a Ricker-density dependence term:

where k₀ is the rate of phenotypic change (as described in [1]) in the first time step, averaged across all simulation trials within a genetic diversity treatment. We also estimated and visualized the following response variables, averaged over the ensemble of all populations within a treatment: mean population genotype, $⟨{\overline{g}}_t⟩$, mean intrinsic fitness, $⟨{\overline{W}}_t⟩$, mean additive genetic variance, $⟨{({\sigma }_{\mathrm{a}}^2)}_t⟩$ and number of loci at fixation for positive and negative alleles across all populations in each generation. To differentiate ensemble means from population means for a single trial, we use angle brackets ($⟨\hspace{0.17em}⟩$) to denote means of population-level means across an ensemble of simulations. Because these variables are undefined for extinct populations (which have no mean genotype, fitness or additive genetic variance), estimating unconditional expectations for each generation was impossible, and thus we estimated these means conditioned on survival or extinction at the end of the simulation and interpret plots with caution (see electronic supplementary material, SE).

To assess the probabilities of extinction and rescue on longer timescales, we ran 1000 simulations per treatment lasting 50 generations. Owing to the computational intensity of simulating large populations, we stopped a trial upon population size exceeding 10 000. Because of this censoring, we did not analyse trajectories of demographic or genetic state variables in these longer simulations, instead focusing only on the probabilities of extinction and rescue. We evaluated extinction in each generation in two ways: (i) instantaneous probability of extinction (the probability of extinction in generation t conditioned on surviving to t) and (ii) cumulative probability of extinction over time (the probability of extinction in or prior to t). We also used two separate criteria of rescue (defined above): (i) fitness-based rescue and (ii) size-based rescue. We also estimated time to rescue across treatments, defining time to rescue as the first of the three generations where the criterion was met (e.g. a population exceeding original size in generations 10–12 would be considered rescued in generation 10).

We estimated treatment effect sizes on the probabilities of extinction and rescue in the 50-generation trials using Bayesian generalized linear models (GLMs) in the R package rstanarm, 2.21.1 [33]. We fitted separate models for estimating probabilities of extinction and both types of rescue, treating each simulated population as a Bernoulli response and each of our three treatments (density dependence, initial size, genetic diversity) as categorical variables in each model. Additionally, we estimated effects of initial degree of maladaptation, an important determinant of extinction risk [1,24], by including the mean population genetic load in the initial generation ($\theta -{\overline{g}}_0$, henceforth, ‘maladaptation') as a continuous explanatory variable. This was possible for our stochastic model because sampling variation in the founding of populations produced random among-trial variation in mean initial genotype. We estimated effects of negative density dependence and its interactions with other variables by fitting a model with a four-way interaction and examining effect sizes, estimated as described in electronic supplementary material, SF. The uncertainty around these estimates was sufficiently low and could be arbitrarily improved by running more simulations. An effect size for a categorical variable (for example, initial population size) gives the effect of changing that variable's parameter on the log odds of extinction (or rescue) while holding all other parameters constant, while an effect size for maladaptation gives the change in the log odds of extinction (rescue) associated with switching five alleles from the positive allele to the negative allele in all individuals.

To test for the presence of an extinction vortex, we first reversed and re-centred time for each extinct population to give τ, the time until extinction, such that populations had a common origin (τ = 0) with consistent meaning across populations [34]. We define τ such that τ = 0 is the first time step where the population is extinct (i.e. population size of 0 or 1), τ = 1 is the last time step in which the population size is above 0 or 1, etc. The extinction vortex is characterized by accelerating declines in population size, genetic variation and/or rate of adaptation as τ approaches 0 [25,34]. Specifically, we calculated the rate of population growth ($r_{\tau }=\ln (N_{\tau }/N_{\tau +1})$), the proportional loss of genetic diversity (${\nu }_{\tau }=1-({({\sigma }_{\mathrm{a}}^2)}_{\tau }/{({\sigma }_{\mathrm{a}}^2)}_{\tau +1})$), and the rate of adaptation ($1-k_{\tau }$). $k_{\tau }$ is defined as in [1] such that $k_{\tau }=(\theta -{\overline{g}}_{\tau })/(\theta -{\overline{g}}_{\tau +1})$; as $({\overline{g}}_{\tau }-{\overline{g}}_{\tau +1})$ approaches 0 (i.e. as the rate of genotypic change approaches 0), $k_{\tau }$ approaches 1 and $1-k_{\tau }$ approaches 0. If adaptation proceeds at a constant rate, then $k_{\tau }$ will be constant over time. We used these normalized rates of change instead of non-normalized state variables (e.g. $r_{\tau }$ instead of $N_{\tau }$) to facilitate comparisons across treatments and to remove potential biases due to time until extinction. We estimated the ensemble means $⟨r_{\tau }⟩$, $⟨{\nu }_{\tau }⟩$ and $1-⟨k_{\tau }⟩$ for all extinct populations in each treatment using the 15-generation simulations, where all state variables were tracked for each population. Following [34], we considered decreasing $⟨r_{\tau }⟩$ (accelerating decline in population size), increasing $⟨{\nu }_{\tau }⟩$ (accelerating loss of genetic diversity) or decreasing ($1-⟨k_{\tau }⟩$) (slowing rate of adaptation) in the lead up to extinction (i.e. as τ approaches 0) to be evidence of an extinction vortex. Likewise, we assessed how the shape of the extinction vortex varies between density-dependent and density-independent populations qualitatively by noting differences in the mean trajectories of the respective experimental treatments. We did not estimate rates at τ = 0 because$\hspace{0.17em}r=-\mathrm{\infty }$ for a population at size 0 and because additive genetic variance and mean genotype are undefined for extinct populations.

(a) Population size

For all combinations of size and genetic diversity, in each generation density-dependent populations were smaller on average than corresponding density-independent populations (compare black and purple lines, figure 2). The minimum average size (bottom of the ‘U') reached by density-dependent populations was smaller than that of density-independent populations across treatments, with larger effects in large populations (two-thirds reductions in minimum size) than small populations (one-third reduction). Density-dependent populations were also smaller on average at the end of 15 generations; populations with low diversity had slightly larger reductions in final size (92–95% reduction) than populations with high diversity (85–90% reduction). Regardless of starting population size, density-dependent populations declined for longer and took longer to recover their original size (dotted lines, figure 2), if at all. These results are robust to the parameter values we use; electronic supplementary material, SG shows qualitatively similar patterns with different parameter values and functional forms of density dependence.

• Download figure
• Open in new tab
• Download PowerPoint

We observed an interaction between density dependence and changing genetic diversity on population size over time. Under density dependence, simulated populations (with changing additive genetic variance) were smaller compared with a model with constant additive genetic variance by the end of simulations (compare thick and thin purple lines, figure 2). However, the opposite result was observed for density-independent populations: populations with changing genetic variance were larger than with constant genetic variance (compare thick and thin black lines, figure 2). These effects of changing genetic variance emerged for density-independent populations after populations reached their minimum size but took longer to emerge in density-dependent populations. The combined effect was that the effect of density dependence on mean size after 15 generations was larger with genetic erosion than with constant genetic diversity.

(b) Extinction and rescue

These differences in population size and growth rates over 15 generations translated into more frequent extinctions for populations with density-dependent growth over the course of the 50-generation simulations (figure 3). For both density-dependent and density-independent populations, most extinctions occurred in the first 20 generations (figure 3a,b). Density dependence increased extinction probabilities primarily in large populations (electronic supplementary material, table SH1): among large populations, density dependence increased the cumulative odds of extinction 5.0-fold for high-diversity populations and 3.6-fold for low-diversity populations (figure 3ci,ii, exponentiated effect sizes, electronic supplementary material, table SH1). Among small populations, density dependence increased the cumulative odds of extinction 1.5-fold and 2.0-fold for high- and low-diversity populations, respectively (figure 3ciii,iv; electronic supplementary material, table SH1). These patterns were also observed for simulations at different parameter values (electronic supplementary material, SG). Greater initial maladaptation (i.e. a mean initial genotype further from θ) was associated with increased extinction risk across all treatments (figure 3d; electronic supplementary material, table SH1). The effect of maladaptation on extinction was stronger for large populations than small populations (electronic supplementary material, table SH1); in large populations with high diversity, only the most strongly maladapted populations had considerable extinction risk (figure 3di). Density dependence increased extinction risk for all genotypes and only marginally affected the relationship between genotype and maladaptation (i.e. raised the intercept while slightly lowering the slope; figure 3d; electronic supplementary material, table SH1). The combined result was that in large populations with density-independent growth, only populations that were highly maladapted were at risk of going extinct, but under density-dependent growth all populations had substantial extinction risk (figure 3d).

• Download figure
• Open in new tab
• Download PowerPoint

The results for probability of rescue were qualitatively identical to those for probability of extinction but with directions switched (table 1; electronic supplementary material, table SH2 and figure SH1). Likewise, the results were qualitatively similar regardless of the type of rescue (i.e. fitness-based rescue or size-based rescue). For both types, the probability of rescue declined with density dependence and density dependence had stronger effects in large populations than small ones. Density dependence increased mean rescue time only slightly (less than 1.5 generations) for all population treatments except for large populations, for which density dependence increased the time to return to initial size by 9–10 generations (electronic supplementary material, figure SH2). Notably, many small populations were rescued instantaneously, i.e. they achieved three generations either exceeding initial size or $\overline{R}>1$ (or both) either instantly or after one time step (table 1 and figure 3a). These ‘instantaneous rescues' happened primarily in small populations, comprising up to 16% of successful rescues observed in small-population treatments. Both the overall frequency and the proportion of total rescues composing these instantaneous rescues decreased with low genetic diversity and with density dependence. Considered across all populations with large initial size, only four populations achieved instantaneous rescue of either type, composing less than 1% of rescue episodes for all large treatments. Extinctions occurred after both fitness- and size-based rescue (table 1). Extinction after rescue was more common for small populations than large ones, for low-diversity populations than high-diversity ones, and under density dependence than density independence for both size- and fitness-based rescue.

(c) Genetic variance, fixation and fitness

Over the first 15 generations, mean additive genetic variance, $⟨{({\sigma }_{\mathrm{a}}^2)}_t⟩$, declined for all treatments (figure 4). Surviving density-dependent populations lost more genetic diversity than density-independent populations (figure 4, solid black versus solid purple lines); this effect was greater in large populations. Across treatments and throughout time, populations that went extinct had lower mean genetic variance than populations that survived (figure 4, solid versus dashed lines). This loss of genetic variation was associated with fixation for both positive (adaptive) and negative (maladaptive) alleles (electronic supplementary material, figure SI1). Populations that went extinct typically had on average more loci at fixation for either allele than surviving populations, although in large populations this difference did not appear until the fifth generation, when populations approached minimum size (electronic supplementary material, figure SI1). Differences in the trajectories over time of genotype and intrinsic fitness between density-dependent and density-independent populations conditioned on survival or extinction were slight (electronic supplementary material, figures SI2 and SI3).

• Download figure
• Open in new tab
• Download PowerPoint

(d) The extinction vortex

Extinct populations in all treatments within the 15-generation simulations had negative mean population growth rate ($⟨r_{\tau }⟩$, figure 5a) and non-zero mean loss of genetic diversity ($⟨{\nu }_{\tau }⟩$, figure 5b) in each time step. Declines in population growth rates and rate of loss of genetic diversity accelerated, as expected for an extinction vortex [25,34] in the three generations prior to extinction (τ ≤ 3), with populations losing an average of 17% of their standing additive genetic variance in the final generation preceding extinction. By contrast, the rate of adaptation ($1-⟨k_{\tau }⟩$) decreased steadily over time in all treatments rather than accelerating immediately prior to extinction (figure 5c). In the three generations prior to extinction, trends in density-dependent and density-independent populations were similar. Four or more generations prior to extinction (τ > 3), large density-dependent populations had faster population decline (i.e. more negative $⟨r_{\tau }⟩$), but recovering growth rate compared to large density-independent populations, for which the rate of population decline was constant over time (figure 5a). By contrast, density-dependent and density-independent populations had similar rates of loss of genetic diversity (all differences in $⟨{\nu }_{\tau }⟩$ between density-dependent and density-independent treatments within 0.02, figure 5b). In the generation prior to extinction, populations in all treatments had a mean adaptation rate of approximately 0 (i.e. $⟨k_{\tau =1}⟩\approx 1$ and $⟨{\overline{g}}_{\tau =1}⟩\approx ⟨{\overline{g}}_{\tau =2}⟩$). Non-normalized state variables conditioned on generation of extinction revealed qualitatively similar results (electronic supplementary material, SJ).

• Download figure
• Open in new tab
• Download PowerPoint

(a) Interactive effects of density dependence and genetic erosion

We found an important interaction between size and density dependence not noted in previous models of rescue. Theory [1] and empirical studies [3–5] predict that in density-independent settings, initial population size confers an advantage to populations experiencing environmental change; our model agrees with this finding (figure 3 and table 1). However, we found that density-dependent growth considerably reduced the advantage that initial size conferred (figure 3; electronic supplementary material, table SH1). Upon environmental change, density-dependent populations at or near the carrying capacity will experience more rapid decline and reach lower densities than density-independent models predict. Natural populations and conservation targets may be below carrying capacity owing to pressures such as predation [35] or inbreeding depression [36]. While populations considerably below carrying capacity are at a disadvantage owing to increased demographic stochasticity, their growth rates will be less affected by density dependence. Thus, populations considered at low risk by classic, density-independent models may face a non-trivial threat of extinction from environmental change under more realistic density-dependent population growth.

Adaptation to environmental change was accompanied by a loss of genetic diversity and potential for adaptation not captured in most models (but see [23,24]). Quantitative genetic models often assume populations are large for analytical tractability [10]. Under this assumption, fixation and loss of alleles are rare as there is only marginal sampling variation and thus minimal drift. However, population decline is by definition part of the rescue process. Our model demonstrates that additive genetic variance is lost even in populations that survive the environmental change (figure 4). This loss of genetic diversity was associated with slowing rates of adaptation, contributing to prolonged time until rescue. This loss was partly driven by fixation of both adaptive and maladaptive alleles, typically occurring when populations were at their smallest size (electronic supplementary material, figure SI1), and thus is likely due to both selection following environmental change and genetic drift, which could also occur to a lesser extent in the absence of environmental change. This matches theoretical results [23] as well as empirical results in a seed beetle adapting to a novel host plant [37]. Likewise, laboratory studies with another beetle, Tribolium castaneum, adapting to novel conditions found that accumulating genetic load likely slows the rate of adaptation [38] and that more genetic diversity likely assists population recovery through higher rates of recombination [4]. We therefore recommend that evolutionary rescue models include genetic erosion where appropriate, as environmental stress [39] and inbreeding due to bottleneck events [11] are known to erode standing genetic variation [12].

Our simulation results also showed a surprising interaction between density dependence and genetic erosion. Under density dependence, the effects of genetic erosion were negative, although small, and did not emerge until nearly the end of simulations (figure 2). Conversely, under density independence, populations that experienced genetic erosion had higher growth rates than populations with genetic diversity held constant. In total, this meant that the effect of density dependence on final population size was considerably larger under genetic erosion than in the analytical model, where genetic diversity was held constant. While theory primarily focuses on the effects of genetic and phenotypic diversity on the rate of adaptation, there is a second effect evident in equation (2.1) that is often ignored in analysis of rescue: a ‘variance load’ [16,40] where mean population fitness is decreased by the increased mean phenotypic distance from the optimum. Our results can be explained by noting that adaptation in response to selection reduces phenotypic variance and thus reduces variance load. In populations reaching low densities owing to density dependence, the positive effect of reduced variance load on population growth is counteracted by drift (including the fixation of maladaptive alleles; electronic supplementary material, figure SI1B) and increased extinctions (figure 3), such that genetic erosion has a net negative effect on expected population size. However, under density independence, populations are less likely to reach low densities and experience drift, meaning the net effect of reduced phenotypic variance is positive, as it is associated with adaptation. This further highlights the need to simultaneously include negative density dependence and effects of stochasticity at small densities in models of rescue.

(b) How an extinction vortex impedes evolutionary rescue

Challenging environments can draw populations into an extinction vortex [25]. Although previous analyses have hypothesized that reduced genetic diversity associated with declining population size may reduce adaptive potential [18,21], to our knowledge our analysis is the first to demonstrate the theoretical existence of positive feedback loops connecting declining genetic diversity, reduced rates of adaptation and reduced population growth that characterize this vortex ([25], figures 1 and 5 in the present study) in the context of evolutionary rescue. Gilpin & Soulé's original conception of the extinction vortex [25] did not assume a population is severely maladapted to its environment, only that small density increases susceptibility to drift, inbreeding, and demographic stochasticity. By contrast, evolutionary rescue is defined by populations facing severe maladaptation that pushes mean fitness below the rate of replacement. Thus, rather than an extinction vortex arising by a population slowly becoming less fit owing to drift, a population experiencing environmental change immediately becomes less fit and also faces an extinction vortex via drift and stochasticity that further reduces its ability to sufficiently adapt to the new environment (figure 5b,c). Negative density dependence, rather than changing the demographic and genetic dynamics in the generations immediately preceding extinction, ‘widens' the opening to the vortex as populations are more likely to quickly decline to small size several generations before extinction (figure 5a). For managers of small, at-risk populations that are susceptible to sudden environmental change and extinction, intervention strategies such as translocation may be needed to pull populations out of the self-reinforcing demographic and genetic processes of the vortex [5,41].

(c) Rapid rescue in small populations

Upon environmental change, some small populations already had a large proportion of phenotypes adapted to the environmental change owing to sampling variation, allowing them to be rescued very quickly (figure 3d and table 1). By contrast, large populations were rarely rescued immediately (table 1) because the average initial phenotype was less variable across populations (figure 3d). These ‘rapid rescues' appear to be a polygenic analogue of the phenomenon demonstrated in a haploid single-locus model [22], where small populations may, under certain conditions, have an over-representation of individuals with genotypes adapted to the novel environment. In our quantitative trait model where adaptive and maladaptive alleles have equal effect on the genotype, small populations are just as likely to have an over-representation of especially maladapted genotypes as they are to have an over-representation of favoured genotypes (figure 3d). Another interpretation of this result is the fact that a sizeable proportion of small populations were rescued instantly owing to a ‘lucky’ initial genetic state (table 1) and underscores the difficulties that small populations without this ‘luck' have in adapting to novel environments.

(d) Operational use of current rescue definitions

Current definitions of rescue are difficult to operationalize in stochastic settings. Rescued is often understood to mean ‘not extinct', but this definition is lacking, as our analysis shows that even on relatively short timescales (50 generations in our simulations), populations can show sensible indications of ‘rescue' and still go extinct (table 1). Likewise, in conservation or management settings, quantifying fitness to determine when mean fitness (intrinsic or realized) exceeds unity may not be feasible [6]. Genetic or demographic stochasticity (to say nothing of environmental stochasticity or other sources of variation) complicates this further, as some natural benchmarks such as population increase may easily be undone, causing populations to revert to an ‘unrescued' state. We defined rescue with two heuristic criteria that provided qualitatively similar (although not identical) results (table 1), but they were not guaranteed to be simultaneously met in the same population. Having more readily operationalizable definitions of rescue could prove to be important both for empiricists testing theoretical predictions and for managers and conservation practitioners seeking to use the theory. For example, rescue criteria can involve surpassing pre-registered quasi-extinction thresholds from population viability analyses [42], although additional work would be necessary to demonstrate sufficient genetic adaptation to the novel conditions.

(e) Caveats and limitations

Our model makes several assumptions for simplicity and tractability that would be fruitful to explore in future work. First, our model shows that density dependence hinders adaptation through genetic erosion, but our model treats density dependence as independent of trait values and fitness. If, alternatively, density dependence acts differently on different trait values, then rates of adaptation will be affected by density dependence [43,44]. For example, competition for scarce resources where more-fit individuals exclude less-fit individuals can hasten the rate of adaptation [45], facilitating rescue [43]. Conversely, such fitness-dependent competition can increase the strength of inbreeding depression [46], which in the case of weak purging can decrease population growth rates [47], possibly hindering rescue. Second, our analysis disregards mutations because we primarily model populations that are small, such that the (re-)introduction of beneficial alleles is relatively rare [48], and the effects of such introductions may be small relative to the effects of other evolutionary processes that happen at small population density, such as drift. Analysis focusing on dynamics at larger population sizes would be improved by including mutations. Third, our model does not consider plasticity, which in the context of evolutionary rescue and environmental tracking is well studied [9,16]. For example, in a model where populations have constant genetic variation in plasticity over time, Chevin & Lande [9] found that the evolution of plasticity can further facilitate rescue [9]. However, our results suggest that effects of genetic erosion (non-constant genetic variation) as populations shrink could reduce the evolutionary potential for plasticity. Finally, our analysis focuses exclusively on the case of a single environmental shift, followed by a constant environment. This is just one of several possible scenarios for environmental change [18], including gradually changing environments and different regimes of temporal environmental variation [16,49]. A natural next step would be to examine the simultaneous dynamics of density dependence and genetic erosion in other regimes of environmental change.

(f) Implications for conservation and management

Our simulations were on short timescales relevant to immediate or near-term conservation concern. However, we note that the reductions in population size or genetic diversity may be important on longer timescales, even after populations have been rescued. An effective population size on the order of hundreds is needed to reduce longer-term extinction risk due to accumulation of deleterious alleles [50–52]. Populations need to spend long durations at large size to offset even brief episodes of drift and genetic erosion [51], like those that might be seen during the ‘trough' or low point in population size following environmental change. Under density dependence, not only are recovered populations held at carrying capacity, potentially precluding this re-establishment of variation, but they also spend more time at smaller densities before being rescued (figures 2 and 3a), leading to greater genetic erosion (figure 4). Indeed, our simulations showed that up to 40% of populations in a treatment that were ‘rescued' under our criteria go extinct within 50 generations, with more extinctions in populations with initially low genetic diversity (table 1). Furthermore, population size and remaining standing genetic variation are important for persisting through future environmental changes. The erosion of genetic variation is similar to a genetic bottleneck that would hinder a population's ability to adapt to additional environmental changes [11,53]. In our simulations, we observed fixation of the positive allele in surviving populations (electronic supplementary material, figure SI1A); in the event of future environmental change favouring another allele, its loss following the first environmental change event would hinder subsequent adaptation [2,20].

Here, we introduce a biologically realistic model for a more complete and applicable understanding of population rescue after environmental change, providing better prediction of potential outcomes and guidance for conservation and management as well as providing testable predictions for empirical validation. We demonstrate that negative density dependence impedes persistence in a novel environment, and that large and well-adapted populations previously predicted to be at minimal extinction risk [1,3,24] do, in fact, face non-trivial risk of extinction under density dependence. For managers overseeing threatened populations, one possible strategy is to increase habitat carrying capacities, for example, by increasing resource availability, which will alleviate reductions in growth rates as populations recover from environmental change. Although our simulated experimental design only featured density dependence as a binary treatment, in natural settings it is a continuous variable that managers may be able to influence. Because expected time to extinction increases exponentially with carrying capacity [17], decreasing strength or degree of competition could drastically improve long-term population fates, particularly for large populations previously perceived as having low extinction risk. A corollary is that, in addition to other predictors of extinction risk under rescue [3,7], the carrying capacity or degree of intraspecific competition should be considered when assessing risks to populations subject to novel habitat alteration, degradation or reduction. Considering how competition and population size may interact to influence adaptation under environmental change could more effectively achieve a management goal of maximizing evolutionary potential [53]. These predictions highlight the effects of density dependence and loss of genetic diversity on populations exposed to severe environmental change, which is crucial in an age of global anthropogenic changes that put populations at unprecedented risk of extinction.