A Figures

# Predicting evolutionary rescue via evolving plasticity in stochastic environments

## Abstract

Phenotypic plasticity and its evolution may help evolutionary rescue (ER) in a novel and stressful environment, especially if environmental novelty reveals cryptic genetic variation that enables evolution of increased plasticity. However, the environmental stochasticity ubiquitous in natural systems may alter these predictions because high plasticity may amplify phenotype-environment mismatches. Although previous studies have highlighted this potential detrimental effect of plasticity in stochastic environments, they have not investigated how it affects extinction risk in the context of ER and with evolving plasticity. We investigate this question here by integrating stochastic demography with quantitative genetic theory in a model with simultaneous change in the mean and predictability (temporal autocorrelation) of the environment. We develop an approximate prediction of long-term persistence under the new pattern of environmental fluctuations, and compare it with numerical simulations for short- and long-term extinction risk. We find that reduced predictability increases extinction risk and reduces persistence because it increases stochastic load during rescue. This understanding of how stochastic demography, phenotypic plasticity, and evolution interact when evolution acts on cryptic genetic variation revealed in a novel environment can inform expectations for invasions, extinctions, or the emergence of chemical resistance in pests.

###### keywords:
Evolutionary Rescue, Phenotypic Plasticity, Baldwin Effect, Environmental Predictability, Stochastic Demography, Environmental Stochasticity, Cryptic Genetic Variation
\cortext

[cor1]Corresponding author

## 1 Introduction

Abrupt environmental change beyond species’ tolerance boundaries occurs both naturally and due to human-driven global change ((1)). Change affecting an entire population (or one unable to disperse) leaves two possibilities for persistence: adapt or acclimate, that is, genetic evolution or phenotypic plasticity ((2)). Adaptive responses after a shift in the environment can prevent extinction if there is sufficient additive genetic variation ((3)). Such evolutionary rescue (ER) takes time, however, and a declining population may go extinct before evolutionary response leads to positive growth and recovery of population size ((4)). Response via phenotypic plasticity may be faster, while also permitting survival in novel environments and time for further evolution.

Evolution and plasticity thus inevitably interact. On one hand, perfectly adaptive plasticity prevents selection on fixed genetic characters ((5)) and more generally may reduce the strength of selection on a trait in predictable environments. On the other hand, partially adaptive plasticity, simply by increasing survival in the new environment, results in more time for selection, and thus evolution, before extinction ((6); (7)). Furthermore, plasticity may itself evolve if it varies genetically (GxE interaction, ((8))). Plasticity, when quantified as the slope of a linear reaction norm, can theoretically evolve to become transiently higher in new environments ((9)). This greater plasticity among surviving lineages requires that the environmental shift causes increased additive genetic variance () of the trait under stabilizing selection due to plasticity (i.e., that stress reveals “cryptic” , which occurs in some cases (reviewed by (10); (11); (12)) although the opposite pattern is also frequently found). Furthermore, empirical observations of heightened plasticity in lineages surviving anthropogenic disturbances like climate change ((13); (14)) or transcontinental introductions ((15)) agree with suggestions from deterministic theory that plasticity facilitates ER ((16)).

Both the evolution of plasticity ((17); (18)) and extinction risk ((19); (20)) depend on environmental variation. For plasticity, variation in the environment favours evolution of plasticity if an environmental cue reliably correlates with the environment that imposes selection ((17); (18)). More precisely, the optimal level of developmental plasticity matches the correlation between the environment of development and that of selection (i.e., environmental predictability), with any mismatch reducing the expected long-term fitness ((18)). For extinction risk, long-run population growth determines long-term persistence, and declines with increasing variance in environmental fluctuations in the growth rate ((19)). When such fluctuations are positively autocorrelated, they increase extinction risk by allowing for many successive generations of negative growth ((21)). In contrast, autocorrelation in phenotypic selection might decrease extinction risk, because it allows closer evolutionary tracking of an optimum phenotype, thus increasing the mean population growth rate ((22)). Importantly, the pattern of environmental stochasticity affects not only the mean, but the whole distribution of population sizes. In the context of ER, this implies that many populations may go extinct, even when the expected population does not ((19); (23)).

These separate influences of stochasticity on plasticity’s evolution and on extinction suggest the potential for stochasticity to reduce, or possibly reverse, the adaptive role of plasticity during ER. For instance, if predictability is low and plasticity is high, environmental variation in mean fitness will be large (because excess plasticity causes overshoots of the optimum; ((24))). Such excessive plasticity can lead to extinction ((25)). For example, if the environment undergoes an abrupt change in predictability, which can happen if its temporal autocorrelation rapidly changes, phenotypic plasticity might become transiently maladaptive, which would not only reduce the expected fitness, but also increase the variance in population sizes across replicates, further increasing extinction risk (as shown without plasticity by Ashander and Chevin in prep). Therefore, considering environmental stochasticity is necessary to understand the conditions under which evolving plasticity enhances or impedes ER. Yet previous analytical treatments (e.g., (3); (16)) have neglected stochastic effects on ER, which has rarely been studied outside of simulation models (e.g., (26)). Furthermore, for plasticity to evolve at all requires GxE interactions, which with linear reaction norms may lead to higher phenotypic variance in novel environments ((18); (7); (9)). Large phenotypic variation in a new environment, for a trait under stabilizing selection, results in standing variance load that reduces population growth, which may prevent long-term persistence and thus impede ER (see ((16)), Fig. 1c at large ) but whether these effects occur in stochastic environments is unknown.

Here, we investigate whether and how stochastic environmental fluctuations, and the variance load induced by expression of cryptic genetic variance in a novel environment, constrain evolutionary rescue with evolving plasticity. To do so, we integrate quantitative genetic theory on evolution of plasticity with stochastic demography. Modelling a large shift in the mean optimum trait, to a value outside the previous range of temporal environmental variation, combined with a change in the environmental predictability of fluctuations in this optimum, we develop an approximation for the population growth rate after the mean trait reaches a stationary distribution around the expected optimum. We also examine, using simulations of the underlying model, the risk of quasi-extinction both in the short term and overall. The approximations predicts long-term persistence in the new environment, quantifying the eco-evolutionary dynamics that emerge with evolving plasticity when a major detrimental environmental shift is combined with random environmental fluctuations. We find that for ER to occur in these conditions the environmental predictability after the shift must be above a critical level.

## 2 Materials and methods

### 2.1 Reaction norm, phenotypic selection, and population dynamics

We assume random mating in a closed population with discrete generations and environmental stochasticity that is “coarse-grained” so that every individual in a generation experiences the same environment. The environment both determines an optimal value for a primary trait , and cues a plastic response from that trait. We assume linear dependence of the optimal trait on the selecting environment , so (where the environmental sensitivity of selection defines the change in the optimum phenotype for a unit change in environment ). We model the genotypic reaction norm (i.e., plasticity), a linear response in the trait to the environmental cue with slope and intercept , such that the phenotype of an individual is . Here the residual environmental variation is independent of the macro environment and has mean zero and variance ((9)). Our model applies to irreversible (non-labile) forms of plasticity such as developmentally-plastic traits.

We model the reaction norm intercept and slope as quantitative traits with means and and additive genetic variances and , respectively ((18); (9)), so the intercept represents the breeding value in a reference environment, . In addition, we assume that the population has evolved in a range of environments centred around zero, so that phenotypic variance in the reference environment is minimal. Then with linear reaction norms as here the slope and intercept have zero additive genetic covariance ((9)), and the additive genetic variance of the expressed trait increases quadratically away from the reference environment (grey band in Figure 1(a)), which implies strong increases in heritability away from the reference environment. The mean and variance of the expressed trait value before selection are

 ¯z(t)=¯a(t)+¯b(t)εc(t) (1a) σ2z(εc(t))=σ2a+σ2bε2c(t)+σ2e, (1b)

which assumes the reference additive genetic variances are constant in time. The expressed trait and the slope have covariance and so, with large , direct selection on the trait results in stronger selection on reaction norm slope ((9)). (This does not hold with an alternative assumption where variance decreases away from the reference environment; see F.)

In each generation, we assume stabilizing selection for an optimum value, which gives an expression for the Malthusian growth rate. Given a fitness function of width and optimum defined above, absolute fitness is . Averaging over the (normal) distribution of phenotypes within a generation, mean maladaptation of the trait from the optimum, , drives mean absolute fitness ((27); (9)) , where is the strength of stabilizing selection experienced by the population. Because the strength of selection depends on the phenotypic variance, it changes with the environmental shift according to eq. (1b); however, in the new environment it is approximately constant (if selection is weak as we assume here) with value as we show below in 2.3.2 (assuming small covariance between reaction norm intercept and environment, see ((28))). Trait change in a generation is the product of the additive genetic variance-covariance matrix and the selection gradient , i.e., , with the vector of reaction norm traits ; following Lande ((27)), the gradient of log mean fitness with respect to gives (see C). The dynamics for population size (assuming very weak density dependent regulation, or a form of density-dependence, e.g., large exponent in a theta-logistic model, where growth trajectories during rescue are similar; see ((16))) are then , and the population’s Malthusian growth rate is

 r(t)=rmax+log√Sδω2−Sδ2x(t)2 (2)

Two terms reduce the growth rate from its maximum . They correspond to “loads” well-known in evolutionary biology. The first is standing variance load due to phenotypic variability in any generation, the second is the lag load ((29)) due to maladaptation . The latter can be further decomposed into a “stochastic load” caused by fluctuations in the optimum, causing mismatch (with average zero) between the mean trait and the optimum and a deterministic “shift load” caused by mismatch of the mean trait (after accounting for fluctuations) relative to the mean optimum.

### 2.2 Environmental stochasticity, shift in the optimum, and resulting dynamics

We consider an autocorrelated stochastic environment, where noise arises from a stationary process with variance and autocorrelation function (here, is the characteristic autocorrelation time and implies the process is white noise). The environment of development determines the trait with a delay of less than one generation, so determines both the environment of selection and that of development , which is the environmental cue ((9)). (In general, environmental variables acting as cue might differ from the environments causing selection, but this is beyond our present scope.) In the reference environment, the correlation between the cuing environment and the selecting environment represents the environmental predictability of selection (or cue reliability; ((30); (25))), which in our case is the autocorrelation over time , given by .

A single discrete shift in the environment changes the mean of both the environment of selection and the environmental cue by the same amount , and also changes the environmental predictability from to . Our analysis assumes that, even accounting for stochastic variance in the environment, the new optimum is very different from the previous one, i.e., . Under this assumption, rescue occurs over two phases at different timescales ((9)). Over the rapid Phase 1, the maladaptation of the mean phenotype is reduced to near zero by evolution of increased plasticity (mean reaction norm slope increases from solid black line to grey line in Figure 1(a), rapidly as in Figure 1(b) left of the dotted lines). During the slower Phase 2, the mean reaction norm height and slope evolve to approximately their optimal values (mean reaction norm slope decreases from solid grey line to dashed black line in Figure 1(a), slowly as in Figure 1(b) right of the dotted lines). Lande ((9)) showed that the rapid increase of plasticity is controlled by the proportion of additive genetic variance in the new environment due to variance in reaction norm slopes, , where assuming a large shift implies is near . In the new environment, the maladapted population declines at first (Figure 1(c)) because it has a negative expected growth rate, which increases over time due to adaptive evolution (eq. 2; Figure 1(d)). Because Phase 1 occurs much faster than Phase 2, the growth rate first increases rapidly during this phase, then effectively stabilizes.

### 2.3 Simulations and analysis

We quantify the risk of extinction before rescue (i.e., before the end of Phase 1) by using simulations to compute the proportion of trajectories that experience quasi-extinction. Additionally, we quantify two components of extinction risk, which reflect the action of differing eco-evolutionary processes. First, the short-term risk reflects both temporally stochastic environment and deterministic effects of reduced shift load during ER. We compute it by using simulated quasi-extinction during the initial population decline before the end of Phase 1. Second, the long-run growth rate reflects effects of stochastic load and standing variance load at stationarity; a negative long-run growth rate implies eventual extinction. We analytically approximate the long-run growth rate at the end of Phase 1, and also we compute it from simulations for comparison.

#### Simulations of short-term and total extinction risk

For all simulations of extinction risk, we drew initial conditions from a stationary distribution of reaction norm intercept and slope generated after 15,000 generations in an environment with mean and predictability (where theory predicts and slope ; see ((18); (9))). Using mean fitness as growth rate, we compute population dynamics as . We also track the reaction norm parameters , and mean fitness , with trait change in a generation given by the product of genetic variance (assumed constant, at an equilibrium between mutation and stabilizing selection) and the selection gradient (see 2.1). We computed means of these quantities at each generation across 250 replicate simulations (except in Figure 1, where we used 50 replicates to ease visualization). We also performed simulations where we relaxed the assumption of constant variance, assuming instead that variance reaches an equilibrium at each population size due to mutation-selection-drift balance, (under a modified stochastic house-of-cards approximation ((23)); see Figure S5). All numerical simulations and plotting of analytical predictions were performed in R ((31); (32); (33)); further details are in E.

We define quasi-extinction probability as the proportion of population trajectories that fall below a critical population size before time . From the simulations, we computed quasi-extinction probability for the short-term, at a time before the end of Phase 1 (, with as half the characteristic timescale of Phase 2). Also from the simulations, we computed the total probability of quasi-extinction, over the whole time period into Phase 2 (, see Figure 1(c); we defined as the characteristic timescale of Phase 2 plus 500 generations).

#### Analysis of the approximate long-run growth rate

A positive long-run growth rate (asymptote of the black line in Figure 1(d)) is necessary for long-term persistence after rescue. The long-run growth rate depends on the expectations of the standing variance and lag loads, taken over stochastic fluctuations, of log mean absolute fitness ((19); (23); (22)). An analytical formula for growth rate follows from two main assumptions, under which we compute the expected loads (derived in B). First, if stabilizing selection is weak, by taking an expectation over stochastic fluctuations the variance load caused by increased phenotypic variance is approximately constant, with value , where , is the average selection strength in the new environment. Also assuming weak stabilizing selection, the expected lag load is approximately and is the only quantity that varies in time. Thus, the lag load determines temporal variation in the long-run growth rate, which is .

A second assumption is that fluctuations in maladaptation achieve stationarity at the end of Phase 1, which permits us to compute the variance in maladaptation . Because the lag load can be decomposed into the mean and variance in maladaptation, and the mean maladaptation goes to zero by the end of Phase 1 (see C), the long-run growth rate depends only on the variance in maladaptation , and is . The variance in maladaptation affects not only the spread in trajectories of growth rates, and thus population size, but also long-term persistence, based on the long-run growth rate. After Phase 1, change in reaction norm parameters is slow which provides some justification for the assumption, which yields an explicit formula for the variance in maladaptation, and thus the growth rate. The formula depends on the value of plasticity is at the end of Phase 1 (with value ), as well as characteristics of the novel environment (see full derivation in D). Persistence occurs when

 Missing or unrecognized delimiter for \right (3)

Because depends on environmental predictability in the new environment, the size of shift in mean , and (through selection strength ) variance in plasticity , eq. (3) defines critical levels of these that are necessary for a population to persist. For comparison to the analytically-predicted critical parameter values ( eq. 3), we also computed (from the simulations described above) the stochastic population growth rate over the cusp between Phases 1 and 2 (i.e., between to Figure 1(d)) from its maximum likelihood estimator, ((34)).

## 3 Results

We found that with evolving plasticity, the combination of a major environmental shift with stochastic fluctuations in the environment alters the eco-evolutionary outcome, as compared to the effects of each of these factors in isolation. In particular, our analysis reveals the importance of environmental predictability for ER with evolving plasticity.

### 3.1 Environmental predictability is critical to ER with evolving plasticity

Evolutionary rescue following a shift in the mean optimal trait requires a critical level of final environmental predictability. Predictability below this critical level both reduces long-term persistence (Figure 2(a-b) and increases extinction risk in the short term (Figure 2(c,d) . Furthermore, this decreased predictability strongly increases total extinction risk across a range of initial genetic variances in plasticity (Figure 2(e)) and sizes of the environmental shift (Figure 2(f)). None of these effects of predictability can be understood from deterministic models, which predict inaccurate trajectories over rescue in the presence of stochasticity (see Figure S4).

Increases in stochastic load, due to increased plasticity during ER, cause this constraint by reducing the long-run growth rate. For any fixed shift size ( Figure 2(b)) or additive variance in plasticity ( Figure 2(a)), decreasing the final environmental predictability increases the stochastic load because plasticity in excess of predictability causes the mean trait to overshoot the mean optimum (see D and Figure S2) resulting in larger mismatch variance and hence decreased expected growth rate after Phase 1. It is important to note that a reduction in predictability (i.e., temporal autocorrelation in the optimum) is expected to increase stochastic load even without plasticity, because it will decrease adaptive tracking ((22)). In our case, however, evolved increases in plasticity causes much higher stochastic load (often more than four times greater, Figure S2). In simulations where the genetic variance changes with population size due to drift, there is still a critical predictability (see Figure S6) but it increases much faster with shift size and the effect of standing variance load disappears.

### 3.2 Effects of genetic variance in plasticity depend on initial plasticity

In populations with low initial plasticity, determined by the environmental predictability with which a lineage has evolved, increasing genetic variance in plasticity can greatly reduce short-term risk of quasi-extinction (for Figure 2(c,d) left column, for all below about 0.025). Effects on the total risk of extinction are similar (Figure 2(c), left column).

On the other hand, in populations with high initial plasticity, genetic variance in reaction norm slope does not affect extinction very much in the short term ( for below 0.025 with high predictability in Figure 2(c,e) right column). Such populations generally have much lower risk of short-term quasi-extinction consistent with earlier results of Chevin and Lande ((16), their Fig. 2) on the effect of initial relative plasticity on (deterministic) extinction. However, in these populations when predictability following the shift is intermediate to low, increasing additive variance in reaction norm slopes decreases long-run growth rate (positively sloped solid lines in Figure 2(a) indicate increased variance moves from ‘+’ to ‘-’, with similar increase in the total risk extinction (Figure 2(c), right column).

In low-plasticity populations transitioning to high predictability environments, the effects of stochastic load are low and the benefits from increased plasticity during ER can outweigh the negative effects of the standing variance load. In contrast, for populations with initially high plasticity the change in plasticity during ER is small and the effects of standing variance load and stochastic load dominate. This can be seen by comparing the analytical persistence threshold (Figure 2(a-b)) to the total extinction probability (Figure 2(e,f)). For high initial plasticity (right columns) the analytical prediction matches the total extinction probability, but this is not the case for low initial plasticity (left columns). The equation predicts the same constraint in both low- and high-plasticity populations (solid lines have similar shape in both columns of Figure 2(a-b)). Numerical simulations of long-run growth rate agree (heatmap and dotted line have same shape in both columns of Figure 2(a-b)). This occurs because eq. (3) depends on plasticity at the end of Phase 1, which is influenced very little by initial plasticity. The condition of positive long-term growth, however, is only necessary for ER, not sufficient, and is sensitive only to effects of variance load and stochastic load. Total extinction risk reflects deterministic reduction of lag load due to increase plasticity during ER; such increases are more important for low-plasticity populations, which incur greater maladaptation initially.

### 3.3 Analytical prediction of the growth rate performs well

Across most parameter values agreement between eq. (3) and simulations is strong (compare solid and dotted lines within Figures 2(a-b)). The exception is high initial plasticity and small environmental shift (solid and dotted lines mismatch for low predictability in Figure 2(b)), where changes in are driven by small shifts and large additive variance in plasticity in the reference environment (see G).

## 4 Discussion

Our analytical results and simulations reveal the eco-evolutionary dynamics that emerge with evolving plasticity when a major detrimental environmental shift is combined with random environmental fluctuations. We find that whether evolving plasticity will enable ER depends on environmental predictability in the new environment. If predictability is moderate to low after an environmental shift, the transient evolution of high plasticity that occurs in the new environment causes a large stochastic load that reduces the likelihood of ER. Even without plasticity, environmental predictability affects the stochastic load (lower autocorrelation reduces adaptive tracking of the optimum by genetic evolution; ((22))), and thus the probability of evolutionary rescue in a stressful new environment (Ashander and Chevin in prep). When ER causes increased plasticity, these effects are stronger (see Figure S2), because frequent mismatches caused by excess plasticity result in large variance in growth rate and negative population growth (Figure S4(d)), even after mean maladaptation has reduced to zero (Figure S4(a)). This parallels findings that non-evolving plasticity can amplify fluctuations in population mean fitness and thus growth rate without an environmental shift ((25); (24); (35)); here, we demonstrate that this process also constrains ER. On the other hand, if predictability is high in the new environment, ER is relatively likely. Our findings accord with other theory that suggests irreversible developmental plasticity is not useful in low predictability environments, which might instead favour reversible plasticity ((36)). In addition, we find that positive effects of increased genetic variance in plasticity for ER in a stochastic environments are limited to situations where lineages with low evolved plasticity experiencing shifts to a more predictable environment.

Due to the trade-off between short-term adaptive benefits and long-term stochastic and standing variance loads, large genetic variance in plasticity can sometimes decrease the chance of ER for lineages where plasticity is initially high that experience shifts to environments with low predictability. Chevin and Lande ((16), their Fig. 1c at large ) noted the effect of the standing variance load, but focused on deterministic environments that do not include random noise. We extend this theory to noisy environments and demonstrate, for low predictability, the effects of genetic variance in plasticity: it lessens exposure to small population size in the short term, but may cost increased variance load that slightly reduces long-term persistence (in population with already-high plasticity that shift to low-predictability environments). (Note, however, that in our model the genetic variance in plasticity is fixed, and so the model cannot produce any direct selection to reduce this load.) These effects of genetic variance, however, are weak compared to the constraint imposed by predictability.

Our findings extend deterministic theories on ER (e.g., (3); (16)) by directly quantifying the effect of environmental stochasticity with evolving plasticity on persistence. Although models of ER on quantitative traits have not typically accounted for environmental stochasticity (but see (23); (26)), demographic stochasticity has been shown to affect ER by de novo mutation or standing variation at a single gene (e.g., (37)) because population trajectories depends on births and deaths during the initial period when advantageous genes were rare. Environmental stochasticity, our focus here, is arguably more important than demographic stochasticity because it operates with equal strength at all population sizes ((20)), and affects the population size distribution during ER or with fluctuating selection on quantitative traits (Ashander and Chevin in prep). Stochasticity’s effect on ER with evolving plasticity may be especially strong, as mismatches of increased plasticity with predictability both increase short-term extinction risk and reduce long-run growth. We also showed that lowered predictability can cause the high plasticity that evolved during ER to eventually be maladaptive, unless the new environment is highly predictable.

### 4.1 Assumptions and caveats

To obtain analytical approximations that predict evolutionary trajectories we used three main assumptions. We assumed first, that baseline additive variances in reaction norm parameters remain constant during the shift, second, that the linear shape of reaction norms extend beyond the reference environment where they evolved and where variance is minimized, and third, that the new environment is far outside the distribution of past environments and causes a density-independent decline in population size.

Constant additive genetic variance, as modelled in our simulations and analytical results, is commonly assumed in models like ours to make analytical progress, and although it is not biologically realistic it can provide a good approximation to more complex dynamics ((38)). Accounting for evolving genetic variance would require added complexity, such as tracking the full distribution of breeding values: (e.g., (39)) or using an approximation like the stochastic house-of-cards (e.g., (23)). More explicit genetics have already been included in some models of ER (e.g., polygenic adaptation, ((40))), but this is more challenging with plasticity and a stochastic environment. With environmental noise in a constant environment, variance in reaction norms is expected to decrease ((41)) which could represent the state of the population in the long run after the phenotype has evolved to become canalized around the new mean environment ((42)) but theory is lacking for the transient change in genetic variance after the shift in mean environment. It is likely that reductions in population size during ER will reduce genetic variance. We investigated the sensitivity of our results to this possibility (see Figure S6) and still found that if predictability in the new environment is below a critical level then ER is unlikely. The critical levels in this case, however are higher than those suggested by our analytical results, which thus can be viewed as a lower bound on extinction risk.

The dynamics we show are all derived by assuming that partially adaptive plasticity with linear reaction norms extending to the novel environment and that variance is minimal in the reference environment. Theory demonstrates that linear reaction norms will evolve within the reference regime over long timescales (e.g., (18)), but it is the assumption that variance is minimal in the reference environment (and that reaction norms extend to the novel environment, ((9))) that implies increased heritable variation in the novel environment. Without increased additive genetic variance , there is no strong covariance between reaction norm slope and the trait expressed, which means no strong selection to increase reaction norm slope and so no transient increase in plasticity (see F). An environmental shift then increases neither variance load nor, necessarily, plasticity. Because we expect quite different results without assuming increases in the new environment, we emphasize that our models will only apply when novel or stressful environments reveal cryptic genetic variation (e.g., (7)). Although this idea finds support in some systems, in meta-analyses the opposite trend (of decreasing heritable variation in rare, stressful or novel environments) is equally frequent (reviews: (10); (11)). There are relatively few studies that obtain clear results either way, however, in part due to the difficulty of replicating an experiment across many environmental values ((12)). In these cases, then, evolution of plasticity should have little influence on ER and we expect dynamics to follow results for non-plastic ER [e.g., classic deterministic theory (3), or its extension to stochastic environments by Ashander and Chevin in prep]. Furthermore, although we assumed partially adaptive plasticity, it can be sometimes be maladaptive (e.g., (43)). Developing theory on this may require modelling non-linear reaction norms (e.g., via function-valued traits, ((44))).

Finally we ignore demographic regulation, effectively assuming that density-dependence is very weak. In practice we assume the novel environment is stressful and initially causes a density-independent population decline because the environment is far outside the previous range of environments. Introduced by Lande ((9)), it is an extreme version of environmental novelty, but one that yields mathematically tractable expressions for the growth rate. As shown previously without stochasticity, the density-independent trajectories we study are close to those under very weak density dependent regulation or a form of density-dependence, e.g., large exponent in a theta-logistic model, where growth trajectories during rescue are similar ((16)). Under stronger density regulation that acts even when the population is far below carrying capacity, we would expect steeper declines in population size.

Despite the limitations mentioned above, our assumptions apply nicely to some systems. For example, compare the implied increase in (-fold; Figure 1(a)) to Husby et al. ((45)) who showed higher temperature increased genetic variance of breeding time of the great tit Parus major (-fold increase in mean ), or to McGuigan et al. ((46)) who showed for three-spine stickleback (Gasterosteus aculeatus) an even stronger increase in genetic variance of body size with low-salinity (-fold increase in mean ). How frequently such increases in occur with environmental novelty remains an open question. Furthermore, even if increases occur, they are not sufficient to guarantee rescue. Among other conditions, heritability and evolvability must also increase, and as we show heritable variance in plasticity may be detrimental for several reasons (standing load, and stochastic lag load in unpredictable environments).

### 4.2 Empirical context and applications

We define a long-term persistence criterion by predicting the stochastic growth rate for hundreds to thousands of generations after the population’s mean trait has adapted to the new optimum. To apply this theory, either for empirical verification or for prediction, several types of data are needed. First, estimates for parameters governing the genetics of GxE interactions (i.e., additive genetic variances in several environments) and the trait’s effect on fitness can be obtained from common garden and other experiments ((14)). Second, information about the environmental sensitivity of selection can be measured using a single episode of selection (e.g., in Parus major: ((47); (48); (49))). Third, parameters for environmental predictability can be characterized statistically (e.g., (50)), but this requires knowledge, or assumptions, about the environment of selection.

These data requirements are challenging but achievable in several field and laboratory systems. Phenological traits, because they are developmental traits under strong selection and cues are often known, may be the best fit. Furthermore, these traits among the most observable biological responses to climate change ((51)). Germination timing of high altitude plants is particularly promising: winter temperature and snow melt cue development that is also genetically-influenced and under stabilizing selection (with risk of frost-killing if too early, or dessication if too late; ((52))), and optimal timing varies with altitude ((53)), so reciprocal transplants shift the optimum. Furthermore, optimal timing is shifting with climate change ((54)). Osmoregulatory traits are another candidate. Studies in the copepod Eurytemora affinis support a role for rapid evolution driven by plasticity in parallel adaptations to freshwater ((55)). For three-spine stickleback additive genetic variation in body size is higher in stressful low-salinity environments ((46)). ER has already been studied in several microorganisms (e.g., Pseudomona fluorescens, ((56))), some of which display phenotypic plasticity (e.g., Escherichia coli, Saccaromyces cerevisiae, reviewed in ((57))).

Quantitative predictions of persistence for populations currently undergoing ER could aid conservation planning, assessment of invasive species, and management of antibiotic or pesticide resistance. Although many examples of the latter are major-gene effects (for analysis of ER via a single gene, see e.g., ((37))), even these cases may include a quantitative contribution from minor genes ((40)). Our theory is relevant for applied contexts where plasticity is thought to aid persistence or invasion, including reintroduction for conservation purposes and invasive species control. For example invasive species are more plastic in response to added resources ((15)), suggesting more plastic species are better invaders. Our findings imply a more subtle prediction: a “filter” against invaders long-adapted to low-reliability cues (where the same cue is maintained from the native to invaded range), and a “shield” for regions where cues used by common invaders are unreliable. Applying this theory to a variety of systems might help resolve observed variation how plasticity changes following a disturbance (e.g., (58)).

### 4.3 Conclusion

Overall, evolving plasticity facilitates evolutionary rescue unless the new environment is unpredictable. If it is not, then large variance in plasticity might help lineages long-evolved to low-predictability environments adjust to novel environments with high predictability. These findings suggest the role of plasticity in longer-term evolution to changing environments is positive, but limited. The rapid increase in plasticity that occurs in our model is an example of the Baldwin effect ((9)) where plasticity increases in species colonizing stressful environments. (Baldwin actually proposed theory for the evolution of plasticity that is much more general than this ((59)).) This effect, and related processes, have often been mentioned as under-appreciated factors in evolution ((60); (61)). In novel environments that fluctuate with low predictability, however, we show that a transient increase in plasticity (Figure 1(b)) can impose a substantial load on average growth, and thus a barrier to ER. For populations whose plasticity evolved in response to low-predictability cues, then, the Baldwin effect (as embodied in the two-phase process studied here where plasticity increases) may have limited importance in adaptation. An implication of these results is that over long timescales where the environment has shifted frequently, we expect phenotypic plasticity (and its genetic variance) to be absent or strongly reduced, unless cues are consistently reliable.

Code Code used to perform the simulations is available ((33)).

Competing interests We have no competing interests.

Author contributions JA designed the study, carried out the analyses, and drafted the manuscript; LMC designed the study, guided the analyses, and helped write the manuscript; MLB designed the study, guided the analyses, and helped write the manuscript. All authors gave final approval for publication.

Acknowledgements Feedback from Swati Patel, Sebastian Schreiber, and Michael Turelli improved an earlier version of the manuscript.

Funding Support from IGERT (JA, NSF DGE-0801430 to P.I. Strauss), ContempEvol (LMC, ANR-11-PDOC-005-01), and FluctEvol (LMC, ERC-2015-STG-678140-FluctEvol).

## Appendix A Figures

Supplemental Information for

“Predicting evolutionary rescue via evolving plasticity in stochastic environments”

Jaime Ashander, Luis-Miguel Chevin, Marissa L. Baskett

## Appendix B Approximate dynamics of growth rate in terms of maladaptation x(t)

Our equation (2) shows the Malthusian growth rate is the sum of the maximum Malthusian growth rate , the variance load , and the lag load with , , and ,

 log¯W(t)=rmax+LG(t)+LL(t). (4)

Below, we derive an approximation to this where the only time dependence is through the squared maladaptation .

### b.1 Strength of selection

The strength of selection, , appears in both the variance load and lag load but depends inversely on the phenotypic variance, which changes due to stochastic variation in the environmental cue, which affects the phenotype. Therefore, we approximate the expectation of the strength of selection using a Taylor series: \linenomath

 E[S(εc(t))] =(E[σ2z(εc(t))]+ω2)−1−Var(σ2z)2E[σ2z]3+O(Var(σ2z)3E[σ2z]4) ≈(E[σ2z(εc(t))]+ω2)−1.
\endlinenomath

Expanding using equation (1) of the main text and taking expectations, we get the approximate expected strength of selection , assuming , which we denote

 Sδ=E[S(εc(t))]≈(σ2a+σ2b(δ2+σ2c)+σ2e+ω2)−1. (5)

Where we use the and the mean environmental shift .

The first term of equation (4) is the phenotypic variance load at time

 LG(t)=1/2logS(εc(t))+1/2logω2.

Taking expectations, and using Taylor series: \linenomath

 E[LG] =1/2logω2+1/2E[logS(εc(t))] ≈1/2logω2+1/2log(E[S(εc(t))])−Var(σ2z)2E[σ2z]2 ≈1/2logω2+1/2log(E[S(εc(t))]) ≈1/2logω2+1/2logSδ
\endlinenomath

where the right hand side of the last line is the approximate expected variance load after using from (5).

The second term of equation (4) is the lag load at time , . Here we show that in expectation, this load has two components, a stochastic load and a shift load. To derive expressions for these, we take the expectation of the entire second term: \linenomath

 E[LL(t)]= E[S(εc(t))2x(t)2], = E[S(εc(t))2]E[x(t)2]+Cov[S(εc(t))2,x(t)2],
\endlinenomath

where we used the identity . The covariance in the final line represents how stochastic changes in environment (and optimum) that cause maladaptation also cause either larger or smaller phenotypic variance, depending on the direction, due to our assumption that genetic variance in plasticity increases with . Furthermore, under weak stabilizing selection, variance in contributes little to the strength of selection . For both these reasons, we expect this covariance to be small.

If we neglect the covariance term and use the approximate variance load from B.2 ( in main text) the expectation is . Removing the expectation, we use this as to approximate the dynamics of the lag load

 LL(t)=Sδ2x(t)2.

### b.4 Full dynamics

Bringing together the approximations developed above, we have an expression for the growth rate, where all components except the maladaptation are averaged over the fluctuations,

 r(t)≈rmax+1/2log(ω2Sδ)−Sδ2x(t)2. (6)

From this, we can obtain both the average long-run growth rate, by taking an expectation to obtain (because ). This shows the lag load consists of two loads. The first, expected load from reduction in log mean fitness due to maladaptation of the mean trait relative to the mean optimum in the new environment, is “shift load” . The second, expected load due to reduction in log mean fitness due to random fluctuations in the environment, is “stochastic load” . We can also analyse the variance in trajectories of , and thus population growth because .

Both of our subsequent analyses require computing the dynamics of mean squared maladaptation and the variance in maladaptation .

## Appendix C Dynamics of shift load depend on mean maladaptation ¯x(t)2

In this section, we derive the dynamics of the shift load under the approximation introduced by Lande ((9)) that separates adaptation into a fast Phase 1 and a slow Phase 2. We demonstrate that the population is approximately perfectly adapted in mean trait value by the end of Phase 1. We first write down the mean trait dynamics without the approximation, then describe the timescales of the two phase approximation, and derive approximate dynamics of the shift load owing to maladaptation in the mean trait.

### c.1 Mean trait

Trait dynamics follow the standard equation , where , is the additive genetic variance-covariance matrix, and is the selection gradient. The selection gradient on reaction norm height and slope obtained by taking the log-gradient of (from the equation for fitness above eq. (2) of main text; (9)) is

 β=−S(εc(t))(¯a(t)−A+¯b(t)εc−Bϵs(¯a(t)−A+¯b(t)εc−Bϵs)εc). (7)

With a constant additive genetic variance-covariance matrix ( matrix), the change per generation in and is given by

 Δ(¯a(t)¯b(t))=(σ2a00σ2b)β.

In the new environment, the expectation of the change per generation conditional on and is \linenomath

 Δ(E(¯a(t))E(¯b(t))) =E[Gβ] ≈−SδG(¯a−A+¯bδ−BδE[(¯a−A)εc(t)]+E[¯bε2c(t)]−BE[εcεs]),
\endlinenomath

where the approximation comes from the treating as a constant. We assume no environmental tracking by phenotypic plasticity, i.e., so , or by reaction norm elevation, i.e., so . Note that the tracking of the environment by the reaction norm elevation could be included, reducing the expected mean plasticity ((28)), but we neglect it here. However, we do allow for environmental tracking when computing the lag load (below). Then, using and , the expectation of the change is approximately

 E[Δ(¯a(t)¯b(t))]≈−SδG[(1δδδ2)(¯a(t)−A¯b(t)−B)+(0(¯b(t)−ρB)σ2)]. (8)

The approximation is exact if . Note that this differs from Lande ((9)), where this relation was treated as exact ((28)).

We solve for the explicit trait dynamics relative to the long-run equilibrium state ((9)). Setting selection gradient in equation (7) to zero, solve for long run trait values

 (¯a∞¯b∞)=(A+Bδ(1−ρδ)Bρδ) (9)

One can show, with some algebra ((9)), the one-generation change (8) is the product -, where is the difference between mean trait values and their long-run equilibrium values computed above and

 ~G=(σa2σa2δσb2δσb2δ2(1+σ2δ2)).

The expected dynamics can then be expressed in terms of the eigenvectors and eigenvalues of the matrix ((9); (16)),

 z(t)=c1e1(1−Sδλ1)t+c2e2(1−Sδλ2)t, (10)

where and are eigenvectors and eigenvalues of respectively and terms are constants determined by initial conditions.

### c.2 Approximation for a large environmental shift

As in earlier work, we consider the case where the shift in the mean environment is very large relative to background noise, . Here, we initially write down the eigenvalues and eigenvectors of to first order in this small term, but thereafter follow Chevin and Lande ((16)) in deriving approximate dynamics to leading order. To first order in , the eigenvalues are

 (λ1λ2)=⎛⎜⎝(σ2a+δ2σ2b)+δ2σ2bσ2δ2ϕσ2aσ2δ2ϕ⎞⎟⎠,

and the eigenvectors are

 e1=(δ(1−ϕ)(1−σ2δ2ϕ)ϕ),e2=(δ(1+σ2δ2ϕ)−1).

These match the calculations of ((9)) up to a constant of (equivalent to in Lande ((9))).

Assuming the population has long evolved in an environment with predictability , the initial trait values are . Using (9), the initial conditions in the re-centred trait are

 z0=(−Bδ(1−ρδ)B(ρ−ρδ))

To leading order, the constants are

 (c1c2)=(−B(1−ρ)−B(ρ(1−ϕ)−ρδ+ϕ))

#### Timescales of phases 1 and 2

If most phenotypic variation in the new environment is due to variance in plasticity, , and the shift in the mean environment is large (relative to background variability as in our approximation above), the trait change takes place in two phases that occur at very different timescales ((9)). When selection is weak, geometric terms in (10) can be replaced by exponential terms , indicating the relative timescales of change along the eigenvectors are given by . The ratio , and when much of the additive genetic variation is due to variation in plasticity so , then , which is small in the approximate case we treat. Then, change along occurs very fast relative to change along ((9)). At the end of Phase 1, while Thus, the approximate state of the system relative to its final state, i.e., (9), is . The trait values, at the end of Phase 1 are, to leading order,

 E[(¯aO(t1)¯bO(t1))]≈(A+Bδ(1−ρ)(1−ϕ)B(ρ+ϕ(1−ρ))) (11)

The effect of the initial environment occurs through predictability , which under our assumption that the population is adapted initially also determines the initial mean plasticity . We see the initial plasticity has a strong influence at the end of Phase 1 only if is small. When is large, the plasticity at the end of Phase 1 is close to “perfect” i.e. . Note also that to first order, the mean phenotype is perfectly adapted .

Where the extinction risk is calculated at half of the characteristic timescale of Phase 2, i.e.,

#### Trait dynamics during Phase 1

Throughout Phase 1, the term , so the dynamics are given by

 z(t)=c2e2+c1e1(1−Sδλ1)t.

We again replace the geometric term with an exponential (valid for weak selection) and re-normalize. To leading order, after some rearranging, the right hand side of the expected dynamics is

 E[(¯a(t)¯b(t))]=(1−e−tSδλ1)B(1−ρ)(δ(1−ϕ)ϕ)+(ABρ), (12)

which is analogous to the result of Chevin and Lande (Supporting Information, eq A6; (16)). As in that paper, we compute the eigenvalue only to leading order in so which is equivalent to the expression used in Chevin and Lande ((16)).

#### Trait dynamics during Phase 2

During Phase 2, the term , so the dynamics are given by

 z(t)=c2e2(1−Sδλ2)t.

Then, using (9) and again replacing the geometric term with an exponential (valid for weak selection) and re-normalizing, the expected dynamics to leading order in during Phase 2 are

 E[(¯a(t)¯b(t))]=(A+Bδ(1−ρδ)Bρδ)−B(ρ−ρδ+ϕ(1−ρ))(δ−1)e−tSδσ2aσ2δ2ϕ, (13)

where we have used . Note that for , this exponential term equals 1 and this equation agrees with (11).

### c.3 Dynamics of expected maladaptation during Phase 1

We derive the expected maladaptation of the mean trait during an initial phase of evolutionary rescue, focusing on a case where the size of the environmental shift is large and much of the additive genetic variance in the new environment owes to genetic variance in reaction norm slope. The shift load (computed in B) is . After we compute the dynamics of the mean maladaptation , we will have an approximation for dynamics of the shift load, \linenomath

 ¯x(t)2=E[x(t)]2= (E[¯a(t)]−A+E[¯b(t)εc(t)]−BE[εs])2 ≈ (E[¯a(t)]−A+δ(E[¯b(t)]−B))2.
\endlinenomath

The last equation comes from assuming the covariance between reaction norm slope and the cuing environment is small relative to the mean value of the new environment. Then, in the new environment. This is reasonable when . After using (12) and some algebra, we obtain

 ¯x(t)2≈B2δ2(1−ρ)2e−2tSδσ2a1−ϕ. (14)

Where we use to represent the proportion of additive genetic variation in the new environment due to variation in plasticity,

 ϕ=δ2σ2bσ2a+δ2σ2b.

Equation (14) indicates the shift load goes to zero as increases.

### d.1 Perceived environment with fixed plasticity

A tactic from Michel et al. ((35)) aids in calculating the variance as a function of fixed plasticity. We define the perceived optimum as the difference between the optimum and the mean trait after accounting for the plastic response so that . Then, the perceived variance in the optimum is

 σ2ψ(¯b∗,ρδ)=σ2(B2+¯b∗(¯b∗−2Bρδ)), (15)

and autocorrelation in the perceived optimum is

 Tψ(¯b∗,ρδ)=−log[ρ1/τδ(1−B¯b∗(ρ−1δ−ρδ)B2+¯b∗(¯b∗−2Bρδ))]−1 (16)

((35)). We can then express maladaptation in terms of the intercept and perceived environment as .

### d.2 Variance at stationarity

We derive an approximation for variance in maladaptation under stationarity, which we denote . In practice, this means we solve for the effect of fluctuations on maladaptation after a long time, we assume the mean maladaptation is zero, and we also assume fixed mean plasticity . We are interested in finding an asymptotic expression for the variance of this term.

Assuming fixed plasticity, all change in the trait occurs through change in the reaction norm height,

 Δ¯z=Δ¯a(t)=−Sδσ2ax(t).

When selection is weak relative to genetic variance in reaction norm height, and the fluctuations in the perceived environment are not large, evolution can be approximated in continuous time ((22); (35)) as

 dxdt+Sδσ2ax=−dψdt,

where . For and constant genetic variance , the solution to this differential equation is

 ¯a(t)=Sδσ2a∫∞0exp(−Sδσ2aτ)ψ(t−τ)dτ. (17)

What remains is to compute . Using our quasi-stationarity assumption, we need only compute . The last of these expectations is simply the variance of the perceived environment . The first and second expectations integrate over time (from eq. 17), \linenomath

 −2E[¯a(t)ψ(t)] =−2Sδσ2a∫∞0exp(−Sδσ2aτ)E[ψ(t)ψ(t−τ)]dτ E[¯a(t)2] =S2δσ4a∫∞0∫∞0exp(−Sδσ2a(τ1+τ2))E[ψ(t−τ1)ψ(t−τ2)]dτ1dτ2.
\endlinenomath

Because is a linear combination of autoregressive Gaussian processes and , we can express the expectations involving in terms of autocovariance and . In both of these expressions, is the characteristic autocorrelation time of the perceived environment .

The first expectation is \linenomath

 −2E[¯a(t)ψ(t)] =−2Sδσ2aσ2ψ∫∞0exp(−τ(Sδσ2a+1/Tψ))dτ, =−2TψSδσ2aσ2ψ(TψSδσ2a+1).
\endlinenomath

The second expectation involves an absolute value term, meaning the integral must be taken in two parts, and evaluates to \linenomath

 E[¯a(t)2]= S2δσ4a∫∞0∫∞0exp(−Sδσ2a(τ1+τ2))σ2ψexp(−|τ1−τ2|/Tψ)dτ1dτ2, = Sδσ2aσ2ψTψ(TψSδσ2a+1)
\endlinenomath

Combining these expressions, \linenomath

 E[(¯a(t)−ψ(t))2] =σ2ψ(−2TψSδσ2a(TψSδσ