Fluctuations in Lifetime Selection in an Autocorrelated Environment

Most natural environments vary stochastically and are temporally autocorrelated. Previous theory investigating the effects of environmental autocorrelation on evolution mostly assumed that total fitness resulted from a single selection episode. Yet organisms are likely to experience selection repeatedly along their life, in response to possibly different environmental states. We model the evolution of a quantitative trait in organisms with non-overlapping generations undergoing several episodes of selection in a randomly fluctuating and autocorrelated environment. We show that the evolutionary dynamics depends not directly on fluctuations of the environment, but instead on those of an effective phenotypic optimum that integrates the effects of all selection episodes within each generation. The variance and autocorrelation of the integrated optimum shape the variance and predictability of selection, with substantial qualitative and quantitative deviations from previous predictions considering a single selection episode per generation. We also investigate the consequence of multiple selection episodes per generation on population load. In particular, we identify a new load resulting from within-generation fluctuating selection, generating the death of individuals without significance for the evolutionary dynamics. Our study emphasizes how taking into account fluctuating selection within lifetime unravels new properties of evolutionary dynamics, with crucial implications notably with respect to responses to global changes.


Introduction
The study of adaptation to changing environments has received renewed attention in the context of current global changes induced by anthropogenic activities. Many theoretical studies on adaptation to changing environments have focused on cases where environmental change follows a deterministic trend (Lynch et al. 1991, Lynch and Lande 1993, Gomulkiewicz and Houle 2009, Cotto and Ronce 2014, motivated by major tendencies such as climate warming. Yet, most variation on short timescales results from stochastic fluctuations, and climatic time series generally include stochastic variations around their main trends (Stocker et al. 2013). The variance and autocorrelation of these fluctuations need to be modeled to accurately predict the impact of climate warming on biosystems (e.g. Katz 1996, Rowell 2005. Global climatic dynamics have recently been suggested to affect This work is licensed under a CC BY 4.0 International license. local selection pressures (Siepielski et al. 2017), and time series of climatic variables are generally autocorrelated in time, notably as a result of thermic inertia (Hasselmann 1976, Rowell 2005. As most climatic (and other ecologically relevant) environmental variables are autocorrelated in time (Vasseur and Yodzis 2004), environmental selective pressures are likely to also be temporally autocorrelated, an aspect that is often neglected in models of adaptation that do include stochasticity in the environment (e.g. Lynch and Lande 1993, Engen et al. 2011.
Theoretical models of adaptation to stochastically changing environments often assume that the fitness of an individual depends on the match between its phenotype and an optimum phenotype influenced by the environment (so-called "moving optimum" models, reviewed by Kopp and Matuszewski 2014). When the environment undergoes stationary stochastic fluctuations, the effect of evolutionary responses on long-term fitness and expected population growth depends on the variance and autocorrelation of changes in the optimum phenotype (Charlesworth 1993, Lande and Shannon 1996, Chevin et al. 2017). If the stochastic fluctuations of the optimum are not autocorrelated, then responses to selection in a given generation might increase maladaptation in the next, with detrimental consequences for population growth (Lande and Shannon 1996). In contrast, temporal autocorrelation of the optimum makes the environment more predictable, such that the evolutionary response in a given generation is more likely to be beneficial in the next, and that any factor improving the response to selection leads to a higher expected long-term growth rate of the population Shannon 1996, Chevin 2013). The autocorrelation of the optimum also affects the variance and overall shape of the probability distribution of population size of an evolving population, with consequences for extinction risk (Chevin et al. 2017). Lastly, temporal autocorrelation of the environment is a major driver of the evolution of plasticity, with higher reaction norm slopes evolving under more predictable (autocorrelated) environments of selection (Gavrilets and Scheiner 1993, Lande 2009, Chevin et al. 2015. On the empirical side, a number of studies undertook to determine the temporal characteristics of selection in natural populations. Temporal variation in selection is observed in classic long-term surveys, like beak shape in Darwin's finches (Grant and Grant 2002), banding patterns in Cepea snails (Cain et al. 1990) or spine number in threespined sticklebacks (Reimchen andNosil 2002, reviewed in Bell 2010). Consistent with these observations, studies analyzing several datasets with long-term selection estimates using a common framework also found that selection varies in strength and direction (Kingsolver andDiamond 2011, Siepielski et al. 2011, but see Morrissey and Hadfield 2012). Further, the development of a statistical framework that models temporal fluctuations in phenotypic selection as a random process, characterized by the variance and autocorrelation of an optimum phenotype, allowed to show that the optimum laying date in a population of Great tits undergoes autocorrelated temporal fluctuations (Chevin et al. 2015, Gamelon et al. 2018). This framework rests on a model of a Gaussian fitness peak, as in many theoretical studies (Kopp and Matuszewski 2014), so that the estimated parameters should have a direct theoretical interpretation.
definitions of fitness and selection may not be easily transposed from one to the other. Even regardless of theoretical considerations about the definition of long-term expected fitness in a stochastic environment, understanding how variation in selection integrates within lifetime is not straightforward. Most theoretical studies that included autocorrelation in the environment neglected the dynamics of selection within generations, by assuming (often implicitly) that selection occurs instantaneously, in discrete non-overlapping generations. In other word, lifetime fitness is assumed to result from a single instantaneous episode of selection; or at least, the unfolding of selection within lifetime is not modeled explicitly, making it ambiguous whether the per-generation optimum represents instantaneous or repeated selection. More realistically, selection occurs in several episodes along a generation (Lande 1982, Arnold andWade 1984b). In fact, the vast majority of long-term estimates of selection gradients focused on a single fitness component (~ 98% among 2819 selection estimates in Kingsolver and Diamond 2011), and therefore correspond to a single episode of selection, among several possible others in a generation or lifetime. This has several related consequences for the interpretation of measurements of selection. First, the strength of selection over a single episode cannot necessarily be compared to predictions based on overall selection per generation, because this selection episode needs to be weighted by its contribution to the lifetime (or per-generation / per-unit-time) fitness. Second, both the magnitude and predictability of temporal variation in selection is probably not well captured by treating each selection episode in isolation. The reason is that in organisms that undergo multiple selection episodes, environmental variability gets integrated over lifetime. This has important consequences, notably regarding the contribution of extreme events to selection, a topic of recent interest (Grant et al. 2017, Marrot et al. 2017).
The study of demography has long questioned the relationship between variation in vital rates that contribute to elements of the transition matrix (e.g. fecundities, survival rates) and various measures of the population growth rate (Tuljapurkar et al. 2003, Morris et al. 2004, Saether et al. 2013. This theory has been used to analyse fluctuating selection in stochastic environments that are not autocorrelated (Engen et al. , 2014. However, how stochastic autocorrelated fluctuating selection mediated by different fitness components integrates over lifetime is an under-investigated topic, limiting our understanding of adaptation to randomly changing environments. Making progress on this question requires extending previous approaches that integrate the multiplicity of selection episodes along life. Of particular note is the approach by Wade (1984b, 1984a), who proposed and applied a decomposition of lifetime selection into a series of selection episodes with multiplicative effects on fitness, resulting in selection gradients (dependent on log-fitness) that sum along a generation (see also McGlothlin 2010). Integrating the effect of successive selection episodes on lifetime fitness is appropriate to organisms with non-overlapping generations (Arnold and Wade 1984a), but has not been used to investigate the effect of environmental stochasticity.
We here build a model of a phenotypic trait exposed to multiple selection episodes, in nonoverlapping generations. We assume that the phenotype of an individual at the focal trait does not change between selection episodes. This could correspond for example to a morphological trait that only influences selection after it is fully developed, or a coloration trait that affects survival probability (due to e.g. predation, or thermoregulation) prior to reproduction. Under these assumptions, we detail and clarify how the temporal structure of selection during lifetime influences the predictability of lifetime viability selection, and the load caused by maladaptation in a stochastic environment. Our model applies directly to organisms that are semelparous with non-overlapping generations, including univoltine insects such as cicadas, and species with economic or cultural significance such as salmons or bamboos. Our analysis further applies to cohort analysis (see discussion) when considering more general life histories with overlapping generations. As such, our study provides a necessary step toward a complete understanding of how life history affects evolutionary dynamics in the context of autocorrelated fluctuating selection. Overall, we find that the way selection operates during a lifetime strongly determines the impact of environmental fluctuations on variation in selection across generations, and the incurred maladaptation load in a stochastic environment.

Selection model
We investigate the evolution of a single quantitative trait z in an organism with nonoverlapping generations. Along their life cycle, individuals undergo several selection episodes that determine their probability to survive between stages before reproduction, with the latter occurring only at the end of each generation. We further assume that there is no selection on fecundity, which we denote B. For simplicity, we focused on a single trait, but our framework can be readily extended to the multivariate case, with different (correlated) traits undergoing selection at different stages (Arnold and Wade 1984b, Arnold and Wade 1984a, Cotto and Ronce 2014. We further assume that the phenotype z at the beginning of a generation can be decomposed as z = z + e, where z is the additive genetic value and e is a residual component of variation independent of the additive genetic value, and with normal distribution of mean 0 and variance E (Falconer and Mackay 1996). We assume that many loci contribute to the additive genetic value, such that z is normally distributed with mean z¯ and variance G. Under these assumptions, z is also normally distributed with mean z¯ and variance P = G + E.
The survival probability of an individual with phenotype z at age i in generation g depends on the match between its phenotype z and an optimum phenotype ν g,i that changes with the environment, where W max,i is the maximum survival rate at age i and ω i measures the width of the Gaussian survival function (inversely related to the strength of selection). We assume that selection occurs through all transitions preceding reproduction, so that the components of fitness are multiplicative (Arnold and Wade 1984a). If reproduction could occur at several stages, then the components of fitness would no longer be multiplicative, and another framework (e.g. overlapping generations, Charlesworth 1994) should be used. The total absolute fitness of an individual with phenotype z is then (where the subscript g has been dropped for simplicity) Under the assumption that selection is Gaussian at each episode (eq. 1), the total fitness function is also Gaussian, where, ω tot 2 = ∑ i 1/ω i 2 -1 is the total strength of selection (that is the harmonic mean of the width of the selection function at each selection episode), θ tot = ∑ i e i θ i is the per-generation effective optimum corresponding to the sum of the phenotypic optima at each selection episode, weighted by the contribution of these episodes to the total strength of selection. We denote this contribution as e i = ω tot 2 ω i 2 , by reference to elasticities in stage-structured population models (Caswell 2001). The maximum absolute fitness over a generation is With the further assumption that the strength of selection (here measured as a width of the fitness function) is the same at all selection episodes, we have ω i = ω and e i =1/n, so the equations above further simplify as ω tot 2 = ω 2 /n (5) where E^[θ] and V^[θ] are the "sample" expectation and variance (respectively) of θ within a generation.
The change in the mean phenotype over an entire generation, after reproduction, is (Lande 1976) where z¯n is the mean phenotype after the last episode of selection, and additive genetic variance G 0 and heritability ℎ 0 2 = G 0 /P 0 are measured at the beginning of a generation, prior to any selection. Indeed, each selection event can generate a covariance between P and E that needs to be accounted for if dealing with selection episode by episode (as in Arnold and Wade 1984a), resulting in complications (see Supp. Mat.) that can be overcome by focusing on the effect on total selection on the phenotype distribution prior to any selection. With Gaussian stabilizing selection (eq. 1) and a normally distributed trait, the mean fitness in the population is where S tot = 1 P 0 + ω tot 2 . From equation 8, it follows that the per-generation change in the mean trait value in the population is where β tot is the total selection gradient, and which has a similar form as in models with a single selection episode (e.g. Lande 1976).

Fluctuations of the environment and timing of selection
We are interested in cases where the optimum phenotype fluctuates because of random fluctuations in the environment, i.e. fluctuations of the optimum phenotype represent those of the environment. For simplicity, we assume that the optimum phenotypes at all selection events in a life cycle respond to the same environmental variable, so that fluctuations in the optimum during a lifetime can be modeled as a temporal sampling of the same stochastic process. We consider the case where the optimum follows a first-order autoregressive process with mean 0 (without loss of generality), defined by where ρ is the autocorrelation of the optimum over one (absolute; e.g. one year) time unit (which we assume is positive), σ θ 2 its stationary stochastic variance, and ζ t a standard normal random deviate. For such a process, the autocorrelation of the optimum over k units of time is ρ k =exp(-k/T θ ), where T θ =-1/ln(ρ) is the characteristic timescale of the autocorrelation of the optimum. Note that the autoregressive process converges to white noise when ρ tends to 0, i.e. when the optimum phenotype is correlated over an infinitesimally small timescale. Cotto and Chevin Page 6 Theor Popul Biol. Author manuscript; available in PMC 2020 October 01.
We assume that n episodes of selection occur before reproduction. These selection episodes all take place over a time-window that spans a fraction α of generation time T. The timewindow for selection remains unchanged across generations. This models the fact that all traits in an organism need not be exposed to selection throughout life. Instead, it is more likely that there is a time-window within a generation during which a given trait may be sensitive to the environment and exposed to selection (e.g. during the first stages). When exactly the time-window of selection occurs within a generation (e.g. at the beginning or at the end) has no effect for the remaining of the analysis. Under our assumption of nonoverlapping generations, only the duration of this window matters. For simplicity, we assume that the episodes of selection are evenly spaced along the time-window for selection, so that they occur every αT/n units of time. Using the index g for the number of generations, we denote as z¯g ,0 the mean trait value before the first episode of selection of generation g.
The timescale of life history relative to that of the environment is represented in figure 1.

Predictability of selection
The predictability of selection can be defined in several ways. For our purpose, we will focus on the autocorrelation of the effective optimum per generation, and that of the total selection gradient, because these relate to measurable quantities (Lande and Arnold 1983, Chevin et al. 2015, Gamelon et al. 2019) known to affect adaptation and population dynamics in a stochastic environment (Lande and Shannon 1996, Chevin 2013, Chevin et al. 2017. We first calculate the autocorrelation of θ tot , before investigating the properties of the fluctuations in the selection gradient β tot , which despite being less directly connected to theoretical predictions, are closer to what is commonly estimated empirically (Lande and Arnold 1983).
Autocorrelation of the effective optimum-The autocorrelation of the effective optimum at a lag of one generation is (see also Supp. Mat.) ρ θtot = ρ T ∑ k = 1 n e k 2 + ∑ i = 1 n − 1 ∑ j = i + 1 n e i e j ρ j − i αT /n 1 + ρ 2 i − j αT /n ∑ k = 1 n e k 2 + 2∑ i = 1 n − 1 ∑ j = i + 1 n e i e j ρ j − i αT /n = ρ T φ, where ρ T is the autocorrelation of the optimum over a generation time, and its autocorrelation between two consecutive selection episodes is ρ αT/n (which differs from the autocorrelation ρ per time-unit). In the following we detail the inflating factor φ corresponding to the complex fraction that appears in the first equality in equation (12). The first sum in the numerator results from the covariance between the optimum in a given episode of selection and the same episode in the next generation. The second sum results from the covariances of the optimum between two different selection episodes in two consecutive generations. The denominator has a similar structure, and results from similar effects (variance per episode vs covariance across episodes) within a generation, contributing to the variance of the effective optimum (see Supp. Mat.). It can be shown (by noting that 0<ρ (j-i)αT/n )≤1, so that ρ -(j-i)αT/n )>1) that the ratio in equation (12) is always larger than 1, and is thus an inflating factor, which always increases the autocorrelation of the effective optimum relative to the autocorrelation of the environment, as a results of the averaging of selection between episodes within generations. Remarkably, this result extends to the autocorrelation over τ>1 generations when the autocorrelation of the environment decreases exponentially with time, as occurs under a Markovian process such as our firstorder autoregressive process. We thus have for the autocorrelation of the effective optimum overτ≥1 generations equation (13) shows that, once the inflation factor has been accounted for, the autocorrelation of θ tot declines exponentially at the same rate as the autocorrelation of θ ( fig.   2A). In principle, then, information on the autocorrelation of the effective optimum between two generations allows to determine its autocorrelation over many generations, providing that the autocorrelation of the environment (here also corresponding to the autocorrelation of θ) driving selection is known.
More insights arise from the simpler case where the strength of selection is the same at all selection episodes (eq. 5-6), wherein φ simplifies to where ψ = ρ αT is the autocorrelation of θ over the time window for selection. equation (14) shows that when variance in the strength of stabilizing selection within generation can be neglected, then the inflation factor does not depend on the strength of selection, but only on environmental autocorrelation ψ over the time window for selection, and on the number of selection episodes. The inflation factor φ increases when the autocorrelation of the  2B). In other words, weakly autocorrelated environments can actually correspond to much larger autocorrelations of variables relevant to selection. For example, when there are many episodes of selection per generation and the autocorrelation of the environment is 0.2, the autocorrelation of θ tot can be increased by more than 50% as compared to the autocorrelation of θ ( fig. 2B). When the number of episodes of selection is large, φ converges to which only depends on environmental autocorrelation ψ over the time window for selection.
We found that equation (15) provides a good approximation for equation (13) as soon as there are more than a few episodes of selection per generation, and even if the actual strength of selection varies across episodes ( fig. 2B). Interestingly, equation (15) does not depend on the strength of selection at each episode, and proposes a simple transformation between the autocorrelation of the environment and the autocorrelation of the effective optimum over a generation time. Substituting equation (15) into equation (13) thus provides, in principle, a simple way to estimate the autocorrelation of the effective optimum in natural populations from a measure of the autocorrelation of the environment, over the duration of selection within a generation on one hand (to compute the inflation factor), and across generations on the other hand. The longer the organism is sensitive to selection within a generation, the larger the inflation factor ( fig. 2C). Lastly, when the environment is not autocorrelated, the effective optimum is not autocorrelated either, showing that withingeneration selection alone does not generate autocorrelation of the effective optimum.

Fluctuations in directional selection.
Variance of the selection gradient: Ultimately, evolutionary change in the mean phenotype in response to selection is directly related to the total per-generation selection gradient β tot (eq. 10). At stationarity after many generations, we find that the variance of the selection gradient is approximately (see Supp. Mat.) Var β tot ≈ 2 Var θ tot S tot The variance of the selection gradient increases with the variance of fluctuations in the effective optimum (and of the environment, see Supp. Mat.) and with the total strength of selection (eq. 16 and fig. 3). equation (16) further shows that the variance of the selection gradient decreases when the inflation factor increases. Correspondingly, Var[β tot ] decreases with α (because the inflation factor increases with α, fig. 2C) and with the number of selection episodes ( fig. 2B and 3A). In other words, recurrent selection along a generation tends to decrease the variance of the selection gradient, as expected since environmental variation is averaged among selection episodes.
Less expected is the effect of environmental autocorrelation. We find that with several selection episodes per generation, the variance of the selection gradient depends nonmonotonically on the autocorrelation of the environment, with a maximum when the environment is moderately autocorrelated ( fig. 3B). This result contrasts with the case with a single selection episode per generation, where the variance in the selection gradient decreases with increasing autocorrelation of the environment, as a result of better adaptive tracking of the optimum Shannon 1996, Chevin andHaller 2014). This difference arises because, with several episodes of selection, positive autocorrelation in the environment not only produces autocorrelation of the total optimum θ tot across generations, but also increases the variance of θ tot (see eq. A1 and fig. A1 in Supp. Mat.), with antagonistic effects on adaptive tracking by genetic evolution, and hence on the variance of the total selection gradient across generations Shannon 1996, Chevin andHaller 2014).
We also investigated the influence of variable selection strength, by comparing the variance of the selection gradient obtained for a fixed strength of selection per episode (eq. 5-6) with the case where the strength of selection is drawn randomly at each episode from a distribution with the same expected total selection S tot (details in fig. 2 caption). We found that when there are only few episodes of selection per generation, variance in total selection per generation increases the variance of the selection gradient relative to what is expected in the deterministic case. However, with many selection episodes per generation, the variance in total selection decreases, resulting in a good fit between the simulations and the deterministic expectation ( fig. 3, compare crosses with dashed lines and circles). The effect of random selection per generation increases when selection is strong ( fig. 3). Lastly, we compared predictions from equation (16) calculated from a discrete generation model to previous results using a continuous time approximation (Lande and Shannon 1996 eq. 7 and table 1, and Chevin and Haller 2014 eq. 5a). The continuous time approximation provides more compact equations but tends to underestimate the variance of the selection gradient when selection increases ( fig. 3, compare dashed and continuous lines).

Autocorrelation of the selection gradient:
For the autocorrelation function of the total selection gradient across generations, we start by analyzing its expression in the case of a single selection episode per generation (i.e. when φ = 1). We obtain (see Supp. Mat.) for selection gradients τ generations apart where κ ≈ S tot G 0 1 + ρ T 2 -ρ T (see eq. A15). equation (17) is similar in form to an expression obtained previously using a continuous-time approximation (eq. 5a in Chevin and Haller 2014), and leads to similar conclusions. The autocorrelation of the selection gradient with one selection episode is a linear combination of two exponentially decreasing functions, corresponding on the one hand to the autocorrelation of the environment (first term in numerator), and to evolutionary responses to deviations from the optimum (second term in the numerator) on the other hand. Consistent with the result of Chevin and Haller (2014), we find that for large autocorrelation in the environment, equation (17) converges to (1-S tot G 0 ) τ , demonstrating that inertia in responses to selection sets an upper limit to the autocorrelation of the selection gradient.
In the general case where selection can occur in several episodes within a generation, we find a more complex expression for the autocorrelation of the selection gradient, which can nevertheless still be written in the form where the detailed expressions of ε, γ and ϑ depend on the per-generation evolutionary potential S tot G 0 , the autocorrelation of the environment over a generation ρ T , and γ and ϑ depend additionally on the inflation factor φ. The details of these terms can be found in the calculations provided in supplementary material (eq. A14). Equations (17) and (18)  Comparing equations (17) and (18) highlights that with multiple selection episodes, the autocorrelation of the environment ρ and the rate of response to selection S tot G 0 do not have in general an additive influence on the autocorrelation of selection gradients across generations, because the factor φ multiplying the former depends on the latter, and reciprocally for γ. This occurs because with multiple episodes of selection, the relevant autocorrelation is that of the effective optimum, which integrates the effect of selection across and within generations. The autocorrelation of the selection gradient increases most with the number of selection episodes when the autocorrelation of the environment is low ( fig. 4A), consistent with what we found for the inflation factor ( fig. 2B and eq. 18). We also find that for a given number of selection episodes per generation, the autocorrelation of the selection gradient slightly decreases when the total strength of selection increases ( fig. 4B).
Over time intervals longer than a generation, the autocorrelation function of selection gradients for a single selection episode (eq. 17) and for multiple episodes (eq. 18) have a similar shape ( fig. 5), because they depend on the autocorrelation of the environment and of the effective optimum (respectively), which have the same exponential rate of decrease ( fig.  2A). Importantly however, the autocorrelation of βtot in the first few generations -where the effect of multiple episodes is most pronounced -weights most in determining the expected maladaptation in a given generation (see e.g. eq. A6 in Supp. Mat.), so that differences in autocorrelation over short intervals of generations can have profound effects. It is worth mentioning that equations (17) and (18) can produce negative autocorrelations ( fig. 5) when the effect of the evolutionary inertia (second term at the numerator) grows larger that the effect of the autocorrelation of environment (see also Chevin and Haller 2014).
Lastly, and similarly to the variance of the selection gradient, we found that when the strength of selection is drawn randomly at each selection episode (see previous related paragraph for variance), the autocorrelation of selection is below its expectation under constant selection strength, when there are only few selection episodes per generation (fig .  4A). This difference is especially marked when the environment is very autocorrelated and when selection is strong ( fig. 4A, B) and vanishes as the number of selection episodes per generation increases.

Maladaptation load
Much of the work on adaptation of quantitative traits to autocorrelated stochastic environments has focused on the resulting expected lag load caused by deviations of the mean phenotype from the optimum, as this directly connects to long-term population growth and extinction risk Shannon 1996, Chevin et al. 2017). From equation (10), the log mean fitness in the population (which relates to the relative increase in population size) is

Europe PMC Funders Author Manuscripts
Europe PMC Funders Author Manuscripts where r max = lnB + ∑ k = 1 n r max, k corresponds to the maximum potential fitness, for a hypothetical population that would be monomorphic for the optimal phenotype in each selection episode. The last two terms describe reductions in fitness caused by maladaptation of the mean phenotype and phenotypic variance in the population, as described in the literature (Lande and Shannon 1996). When considering several selection episodes per generation, there is an additional load that is directly linked to the fluctuations in the optimum (second term in eq. 19). Indeed, this load vanishes when the optimum does not fluctuate across episodes (θ k = θ for all k). It thus corresponds to the deaths of individuals induced by environmental fluctuations within a generation, which reduce population growth without evolutionary significance. Note that when the trait has phenotypic variance, the mean phenotype does change under selection within a generation, tracking the optimum phenotype to some extent following each selection episode. But each of these selection episodes requires a number of selective deaths, and the resulting total cost of natural selection (sensu Haldane 1957) over a generation is unaffected by phenotypic tracking: the load is exactly the same if there is no phenotypic variation, and no within-generation phenotypic tracking. When the strength of selection at each episode ω is constant (corresponding to eq. 5), we further have such that the intra-generational load is proportional to the sample variance V¯[θ] of the optimum across selection episodes within a generation.
In models of selection with a Gaussian fitness peak, there is a simple connection between the temporal distribution of the directional selection gradients and the expected maladaptation affecting populations growth, since both depend on the expected squared mismatch of the mean phenotype from the optimum ( A1). However, note that this decreased variance of the effective optimum come at the expense of an increased load of fluctuating selection within a generation, which cannot be overcome by adaptive tracking.

Discussion
Most theoretical studies investigating the effect of an autocorrelated environment on evolutionary dynamics and their consequences for population dynamics assumed simple life cycles, where the timing of selection within a generation was left implicit (e.g. Charlesworth 1993, Lande and Shannon 1996, Chevin 2013. In models of non-overlapping generations, the timing of life-history events is generally not addressed, and selection is most often modeled as a single, effectively instantaneous process. This contrasts with the fact that many traits are likely to affect fitness at different times in life, and thus undergo several episodes of selection Wade 1984b, Marshall andMorgan 2011). Another difficulty at the interface of theory and empirical work on these questions is that measuring fitness is in practice challenging, such that natural selection is mostly estimated in the form of individual selection episodes acting via specific fitness components in the life cycle (Kingsolver andDiamond 2011, Siepielski et al. 2017). Recent refinements (e.g. McGlothlin 2010) of theory on the measurement of selection gradient (Lande and Arnold 1983, Arnold and Wade 1984a, Wade and Kalisz 1989 investigated how estimates of selection at different episodes within a generation can be integrated to infer overall selection, and the resulting evolutionary response. The present study extends this theory by investigating how fluctuating selection resulting from changes in the environment within and across generations impacts evolutionary trajectories, and the potential for population persistence. To this aim, we developed a model of an organism with non-overlapping generations, where several transitions between stages (survival) prior to reproduction depend on the match between a phenotypic trait and a fluctuating optimal value directly related to the environment, consistent with recent empirical estimates of fluctuating selection in the form of a moving phenotypic optimum (Chevin et al. 2015, Gamelon et al. 2019. Our analysis of this model demonstrates that knowledge of fluctuation patterns for the environment affecting selection is not sufficient to infer the predictability of selection across generations, when there are multiple selection episodes per generation. Selection across generations depends on an integrated optimum that includes the relative contribution of each selection episode to total selection along the life cycle, similar to recent formulations for age-structured populations but without autocorrelation (Engen et al. 2011, Cotto et al. 2019. Therefore, fluctuations of selection across generations depend on the pattern of variance and autocorrelation of the effective optimum, rather than of the environment per se. One of our main results is that, when the environment undergoes red noise (first-order autoregressive process), fluctuating selection within a generation causes the autocorrelation of the effective Cotto and Chevin Page 13 Theor Popul Biol. Author manuscript; available in PMC 2020 October 01. optimum to be inflated relative to the autocorrelation of the environment. We show that as the number of selection episode per generation increases, the inflation factor rapidly converges to an asymptotic value that depends neither on the number of episodes of selection nor on the strength of selection, but only on the autocorrelation of the environment over a time window of selection. This asymptotic value can be used in practical situations to estimate the autocorrelation of the effective optimum directly from measurement of the autocorrelation of the environment. Importantly, however, this result depends on the assumption that total selection (or its expectation) does not depend on the number of selection episodes per generation. This assumption allows us to disentangle the effect of selection acting repetitively along a generation from the effect of stronger total selection. We emphasize that this independence may not be satisfied in natural populations, where more selection episodes may often cause stronger overall selection.
The inflation factor for the optimum directly influences the autocorrelation of directional selection gradients, which relates more directly to empirical measurements of the predictability and consistency of directional selection through time Grant 2002, Morrissey andHadfield 2012). However, an important point when considering the connection with empirical work is that, instead of tracking phenotypic change under selection episode per episode as done empirically Wade 1984b, McGlothlin 2010), we combined all episodes to derive a total per-generation individual fitness function. This allowed us to elude the difficulties arising from the dynamics of variances and covariances of breeding and phenotypic values at each successive episode of selection (see Supp. Mat. for the calculations per-episode). Indeed, the quantitative genetics framework assumes that genetic and environmental values are independent prior to selection. It is no longer the case after the first episode of selection within a generation (prior to reproduction), which generates a covariance between the above components. Tracking changes in the mean phenotype and breeding values across selection episodes therefore requires taking into account the change in covariance between breeding value and environmental values (See. Supp. Mat. for a demonstration). The buildup of this covariance complicates the investigation at the scale of individual selection episodes, from both a theoretical and empirical point of view (Wade and Kalisz 1989, Cam et al. 2002, McGlothlin 2010, because the selection gradients are no longer sufficient to describe accurately the change in the mean breeding value.
We also found that allowing for multiple selection episodes complicates the connection between environmental autocorrelation and the variance of selection gradients, which is directly linked to the expected maladaptation load reducing population growth (Chevin et al. 2017). With multiple selection episodes, both these metrics vary non-monotonically with the autocorrelation of the environment, with a maximum at an intermediate value. This contrasts with previous results from models with a single episode of selection per generation, where the expected maladaptation load decreases monotonically with increasing autocorrelation of the environment, because this allows for closer adaptive tracking of the moving optimum via genetic evolution (Chevin 2013, Chevin et al. 2015. That this may not hold with multiple episodes of selection crucially revises previous understanding of the role of environmental autocorrelation on maladaptation: maladaptation can actually increase with increasing autocorrelation of the environment. Lastly, we found that selection that fluctuates within a generation results in a specific load, which adds up to the maladaptation and variance loads classically described in the literature on adaptation to changing environments (see Lynch and Lande 1993, Burger and Lynch 1995, Lande and Shannon 1996, Kopp and Matuszewski 2014. This load only exists when the direction and/or strength of selection varies within a generation, and also increases with the number of selection episodes and the total per-generation strength of selection. The reason for this load is that each selection episode results in selective deaths that reduce the demographic output of the population, as previously discussed in extensions of Haldane's (1957) famous "cost of natural selection" to the context of adaptation to a changing environment, with a single selection episode (Nunney 2003). Within-generation fluctuations in selection increase this number of selective deaths without any net beneficial effect on adaptation: from an evolutionary perspective, individuals are essentially wasted by withingeneration fluctuation in selection. An important implication of this result is that disturbances that increase temporal variation in the environment can lead to the decline of a population, regardless of its phenotypic state (mean phenotype, additive genetic variance). This bears particular importance in light of the increasing frequency of extreme climatic events (Coumou and Rahmstorf 2012), and accumulating evidence of their effect on evolutionary and demographic dynamics (Moreno andMøller 2011, Marrot et al. 2017). This load caused by fluctuating selection within a generation may only be alleviated by mechanisms that cause within-generation adaptive changes in the mean phenotype without requiring any selective death. This suggests that a powerful mechanism by which labile, reversible plasticity can evolve (as investigated theoretically in a few studies, Gabriel 2005, Lande 2014, Ratikainen and Kokko 2019) is by reducing the load caused by fluctuating selection within a generation.
Our non-overlapping generation model applies directly to univoltine species. An interesting extension would be to investigate the case with overlapping generations. The most closely related life cycle with overlapping generations is that of perennial semelparous organisms, encompassing a large diversity of organisms (Young and Augspurger 1991) spanning from microbial parasites (Ebert and Weisser 1997) to long lived plants (e.g. semelparous yucca) and vertebrates (e.g. pacific salmon). This type of life cycle can lead to asynchrony in reproduction among age classes (Caswell 2001 p.81-88), such that the population at any time is composed of cohorts with different ages that rarely interbreed, in which case our conclusions directly apply to each of these cohorts. Perennial semelparous life cycles can also produce instability in age structure, which has been proposed to contribute to periodicity in species such as cicadas (Bulmer 1977) or bamboos (Keeley and Bond 1999). This instability may even lead to a collapse of age structure and to a synchrony in reproduction (Mjølhus et al. 2005), in which case our model would apply to the entire population. In contrast, iteroparity in perennial organisms would add considerable complications to the investigation of the effects of fluctuations in the environment. Indeed, each cohort of offspring is contributed from cohorts of reproductive individuals having experienced different history of selection, possibly generating deviations from the stable (st)age structure (see Lorimer 1980 for an example). We leave such an investigation for a future study. Cotto and Chevin Page 15 Theor Popul Biol. Author manuscript; available in PMC 2020 October 01.

Figure 1.
Schematic representation of the timescales and autocorrelation measures. Generations last T units of time, within which selection occurs in n selection episodes spread over αT units of time. The autocorrelation of the optimum phenotype (which is identical to that of the environment) over the within-generation time of selection is ψ = ρ αT , and over a generation time is ρ T . The per-generation effective optimum integrates the effect of selection within generations. We denote ρ θtot the autocorrelation of this measure over a generation. The schema highlights that the variance and autocorrelation differ between the optimum phenotype at each episode of selection (i.e. the environment) and the effective optimum relevant for the evolutionary dynamics across generations. Cotto and Chevin Page 19 Theor Popul Biol. Author manuscript; available in PMC 2020 October 01. Theor Popul Biol. Author manuscript; available in PMC 2020 October 01. expectation assuming constant selection strength with same mean (full lines: exact prediction using eq. 12, dashed lines: prediction for continuous selection along each generation -infinite number of selection episodes, eq. 15). Numerical simulations are performed using eq. (10), with θ tot and ω tot 2 calculated as described in eq. 3 from phenotypic optima θ i and corresponding strengthes of selection ω i 2 corresponding to each episode of selection i per generation as described in the following: for each episode of selection, a phenotypic optimum is drawn from equation 11 and the strength of selection is drawn independently in a Gaussian with mean ω tot 2 /n -1 and variance 1. Total selection per generation thus follows a non-central χ 2 with mean ω tot 2 and variance 4ω tot 2 -2n. For simulations, we kept the expected total strength of selection constant (independent on the number of selection episode) to isolate the effect of the number of selection episodes from that of the total strength of selection. For A and B, α = 0.8. For all panels: light gray to black: ρ T = 0.2, 0.5 and 0.9 respectively, ω tot 2 = 25, σ θ 2 = 1.

Cotto and Chevin Page 21
Theor Popul Biol. Author manuscript; available in PMC 2020 October 01. The dashed lines correspond to the predictions from eq. (16). Continuous lines represent predictions derived from a continuous time approximation as in Lande and Shannon (1996, eq. 7), resulting in equation 5a in Chevin and Haller (2014): V ar β tot = V ar θ tot S tot 2 / 1 + S tot G 0 T θ , where T θ is the characteristic time scale of the autocorrelation of the environment (from eq. 13, and see eq. 11 and below). For both panels: α = 1 and σ θ 2 = 1.

Cotto and Chevin Page 22
Theor Popul Biol. Author manuscript; available in PMC 2020 October 01.   (18). Continuous lines: expectation for a single selection episode (eq. 17). Note that for parameters shown here (number of selection episode and strength of selection) there is almost no difference between simulations with deterministic (circles) and random (square) strength of selection per episode. From light gray to black, ρ T = 0.1, 0.5 and 0.9, ω tot 2 = 10, α = 1 and σ θ 2 = 1.

Cotto and Chevin Page 24
Theor Popul Biol. Author manuscript; available in PMC 2020 October 01. Theor Popul Biol. Author manuscript; available in PMC 2020 October 01.

Cotto and Chevin Page 26
Theor Popul Biol. Author manuscript; available in PMC 2020 October 01.