Is Evolution in Response to Extreme Events Good for Population Persistence?

Climate change is predicted to increase the severity of environmental perturbations, including storms and droughts, which act as strong selective agents. These extreme events are often of finite duration (pulse disturbances). Hence, while evolution during an extreme event may be adaptive, the resulting phenotypic changes may become maladaptive when the event ends. Using individual-based models and analytic approximations that fuse quantitative genetics and demography, we explore how heritability and phenotypic variance affect population size and extinction risk in finite populations under an extreme event of fixed duration. Since more evolution leads to greater maladaptation and slower population recovery following an extreme event, greater heritability can increase extinction risk when the extreme event is short. Alternatively, when an extreme event is sufficiently long, heritability often helps a population persist. We also find that when events are severe, the buffering effect of phenotypic variance can outweigh the increased load it causes.


Introduction
Globally, humans are causing substantial environmental perturbations, and these perturbations are likely to become more severe in the future. In particular, climate change is projected to lead to more extreme weather events, including droughts and major storms (Ummenhofer and Meehl 2017). With more severe events comes the potential for dramatic demographic and genetic consequences.
In the process of causing mass mortality, extreme events can act as catalysts of evolutionary change. In fact, there are many examples of rapid evolution in response to extreme events (reviewed in Grant et al. 2017), such as droughts causing finches to evolve larger beaks (Grant and Grant 2014) and hurricanes causing lizards to evolve larger toepads (Donihue et al. 2018). Extreme events are sudden transient changes in the environment, that is, pulse disturbances (Bender et al. 1984). In the ecological literature, pulse disturbances lie at the crossroads of two other forms of environmental change: press perturbations, a sudden long-term change in the environment (Yodzis 1988;Ives and Carpenter 2007;Kéfi et al. 2019), and continuously fluctuating environments (Lande 1993;Ozgul et al. 2012). Despite their transient nature, pulse disturbances can have strong and diverse impacts on ecological systems, ranging from transient ecological dynamics to permanent shifts in ecological states (Holling 1973(Holling , 1996Fox and Gurevitch 2000;Ives and Carpenter 2007;Holt 2008;Hastings et al. 2018). Here we study the most extreme form of a permanent shift-species extinction-examining how extinction risk depends on the length of a pulse event (see also fig. 1 in Holt 2008). Understanding short-term extinction risk after a single disturbance is critical for conservation and management.
Previous work on evolution in changing environments can provide intuition for how evolution might affect extinction risk during or after a pulse disturbance. One focus of the evolutionary rescue literature has been on understanding the consequences of phenotypic change in the context of a sudden long-term or permanent environmental shift (a press perturbation). These studies, some of which account for demographic stochasticity, underline the importance of genetic variance for increasing the probability of rescue (Gomulkiewicz and Holt 1995;reviewed in Alexander et al. 2014;Bell 2017). Similarly, studies of adaptation in fluctuating environments suggest that if an environment is predictable, such as the case of positively autocorrelated fluctuations, genetic variation reduces lag load (Charlesworth 1993;Lande and Shannon 1996;Chevin 2013). This reduction in the lag load leads to higher population per capita growth rates and, consequently, is expected to reduce extinction risk. Whereas when the environment is unpredictable, genetic variance typically increases the lag load. These studies calculated lag load and per capita growth rates in a number of environmental contexts, including a press perturbation, randomly fluctuating environments, and cyclic environments. However, they did not account for demographic stochasticity, the ultimate cause of extinction. As higher long-term growth rates need not imply lower extinction risk (Yahalom and Shnerb 2019;Ellner et al. 2020;Pande et al. 2020), it is unclear whether intuition provided by these earlier studies extends to extinction risk following a pulse disturbance.
To understand extinction risk during and following a pulse disturbance, we introduce an individual-based model that fuses population demography with quantitative genetics. Using a mixture of computational and analytical methods, we examine how phenotypic variation and the heritability of this variation influences population growth, lag load, and extinction risk during and following a pulse perturbation. Moreover, we examine how the magnitude and direction of these effects depend on the duration and intensity of the pulse perturbation.

Model
We use an individual-based model that combines the infinitesimal model of an evolving quantitative trait with density-dependent demography. To gain insights beyond simulating the model, we derive analytical approximations of the probability of extinction using a mixture of deterministic recursion equations and branching process theory (Harris 1964). We assume discrete, nonoverlapping generations. The life cycle starts with viability selection. In each generation t, we impose stabilizing selection around some optimal trait value v t , which is set by the environment in that generation, by making the probability of survival a Gaussian function of phenotype z, with a strength of selection proportional to 1=q 2 . Following viability selection, survivors randomly mate and produce a Poisson number of offspring with mean 2l. The habitat supports at most K individuals. If more than K offspring are produced, only K are randomly chosen. An offspring's breeding value is a draw from a normal distribution centered on the mean of its parents' breeding values and with a fixed segregation variance V 0 , that is, the infinitesimal model (Fisher 1918;Turelli 2017). Its phenotype, z, is this breeding value, g, plus a random environmental component, e, which is a draw from a normal distribution with mean 0 and variance V e . We ignore dominance and epistasis; thus, the phenotypic variance in generation t is the additive genetic variance plus the environmental variance, V p,t p V g,t 1 V e . At equilibrium,V p pV g 1 V e .
Prior to experiencing an extreme event, the populations in the individual-based simulations start with a 100generation burn-in from an initial state, where all N p K individuals have breeding value v p 0 corresponding to the optimal trait value v t during this period ( fig. S1; figs. S1-S4 are available online). To model the extreme event of length t, the optimum trait value increases by Dv and reverts back to its original value after t generations ( fig. 1). Unless otherwise stated, we use the parameter values q 2 p 1, V p p 1, Dv p 2:5, l p 2, and K p 500. These q 2 and V p values represent strong selection and large phenotypic variance relative to those estimated in Turelli (1984). Reducing the strength of selection or phenotypic variance reduces lag load but does not otherwise change our qualitative results. For this set of parameter values, the optimum shift (Dv p 2:5) corresponds to two and a half standard deviations beyond the mean of the trait distribution; consequently, we expect roughly 85% of the population to die in the first generation. We have chosen a high growth rate, l, to reduce extinction from demographic stochasticity in the absence of disturbance. We chose a large enough starting population size and carrying capacity, K p 500, to make approximations reasonable (e.g., normal distribution of traits).

Approximating the Evolutionary and Population Size Dynamics
In appendix A, we derive deterministic approximations for the dynamics of the mean breeding value g t , genetic variance V g,t , and population size N t . If we assume that the distribution of breeding values remains normally distributed, then we show where V t p V g,t 1 V s , with V s p q 2 1 V e the inverse of the effective strength of selection. We also show that the population size in the next generation is where the mean survival probability, s t , is (q 2 =V t ) 1=2 exp[2(v t 2 g t ) 2 =(2V t )]. Regardless of the trait or environmental dynamics, the genetic variance approaches an equilibriumV g p[2V 0 2 V s 1 (4V 2 0 1 12V 0 V s 1 V 2 s ) 1=2 ]=4. In a constant environment with v t p v for all t, the mean breeding value approaches the optimum,ĝ p v, and provided that l 1 1, N 0 s 0 is large enough, andV p pV g 1 V e is small enough, the population size reaches carrying capacity,N p K. Starting from this equilibrium, we can then approximate the response of the population to a shift in the optimum using equations (2)-(4).

Approximating Extinction Risk
We next approximate the probability of extinction using branching processes (Harris 1964). The probability generating function for an individual with average survivorship in generation t is f t (x) p 1 2 s t 1 s t exp[2(1 2 x)l]. Assuming that the effects of density dependence are negligible, we can approximate the probability of extinction by the end of generation T, given an extreme event of length t began in generation 1, as To calculate s t , we assume V g,t pV g and use equation (2) to get g t , which together give s t (eq. [1]).

Demographic Recovery
We first explore extreme events lasting a single generation. To characterize the impact of phenotypic variance and heritability on population size, we compare the demographic response of populations with low or high phenotypic variance,V p , across a range of heritabilities, h 2 pV g = V p . During the event, heritability has no effect on population size ( fig. 1A) but phenotypic variance does ( fig. 2A). A population with higher phenotypic variance has a smaller population size immediately following a low-severity extreme event but a larger population size following a highseverity event. This pattern stems from the dual role of phenotypic variance, in that it both increases variance load and contributes individuals with extreme traits who are able to survive an extreme event. High phenotypic variance therefore reduces both mean fitness within a generation and the variance in fitness across generations-a form of short-term bet hedging that can increase the geometric mean of fitness in the generations during and after the disturbance event.
While heritability has no effect on survival during the event, it has a strong effect on population recovery in subsequent generations. In particular, heritability dampens the growth rate in subsequent generations ( (2) and (4). Parameters: q p 1, l p 2, Dv p 2:5. Red curves: the generation after the event, generally slowing population recovery. For longer events, evolution can help population recovery ( fig. 1B). While evolution still increases maladaptation after the event, the increase it causes in growth rates during the event can more than compensate.

Extinction Risk
When a single-generation extreme event is severe enough, increasing phenotypic variation lowers extinction risk both during and after the event ( fig. 2B). The biological intuition behind this pattern is the same as in figure 2A, where increased variance means that more individuals survive the extreme event. However, at such large population sizes, the extinction risk is essentially zero during a mild event. In other words, while having too much variance leads to considerable reduction in population size when events are mild, it is very unlikely to lead to extinction unless there is extremely high phenotypic variance or if carrying capacity is very low. In the former case, load will cause extinction in the absence of extreme events ( fig. S2).
To isolate the effect of evolution, we next compare extinction risk for populations with the same phenotypic variance but different heritabilities. When the extreme event lasts only one generation ( fig. 3A), heritability increases shortand long-term extinction risk. However, for two-generation events, long-term extinction risk is lowest at intermediate heritabilities ( fig. 3B). For three-generation events (or longer), long-term extinction risk decreases with heritability ( fig. 3C ). These patterns hold for milder (Dv p 2:5) and more severe (Dv p 4:5) extreme events (figs. S3, S4).
While equation (5) gives a good approximation of extinction risk, the function itself is too complex to give us intuition. Next, by writing down the geometric mean fit-ness of a population, we reproduce the general trends in long-term extinction risk but with added clarity for how maladaptation contributes to these outcomes.

Contribution of Lag Load
To better understand how evolution affects the probability of extinction, we can approximate the geometric mean fitness of a population as where P T tp1 (v t 2 g t ) 2 =(2V) is the cumulative lag load. This cumulative lag load over T 1 t generations is where v pV g =V is a measure of evolvability (Charlesworth 1993; see app. C). In the long term, the cumulative lag load is Equation (8) generalizes a result from Chevin (2013), who considered the special case of weak selection and a press perturbation (t → ∞).
Noting that v p h 2V p =(V p 1 q 2 ), we find that equation (8) determines how heritability affects the long-term cumulative lag load ( fig. 3D). When the extreme event lasts only one generation (t p 1), the cumulative lag load simplifies to Hence, increasing heritability increases the cumulative lag load (solid blue curve in fig. 3D), a trend consistent with the extinction probabilities for t p 1 ( fig. 3A). Alternatively, when the extreme event lasts two generations (t p 2), the cumulative lag load becomes Dv 2 V p 1 q 2 and is independent of heritability (solid green curve in fig. 3D). Finally, when the extreme event lasts for more than two generations, the cumulative lag load decreases with heritability (solid red curve in fig. 3D), a trend consis-tent with extinction probabilities decreasing with heritability when t ≥ 3 ( fig. 3C).

Discussion
Although it has long been recognized that evolution may affect a population's response to a changing environment, previous studies have primarily focused on understanding this effect over the long term following a nonreversing environmental shift (a press disturbance) or a continuously fluctuating environment. Here, we are concerned with the shortterm effect of a pulse disturbance on population growth and extinction risk. By allowing pulses to be of any duration, we connect our results with this existing literature while providing new insights into the transient dynamics following a single disturbance. Our results provide two general conclusions about the effect of trait variation and its heritability on population growth and extinction risk during and following a pulse disturbance. First, trait variance, whether heritable or not, is a double-edged sword: it adds a variance load due to stabilizing selection yet also produces individuals with extreme traits who can survive large shifts in the environment. Second, while variance can be useful in the generation of a severe event, if heritable it can induce maladaptation that slows demographic recovery and therefore increases extinction risk in the generations after the event.

Phenotypic Variance
Phenotypic variance, whether heritable or not, can be beneficial or deleterious. A simultaneous reduction in the mean and variance in fitness before and during an extreme event can increase the short-term geometric mean of fitness ( fig. 2). This increase occurs only when disturbances are sufficiently severe. Furthermore, variation in survival rates reduces variation in the total number of offspring produced by the population (Kendall and Fox 2002) and thereby lowers extinction risk (Lloyd-Smith et al. 2005). Prior studies of evolutionary rescue have emphasized the beneficial aspect of genetic variance (Charlesworth 1993;Gomulkiewicz and Holt 1995;Bell and Collins 2008;Alexander et al. 2014;Barfield and Holt 2016)-but not nonheritable phenotypic variance-in rescuing a population from an abrupt shift in environment. Here, by teasing out the effects of heritability and phenotypic variance, we emphasize the costs and benefits of each.

Heritability
Contrary to evolutionary rescue of populations experiencing a press perturbation (Gomulkiewicz and Holt 1995;Barfield and Holt 2016), we find that heritability increases extinction risk for short pulse perturbations. We can gain some intuition for why this is by considering the limiting cases of traits not evolving versus tracking the optimal trait perfectly with a one-generation lag. When the population is adapted to the original environment but does not evolve in response to the extreme event, it experiences a reduction in fitness for the duration t of the extreme event. In contrast, when selection tracks the optimal trait with a onegeneration lag, the population experiences a reduction in fitness only in the first and last generation of the extreme event. Hence, when the extreme event lasts one generation, extinction risk is higher for the evolving populations, and when the extreme event lasts more than two generations, extinction risk is higher for the nonevolving populations. A similar understanding can be gained by adapting a classic population genetic model of allele frequency change with time-varying selection (Felsenstein 1976;Dempster 1955; see appendix D, available online).
In general, the trends in short-term extinction risk are parallel to the lag load predictions (fig. 3). However, they differ in two ways. First, when the extreme event lasts exactly two generations, the nonevolving population experi-ences the reduction in fitness in successive generations while the evolving population experiences this reduction in alternate generations. Hence, the evolving population is slightly less likely to become extinct (see app. B). Second, when a population exhibits an intermediate amount of tracking of the optimum, the variance in survival from year to year is reduced and therefore can lower the overall extinction probability (fig. 3B).
While previous studies of temporally variable selection have focused on large populations in the long term, calculating lag load and growth rates when rare, they provide intuition for our results on short-term extinction risk after a one-time event. A single-generation extreme event functions like a negatively autocorrelated fluctuating environment, where a strong genetic response to selection in one generation is likely maladaptive in the next generation (Charlesworth 1993;Lande and Shannon 1996;Chevin 2013;Benaïm and Schreiber 2019). However, an extreme event lasting three or more generations acts like a positively autocorrelated environment, where evolution tends to reduce maladaptation.

Future Challenges and Directions
Our models include a number of simplifications to both evolutionary and demographic processes. First, we do not model the erosion of genetic variance with decreasing population size (Lande and Barrowclough 1987;Barfield and Holt 2016). Furthermore, we have limited our analysis to truly continuous traits, but different genetic architectures, such as a few loci of large effect, likely will respond differently (Barghi et al. 2020). Second, our model ignores the potential for phenotypic plasticity, which has variable effects on evolution and extinction risk (Kopp and Matuszewski 2014;Lande 2015). Third, we used the simplest possible model for density dependence-the ceiling model-as used in previous evolutionary rescue studies (e.g., Bürger and Lynch 1995). For other models of compensating density dependence, we expect similar results. However, overcompensatory density dependence can result in oscillatory population dynamics for which the timing of the extreme event may play a subtle role.
An important next step will be to understand evolution and extinction risk under repeated extreme events. Extreme events, or catastrophic events, can be characterized by causing abrupt, infrequent, and large reductions in biomass or population size. Hence, prior work on adaptation and persistence using autoregressive processes to model environmental fluctuations (e.g., Charlesworth 1993; Lande and Shannon 1996;Chevin 2013;Benaïm and Schreiber 2019) do not accurately capture the nature of repeated pulsed disturbances of the type considered here. We hope future studies exploring the impact of disturbance regime on evolution and extinction risk will benefit from the detailed understanding, like that provided here, of an evolving population's response to a single extreme event.

Dynamics of the Breeding Value Distribution and Population Size
Let the trait value of an individual be the sum of a genetic component (breeding value) and an environmental component, z p g 1 e. Assume that we start with a population in generation t that has a normal distribution of breeding values, p g ( g, t), with mean g t and variance V g,t .
And assume that each environmental component is independently chosen from a normal distribution, p e (e), with mean 0 and variance V e . The joint distribution of g and e, p g,e ( g, e, t), is then initially multivariate normal with mean ( g t , 0), variances V g,t and V e , and no covariance. Let the probability of survival for an individual with trait value z in generation t be where v t is the optimum trait value in generation t and 1=q 2 is the strength of selection. The joint distribution of g and e following viability selection is p 0 g,e (g, e, t) p s(z, t)p g,e ( g, e, t) where s(z, t)p g,e (g, e, t)dgde is the expected fraction of the population that survives in generation t (i.e., the population mean survival probability), with V t p V g,t 1 V s and V s p q 2 1 V e the inverse of the effective strength of selection. Integrating over environmental effects then gives the distribution of breeding values among the survivors, which is normal with mean g t (1 2 V g,t =V t ) 1 v t V g,t =V t and variance V g,t (1 2 V g,t =V t ). The mean breeding value is thus shifted toward v t with a weight of V g,t =V t , and the genetic variance has been reduced by this fraction. We next assume that the breeding value is determined by a large number of small effect loci, such that the distribution of breeding values among siblings, p g,sibs ( gjg mid ), is normal with a mean equal to the midpoint of the parental breeding values, g mid , and a variance, V 0 , that does not depend on the parental genotypes or trait values (i.e., the infinitesimal model ;Fisher 1918;Barton et al. 2017). The distribution of breeding values among the offspring is then which is normal with mean and variance That is, the mean breeding value remains constant through reproduction, while the variance before reproduction is first halved (because of essentially blending inheritance between the parents) and then increased by segregation, V 0 . So we see that given that the initial distribution of breeding values is normal, with Gaussian selection the breeding value distribution remains normal, allowing us to track the entire distribution of breeding values (and therefore phenotypes) across generations by keeping track of only its mean and variance. The variance dynamics are independent of the environment (v t ) and the breeding values; solving equation (A7) gives the genetic variance in generation any t. This expression is rather complicated (see Mathematica file, available online); however, it reaches an equilibrium Holding genetic variance constant at its equilibrium (which is reasonable, given that the variance is not expected to change with the environment or breeding values), in a constant environment, v t p v, the mean breeding value in any generation t is found by solving equation (A6), implying a geometric approach to g p v that becomes faster withV g =(V s 1V g ). We assume that each individual that survives viability selection produces l offspring and that if more than K offspring are produced, then K of these are randomly chosen to start the next generation. If the population size in generation t was N t , then the population size in generation t 1 1 is expected to be N t11 p min(N t s t l, K): ðA10Þ

APPENDIX B Extinction Risk in One-and Two-Generation Events
In this appendix, we examine the effect of long-term extinction risk when populations are either not evolving or are perfectly tracking, with a one-generation lag behind the optimal trait value. Let s o and s m be the survivorship of individuals with the optimal trait or the maladaptive trait.

APPENDIX C Cumulative Lag Load
Here we show how to derive equation (7). Our goal is to develop a formula for the cumulative squared displacement, C(t, T) p P T tp1 (v t 2 g t ) 2 , given event length t. First, note that equation (2) implies that with constant genetic variance V g , the mean trait displacement in the next generation is g t11 2 v t11 p ( g t 2 v t )(1 2 v), where v pV g =V is a measure of evolvability. Thus, if the optimum is fixed at some arbitrary value for t generations, then the displacement in generation t, d t p g t 2 v t , is d t p d 0 (1 2 v) t , and the cumulative squared displacement over those t generations is d 2 0 P t21 tp0 (1 2 v) 2t . If the optimum then reverts to its original value for a further T 2 t 1 0 generations, then the initial displacement is d 0 (1 2 v) t 2 d 0 , and the cumulative squared displacement over this period is d 2 0 [(1 2 v) t 2 1] 2 P T21 tpt (1 2 v) 2(t2t) . Combining these two sums, we get (1 2 v) 2(t2t) : Evolution in Extreme Events 51