The effect of temporal environmental autocorrelation on eco-evolutionary dynamics across life histories

,


INTRODUCTION
Stochastic variation in the environment is one of the most important drivers of population dynamics, that is, changes in population size and structure (Tuljapurkar 1990, Morris et al. 2008).This variation has usually been modeled with the assumption that consecutive values of environmental states are independent (Ruokolainen et al. 2009, Buckley et al. 2010).However, natural environmental conditions are often temporally autocorrelated (Halley 1996), and the sign of this autocorrelation can determine the establishment of invasive species (Cuddington and Hastings 2016), development of different lifehistory strategies (Metcalf and Koons 2007), and the persistence of phenotypes and populations (Morris andDoak 2002, Wieczynski andVasseur 2016).For instance, negative autocorrelation has been shown to delay extinction compared to random fluctuations as years in which population size declines are more likely to be followed by a year where population size increases, which allows populations to recover (Morris and Doak 2002).For positive autocorrelation, a good year is more likely to be followed by another, and vice versa, leading to larger stretches of population growth and shrinkage.Long stretches of adverse conditions can lead to higher risk of population extinction (Petchey 2000, Laakso et al. 2003, Pike et al. 2004, Tuljapurkar and Haridas 2006, Engen et al. 2013, Chevin et al. 2017), especially in populations with weak density feedbacks (Morris and Doak 2002).Under global environmental change, which is expected to alter temporal autocorrelation in natural environments (Boulton andLenton 2015, Lenton et al. 2017), extinction risk may be exacerbated (van der Bolt et al. 2018).Therefore, understanding the myriad of pathways through which temporal autocorrelation drives population persistence is becoming increasingly important, but these pathways are typically not well integrated into a common framework (Ruokolainen et al. 2009).
Temporal autocorrelation in environmental states affects populations through both evolutionary and life-history processes, which are rarely assessed jointly.Shifts in autocorrelation can change selection pressures and alter the genetic variation among individuals (Ranta et al. 2008), leading to changes in phenotypic trait distributions within populations (Michel et al. 2014) that can affect population fitness (Pelletier et al. 2007, Tuljapurkar et al. 2009, Ozgul et al. 2010, 2012).For example, the prevalence of similar environmental conditions under positive autocorrelation allows for adaptation to these conditions (Lande and Shannon 1996).This may lead to maladaptation when individuals are exposed to another environment (Frank et al. 2017).This (mal)adaptation can, in turn, affect vital rates such as survival and reproduction and therefore lead to changes in the mean population structure and growth (Chevin et al. 2017).At the same time, it is not well known to which degree traitmediated demographic effects of autocorrelated environmental conditions are determined by the life-history strategy of a species.Different strategies, such as the pace of life (slow-fast continuum) or reproductive strategy (lifetime reproductive output; Salguero-G omez et al. 2016), are characterized by different sensitives of vital rates to environmental effects (Stearns 1983, Oli 2004) and are likely to show different potential for adaptation.For instance, species with fast life histories and a high lifetime reproductive output appear to be among the most sensitive to temporal autocorrelation, in terms of population growth (Paniw et al. 2018).These life histories are expected to evolve faster than species with slow life histories under changes in environmental variation (Reed et al. 2011).
Aside from differences in adaptive potential among life histories, differences in phenotypic plasticity can also affect their relative responses to temporal autocorrelation (Lynch et al. 1998).By altering their phenotype, for example, fur color or body size, in response to an environmental perturbation, individuals can respond to novel conditions without genetic change.This adaptive phenotypic plasticity is likely to decrease the effects of autocorrelation as it buffers the effects of natural selection and can make evolutionary changes redundant (Thompson 1991, Oostra et al. 2018).It has been suggested that adaptive phenotypic plasticity plays a more important role for slow life histories, due to their decreased evolutionary potential (Vedder et al. 2013).Conversely, slow life histories may also be more sensitive to nonadaptive phenotypic plasticity, in which phenotypes are expressed that are further away from the environmental optimum, thereby decreasing the fitness of an individual (Ghalambor et al. 2015).One of the reasons that this nonadaptive plasticity can occur is due to the delayed nature of most plastic responses (Wennersten and Forsman 2012).These time lags may occur, for instance, when the phenotype of an individual is determined by environmental cues during early development, but these cues change later in life resulting in a maladaptive plastic response (Helle et al. 2012, Bateson et al. 2014).Time lags are especially important under different autocorrelation levels as environmental cues under which phenotypes develop may prove to be unreliable under increasingly negative autocorrelation, leading to maladaptive responses (Reed et al. 2010).This unreliability can decrease fitness of phenotypically plastic individuals and can lead to an acceleration in evolutionary adaptation (Ghalambor et al. 2007, Schmid andGuillaume 2017).Nevertheless, populations might still suffer if the genetic evolution cannot fully compensate the nonadaptive plastic response.
Understanding how different forms of plasticity interact with temporal autocorrelations for different life histories can help us further explain the differences in the effects of autocorrelation on population dynamics across taxa.While selection pressure, genetic variation, phenotypic plasticity, and demographic sensitivities have been investigated in autocorrelated environment (Ruokolainen et al. 2009, Ezard et al. 2014), no studies have integrated all of these aspects into a single model in order to determine the differences between life histories in evolutionary and demographic responses to different autocorrelation patterns.
Here, we used a novel approach, combining an evolutionary explicit model with demographic models to investigate how temporally autocorrelated changes in the environment simultaneously affect the distribution of a phenotypic trait, population fitness, and population size in different life histories.We used a modified version of the individual-based modeling framework Nemo (Guillaume andRougemont 2006, Cotto et al. 2017).This framework allowed us to simulate the evolutionary dynamics of a single quantitative trait, while also including a demographic model in the form of a matrix population model that represents vital rates of different life histories.We simulated autocorrelated amongyear environmental fluctuations by shifting the optimum trait value.This environmental patterning then changed survival patterns of different life-history stages.Furthermore, to investigate how phenotypic plasticity influences the response to autocorrelated environments, we included both (non)adaptive and time-lagged (non)adaptive phenotypic plasticity in trait responses to temporal autocorrelation.

MATERIALS AND METHODS
For our simulations, we used version 2.3 of the individual-based, genetically explicit, evolutionary modeling program Nemo (Guillaume and Rougemont 2006), which was modified to include age structure and phenotypic plasticity (Cotto et al. 2017).We simulated demographic and trait responses to environmental autocorrelation for a total of 150 life histories.

Life histories
We simulated different life histories using a two-stage matrix population model.The life-cycle stages corresponded to juveniles and adults.Transitions between these stages were determined by three vital rates: juvenile survival (S 1 ), adult survival (S 2 ), juvenile transition to adults (G 1 ), and fecundity (F 1 ).The different life histories were then created by perturbing S 1 , S 2 , and G 1 independently and solving F 1 for an asymptotic growth rate k = 1.2, creating a set of viable life histories with similar population growth in the absence of density dependence (Fig. 1A).Choosing k = 1.2 instead of k = 1, which defines a stable population not regulated by density dependence, ensured that populations did not go extinct instantaneously due to selection and regulation (see Density-dependent regulation).S 1 and G 1 were set to 0.11, 0.33, 0.55, 0.77, and 0.99, while S 2 was set to 0.0, 0.25, 0.5, 0.75, 0.9, and 0.99.We used all possible combinations between values of S 1 , S 2 , and G 1 to create 150 stage-structured populations with life histories.All populations were sexually reproducing with random mating between adult male and female individuals with a 1:1 sex ratio in offspring.
For all of the matrix models, we calculated five commonly used life-history traits: generation time, net reproductive rate, age at sexual reproduction, degree of iteroparity, and annual sexual reproduction (Caswell 2001, Salguero-G omez et al. 2016, Paniw et al. 2018, Appendix S1: Fig. S1).The populations ranged, for instance, from 2 to 22 yr in generation time and from almost completely semelparous to iteroparous when assuming asymptotic dynamics (i.e., not considering density dependence).Reproductive rate was positively, nonlinearly related to generation time, as net reproductive rates and annual reproduction were constrained in our simulations due to the constrain on F 1 .Iteroparity varied independently from generation time and reproductive rate (Appendix S1: Fig. S1).
The matrix models for each of the 150 life histories were then used as the demographic model in the Nemo simulations (Fig. 1A).In this simulation, all individuals followed a life cycle (Fig. 1B) in which they were first subjected to a density-dependent survival based on a simple Beverton-Holt model (Fig. 1C).Conditional on surviving density dependence, individuals were subjected to selection in which survival was based on an individual's phenotypic trait and the environment (Fig. 1D).Individuals that survived selection then bred with the probability to produce offspring based on the demographic matrix models (Fig. 1E).Finally, individuals aged, that C Fig. 1.Individuals (points) transition through life-cycle stages (different colors) from time t to time t + 1 and change associated phenotypic trait values according to the experimental setup of Nemo.In this setup, demographic matrix models describe the different life histories (A); and each individual must pass through five life-cycle events (ovals) to progress to the next time step (B).Within these events, density-dependent regulation (C) is based on the Beverton-Holt equation, where the chance of surviving to the next time step is dependent on the total number of individuals (N tot ) and a competition parameter M; individuals that pass density-dependent regulations are then selected based on their phenotypic trait value, which determines survival and can take on two optima depending on the environmental condition (D).At each time step, one of the two conditions, determining selection on the trait, is drawn depending on the autocorrelation.After selection, individuals breed, which corresponds to drawing a number of newborns randomly from a Poisson distribution with a mean of F 1 (E).Lastly, individuals age following the structure of the matrix population model in A, which was modified to include the newborns (F).
is, transitioned stages based on the matrix parameters (Fig. 1F).

Density-dependent regulation
When projecting the dynamics of each simulated population through time, we first exposed all individuals in a population to a densitydependent regulation (Fig. 1B).This regulation mimicked within-species competition to limit the maximum population size and therefore simulate biologically realistic population dynamics (Cotto et al. 2017).The regulation occurred via a simple Beverton-Holt model, where a competition-based survival probability (S bev ) was calculated using the total number of juveniles and adults (N tot ) and the weight of each individual in the population set at M = 0.0001 (S bev = 1/ (1 + N tot 9 M)) (Beverton and Holt 1957; Fig. 1B).Juvenile and adult individuals then survived based on S bev .As density dependence acted on both life stages, population sizes after density-dependent survival were the same across life histories.

Individuals' genetic composition and phenotype expression
Each individual in a simulated population was characterized by 10 unlinked genetic loci (Fig. 1C) that additively contributed to a single quantitative trait g 0 .We used a continuum-of-alleles model of mutation which assumes that the distribution of breeding values at each locus is Gaussian.Mutational effects were drawn from a normal distribution with a = 0 and r 2 = 0.05 and a mutation rate of m = 0.001 (Appendix S1: Fig. S2).The phenotype z of an individual was determined by a linear reaction norm z = g 0 + g 1 e with g 0 as the (evolving) quantitative trait of the individual and slope g 1 as the degree of plasticity attributed to the phenotypic optimum e similar to Lande (2009).The absence of phenotypic plasticity was achieved with a constant reaction norm slope of g 1 = 0 such that the phenotype expression was independent of the phenotypic optimum e and solely depended on g 0 .We simulated constant plasticity such that g 1 was not evolving over time.
Adaptive and nonadaptive phenotypic plasticity were installed with g 1 6 ¼ 0, ranging from À0.5 to 0.5.Positive g 1 values indicate adaptive plasticity while negative indicate nonadaptive plasticity.Phenotypic plasticity was implemented in two different ways: one in which the phenotypic optimum of the current year e t was perceived during development and one in which the phenotypic optimum corresponded to the environment of the previous year e (tÀ1) allowing for time-lagged effects.

Selection
Our environment could fluctuate between two environmental states, characterized by phenotypic optimums (e) of e = +1 and e = À1.The phenotypic optimum described the phenotypic values with a survival probability of 1 (i.e., maximum fitness at selection; Fig. 1D).This fitness declined when moving away from the optimum following a Gaussian distribution with variance = x 2 .This variance depicted the selection strength.During each iteration of a simulation, the environmental condition and therefore the phenotypic optimum were sampled from a sequence generated using a two-state Markov chain: In this chain, p x and p y are the probabilities of transitioning to the environmental state with e = À1 at year t + 1 when the environment is e = À1 or e = 1, respectively, at year t.While keeping the mean of e = 0, we created autocorrelated sequences of the two environmental states by varying p x and p y from 0.1 to 0.9 (as in Paniw et al. 2018).This resulted in an autocorrelation range from a negative autocorrelation of À0.8 to a positive autocorrelation of 0.8 with 0.2 increments.
Each year, the fitness of each individual that survived density regulation was calculated with a Gaussian fitness function with mean l (i.e., the phenotypic optimum e; see Individuals' genetic composition and phenotype expression) and variance x 2 = 5.Fitness corresponded to surviving selection, where survival = 1 at e for both juveniles and adults, and declined following the Gaussian fitness distribution (Fig. 1C).Fitness declined with increasing distance between z and l, corresponding to a survival rate of 0.63 when e corresponded to the opposite phenotypic optimum (i.e., z = À1 when e = 1).In a separate set of simulations, selection was applied to a single life stage only (juvenile or adult) together with an increased selection strength of x 2 = 3 (see Appendix S1: Fig. S3 for different variances of the optimum curve).

Breeding and aging
Sexual reproduction was implemented as the number of offspring per female that survived selection using a Poisson distribution with l = 2 9 F 1 .This allowed us to include individual stochasticity in reproduction.For each offspring of a female, a male was selected randomly to provide genetic material.These offspring were stored in a separate stage until the life-cycle event aging.In this event, all offspring transitioned to the juvenile stage.The inclusion of an offspring stage is a modification to Nemo which ensures that all of the offspring recruited in the breeding step end up in the juvenile stage at t right after aging.In the absence of this modification, individuals would be able to be born into the adult life stage and thus our simulation would not represent the two-stage life-cycle transitions depicted in the demographic model.Juveniles could then remain juveniles in next year with a probability of S 1 9 (1 À G 1 ) or transition to adult with a probability of S 1 9 G 1 .Adults either survived (S 2 ) and remained adults or died.

Simulations
We ran four different types of simulations on all life history 9 autocorrelation combinations: (1) Individuals showed no phenotypic plasticity, and both juvenile and adult survival were under selection; (2) no plasticity, while selection acted on either juvenile or adult survival; (3) individuals showed five different levels of phenotypic plasticity, and selection acted on both adults and juveniles; or (4) time-lagged phenotypic plasticity in the trait, selection on juveniles and adults.Every year, the following sequence of life-cycle events occurred in the simulation: regulation (removal of individuals based on density dependence), selection (survival of individuals based on the expressed phenotype and the current environmental condition/phenotypic optimum), breeding (sexual reproduction), and aging (stage transitions based on the MPM; Fig. 1).To ensure all projections of trait and population dynamics started from a stable state and to eliminate any unwanted transient dynamics, a burn-in of 7000 time steps was run for each combination of the autocorrelation level, life history, and simulation type.Simulations began with the population structure at the end of the burn-in and lasted for another 300 time steps.For each different life history 9 autocorrelation combination, we generated five different sets of environmental conditions using the Markov chain, and for each set, population projections were replicated five times.

Analysis of simulation results
For each replicate simulation run, the coefficient of variation (CV) of the total population size and mean population fitness were calculated.For the simulations including phenotypic plasticity, we also calculated the extinction probability as here population extinctions occurred due to the negative fitness effects of maladaptive phenotypic plasticity (Appendix S1: Figs.S5, S6).We used CV as a metric of population and fitness variation as it accounts for possible differences in population sizes.Since populations with higher variability and lower means are more likely to reach critical population size and this variability can occur independent of the mean, CV can be used as an important indicator of population persistence (Trotter et al. 2013).We tested for the effect of generation time and autocorrelation and their interaction on the CV of population size, fitness, and the mean additive genetic variance within populations using generalized linear models (GLMs).Generation time, which correlates strongly with the fast-slow continuum, was used to depict variation in responses among lifehistory strategies (Lebreton andClobert 1991, Smallegange et al. 2014).We used an inverse gamma GLM to fit the model, as our data were only positive and continuous.In addition, we also tracked changes in mean population size across simulation scenarios (Appendix S1: Fig. S4).Average population sizes were approximately normally distributed, and we therefore used a Gaussian error distribution to model them.Extinction probability was modeled using a quasibinomial distribution.We fit a full model, including main effects and interaction of predictors to estimate effect sizes (Appendix S1 for model coefficients).We removed a life history with a generation time of 2.35, as it showed large between-year variation in population size.This occurred due to high juvenile transition rates and low adult survival leading to a 2-yr cycle in which a year either had a very small or very large population size.These cycles were not related to selection but pure artifacts of parameterizing the matrix model.Removing them did not alter the results.
All data analysis was carried out in R 3.4.2(R Core Team 2017).

Variation in population fitness and size in the absence of plasticity
When selection acted on adults and juveniles, both the temporal patterning in the environment and the life-history strategy of the simulated populations determined the CV in mean population fitness.Positive autocorrelation substantially increased the CV of mean fitness (measured as survival probability during selection; Appendix S1: Table S1).Slow life histories with a long generation time had a lower CV in mean fitness compared to fast life histories with shorter generation times across autocorrelations (Fig. 2A; Appendix S1: Table S1).The CV was highest at positive autocorrelation and was increased more for long-lived populations than for short-lived ones, although variation remained highest for fast life histories.
The effects of autocorrelation and life history on variation in fitness carried over to population dynamics.At higher positive autocorrelations, the CV of population size increased for both fast and slow life histories (Fig. 2B; Appendix S1: Table S2).Furthermore, slow life histories had a lower amount of additive genetic variance compared to fast life histories (Fig. 2C; Appendix S1: Table S3).For the fast life histories, we also found lower mean population fitness and mean population size (Appendix S1: Fig. S4).Variation between replicates was minimal, and the patterns described above occurred in all replicates.
The relationships between temporal autocorrelation and population-size variation across life histories changed when selection affected only one of the stages.Fast life histories showed more variation in population sizes when selection only acted on juveniles at positive autocorrelations (Fig. 3; Appendix S1: Table S4).The opposite was true for slow life histories, which showed more variation at positive autocorrelation when selection only affected the adults (Fig. 3).

Effects of adaptive plasticity on trait and population dynamics
When including adaptive plasticity, autocorrelation displayed a smaller effect on variation in population size compared to simulations omitting plasticity.Nonadaptive plasticity, on the other hand, amplified the effect of autocorrelation on CV in fitness and population size.Under nonadaptive plasticity, population extinction rate decreased with generation time (Fig. 4; Appendix S1: Table S5).When introducing a time lag, the variation in population dynamics in response to autocorrelation changed compared to non-lagged plasticity.For negative autocorrelations, time-lagged adaptive phenotypic plasticity increased the CV of population size (Fig. 4) and increased the probability of extinction (Appendix S1: Fig. S5, Table S6) across fast and slow life histories.Nonadaptive phenotypic plasticity had the exact opposite effect as it increased extinction and CV of population size at positive autocorrelation for fast and slow life histories alike (Fig. 4; Appendix S1: Fig. S6).

DISCUSSION
Population dynamics of species with different life histories can respond differently to changes in the temporal autocorrelation of the environment (Paniw et al. 2018).Moreover, temporal autocorrelation affects genetic variation and selection strength (Ranta et al. 2008, Michel et al. 2014).Our results show that life histories differ in simultaneous population and genetic-trait responses to changes in autocorrelation patterns.Fast life histories had larger genetic variation in their populations compared to slow life histories independently of the level of environmental autocorrelation they were exposed to.This was due to the higher annual reproductive rate of the faster life histories, which allowed them to adapt more rapidly to changes in environmental optima than slow life histories.When the trait mean of a population was closer to the environmental optimum, the mean population fitness, here measured as survival, increased.However, due to the Gaussian shape of the optimum curves, which allowed quick adaptation to one optimum, shifts in the optimum led to equally quick decreases in the mean population fitness and mean population sizes.On the other hand, populations of slow life histories, which adapted slower to an environmental optimum due to their lower genetic variation, experienced smaller differences in fitness compared to fast life histories when the optima shifted, and consequently fluctuated less.These patterns became weaker at more positive autocorrelation as there was more time to adapt, leading to a larger increase of variation in fitness and population size for the slow life histories compared to the fast life histories, although the latter still had higher variation.These results point to potential implications of omitting the effect of autocorrelated environments when assessing the performance of species under environmental change.Previous studies suggested that positive autocorrelation can play an important role in the successful establishment of invasive species (Cuddington and Hastings 2016), which tend to have relatively short generation times (Salguero-G omez et al. 2016).We provide a mechanistic understanding of the evolutionary and demographic drivers that underlie the success of some species under environmental change.

Selection
Our study provides a better understanding of the evolutionary processes contributing to demographic buffering.Long-lived species have been predicted to be more sensitive to changes in the environment due to lower adaptability (Vedder et al. 2013).However, our results, showing lower variation in population size for slow life histories, were more in line with studies showing demographic buffering against temporal autocorrelation in populations with long generation times and low reproduction (Morris et al. 2010, Engen et al. 2013, Paniw et al. 2018).In an evolutionary sense, slower life histories may buffer populations from extinction because more genetic variation can be maintained in adult individuals (Hailer et al. 2006).When selection acted on just the juveniles, our results seem to corroborate this.However, when selection acted on both juveniles and adults our results show the opposite, as there was lower genetic variation present in slow life histories compared to fast life histories.This may occur when reproduction is constrained (as imposed by our models), resulting in the birth of more individuals in each time step for fast species leading to more frequent mutation and recombination events, while the slower life histories are unable to maintain genetic variation in the adult life stage.Aside from the differences in genetic variation between life histories, differences in demographic sensitivity also play a role in the response to selection.This is because longer lived species with low annual fecundity are more sensitive to changes in adult survival compared to faster species (Heppell et al. 2000).This increased sensitivity leads to a stronger selection pressure, leading to a larger variation in the population size because the average phenotype more strongly follows the environmental conditions.Faster life histories are more sensitive to changes in juvenile survival and fecundity.Since we did not perturb fecundity, we only find minor effects of autocorrelation when selection acts upon juvenile survival.These results indicate that understanding the way selection acts upon populations is important in both an evolutionary and demographic sense to take suitable conservation measures.

Phenotypic plasticity
Our simulations provide novel insights on the effects of the interaction between environmental patterning and different forms of plasticity on trait and population dynamics.One of the mechanisms species use to cope with environmental change is phenotypic plasticity (Gienapp et al. 2008, Hoffmann and Sgro 2011, Meril€ a and Hendry 2014).For example, great tits (Parus major) have been shown to plasticly change their behavior allowing closer tracking of the environment (Charmantier et al. 2008).In our study, immediate adaptive plasticity (i.e., occurring at the same time as the environment changed states) buffered against changes in the environment.Population fitness increased leading to decreased differences in fitness and population sizes between autocorrelation levels across life histories.This conflicts theory predicting that phenotypic plasticity is more important for slow life histories due to their inability to genetically adapt quickly enough to changing environments (Vedder et al. 2013, Camarero andFajardo 2017).In fast life histories, plasticity may be equally important as the rate of environmental change may outpace the adaptive potential of these life histories.However, these findings need to be corroborated in future empirical studies.
Under time-lagged phenotypic plasticity (i.e., individuals respond plastically to environmental conditions at time t À 1), the influences of autocorrelation on variation in population size became more apparent.At more positive autocorrelations, conditions remain more stable so time-lagged phenotypic plasticity is mostly adaptive.This corresponds to other theoretical studies showing that plastic responses lagging behind environmental cues benefit from more positive autocorrelated and predictable environments (Reed et al. 2010, Fischer et al. 2011, Ezard et al. 2014, Ashander et al. 2016).Our results therefore highlight that increases in autocorrelation under environmental change can prove advantageous for traits that are affected by timelagged phenotypic plasticity.

Implications and future research
The individual-based model we used effectively combined trait dynamics with population dynamics to assess complex relationship between evolutionary processes and life histories in autocorrelated environments.Using this framework, important eco-evolutionary questions can be addressed, for example, on differences in adaptability between species, evolution of plasticity, and the role of life-history strategies and evolutionary feedbacks in invasion processes.Our approach substantially advances ❖ www.esajournals.orgprevious studies describing the effects of autocorrelation on population dynamics across life histories (Paniw et al. 2018) by integrating phenotypic plasticity into a joint evolutionary-demographic model, the latter accounting for density dependence.Density dependence has been shown to influence population responses to autocorrelation (Levine andRees 2004, Engen et al. 2013).Here, we used a simple density model which affected all life histories in a similar fashion.That is, both juveniles and adults in each population experienced negative density dependence of equal strength, which was implemented by randomly culling a certain percentage of the population.However, this does not necessarily reflect the conditions in natural populations.For instance, in plants density dependence usually affects recruitment of seedlings (Comita et al. 2014).In animals, density dependence can take on many forms.For example, density dependence in the form of intraspecific competition affects juvenile survival or growth in numerous species, including wolves (Cubaynes et al. 2014), great tits (Reed et al. 2013), and mackerel (Jansen and Burns 2015).In other species, for example, cormorants, adult survival is mostly affected by population density (Frederiksen and Bregnballe 2000).Further investigation into how density dependence influences the adaptability of different life histories across taxa can improve our understanding of the eco-evolutionary dynamics.
Aside from assumptions of density dependence, we described phenotypic plasticity using a linear reaction norm, in which all life histories were equally affected by phenotypic plasticity.However, plasticity can also be modeled as a trait that can evolve in response to environmental variability.Since there are costs associated with plastic responses, different patterns of environmental change can lead to fitness advantages and disadvantages (DeWitt et al. 1998, Pigliucci 2005).Describing the evolution of phenotypic plasticity may provide further insight into how different life histories respond to stochastic environments (Franks et al. 2014, Camarero andFajardo 2017).
Lastly, although the two-stage matrix population models we developed could generate a wide range of life histories and have been used in numerous comparative studies (Koons et al. 2009(Koons et al. , 2016)), including more complex life histories will provide more detailed insights into evolutionary dynamics.For example, species with a more complex life cycle such as most herbaceous plants may have seed banks that can act as a storage for genetic material, allowing for better adaptation to changes in the environment (Ellner and Hairston 1994, Sasaki and Ellner 1997, Tanksley and McCouch 1997).This can help reduce the fluctuations in population size due to autocorrelated environmental change (Nunney and Nunney 2002).

CONCLUSIONS
Predicted global changes in environmental patterning can have a large influence on the evolutionary and demographic processes across life histories.Combining life-history differences in adaptability to new conditions with stage-or age-specific differences in sensitivity to selection will allow for better estimations of population responses to environmental change.We emphasize that phenotypic plasticity can play an important role in mediating these responses by either amplifying the amount of genetic and population-level variation in the case of nonadaptive plasticity or buffering in the case of adaptive plasticity.Therefore, we suggest that further research is needed which incorporates differences in the evolution of plasticity between life histories.

Fig. 2 .
Fig. 2. Temporal autocorrelation in trait selection optima alters the coefficient of variation of the mean population fitness (measured as survival probability under selection) (A) and of population size (B), and the additive genetic variance (C) for the different life histories described using generation time.These metrics were calculated over 300 generations and averaged over 25 simulation replicates.Points and lines indicate simulated valued and average predictions from generalized mixed models, respectively.

Fig. 3 .
Fig.3.Selection acting on two different stages (juveniles and adults) separately differently affects the coefficient of variation of the mean population size for the different life histories, characterized by generation time, and temporal autocorrelations in trait selection optima.These metrics were calculated over 300 generations and averaged over 25 simulation replicates.Points and lines indicate simulated valued and average predictions from generalized mixed models, respectively.

Fig. 4 .
Fig.4.The effect of immediate (no time lag) and time-lagged phenotypic plasticity on the coefficient of variation of the mean population size across temporal autocorrelation patterns and for the different life histories described using generation time.Facets indicate different levels of adaptiveness with g 1 = À0.5 being strongly maladaptive, and g 1 = 0.5 being strongly adaptive, as described by the reaction norm function (z = g 0 + g 1 e).Only simulations in which the population did not go extinct during the burn-in are shown.These metrics were calculated over 300 generations and averaged over 25 simulation replicates.Points and lines indicate simulated valued and average predictions from generalized mixed models, respectively.