The inclusion of immediate and lagged climate responses amplifies the effect of climate autocorrelation on long‐term growth rate of populations

Climate variability will increase with climate change, and thus it is important for population ecologists to understand its consequences for population dynamics. Four components are known to mediate the consequences of climate variability: the magnitude of climate variability, the effect size of climate on vital rates, covariance between vital rates and autocorrelation in climate. Recent studies have pointed to a potential fifth component: vital rates responding to climate in different timeframes, with some responding more immediately and some having lagged responses. We use simulations to quantify how all five components modify the consequences of climatic variability on long‐term population growth rates across a range of life histories defined by life expectancy and iteroparity. We use an established method to compose Matrix Population Models for 147 life histories. Our simulations show that including different timeframes for vital rates responses to climate can either reduce or amplify the negative influence of climate variability on long‐term population growth rates. The negative effect of different timeframes for vital rates responses on population growth is amplified when climatic autocorrelations are negative, and when species are long‐lived. Synthesis. The existing literature shows that vital rates often respond to climate in different timeframes, and that studies often ignore climate autocorrelation. Our results show that simultaneously including both of these factors can substantially increase or decrease a population's expected growth rate. Moreover, the relative magnitude of this change increases with the generation time of a life history. Our results are relevant to conservation, population forecasts and population modelling in general.


| INTRODUC TI ON
In recent years, the threat of climate change to both plant and animal populations has become a central topic in ecology (Clark et al., 2001;Urban et al., 2016). Climate variability is projected to increase in the future (IPCC, 2014), and studies suggest that this could pose a larger threat to populations than changes in mean climate (e.g. Vasseur et al., 2014).
Thus, it is important to understand how climate variability influences the vital rates (survival, reproduction, etc.) of species, their annual population growth rate (λ) and their long-term stochastic population growth rate (λs; Barraquand & Yoccoz, 2013;Lewontin & Cohen, 1969). The literature has examined four components that influence the effect of climate variation to long-term population growth: (1) the magnitude of climate variability (Boyce et al., 2006), (2) the susceptibility of vital rates to climate, in particular vital rates with high sensitivity (Morris et al., 2008), (3) the covariances among vital rates (Iles et al., 2019) and (4) the environmental autocorrelation in the climate (Fey & Wieczynski, 2017).
Recent studies identified a fifth component that could affect long-term population growth: the timeframe in which vital rates respond to climate drivers (Evers et al., 2021). For example, some vital rates might respond almost immediately while others have lagged responses to climate drivers. We do not know how much temporally varied responses (TVRs; i.e. the fifth component) influence long-term population growth.
It is also unclear what the relative effect and importance of the first four components are in the presence of TVRs.
Climate variability (component 1) and the susceptibility of a species' vital rates to climate (component 2) play a large role in determining the interannual variation in λ. The higher the interannual variation in λ, the lower the λs (Lewontin & Cohen, 1969;Tuljapurkar, 1990). As a result, populations for which sensitive vital rates (vital rates that strongly influence λ, Caswell, 2001) respond strongly to climate drivers are expected to change the most from increases in climate variance (e.g. Boyce et al., 2006). Covariation among vital rates (component 3) can mediate climatic effects. In particular, positive covariation increases, while negative covariation dampens, interannual variation in λ, and thus decreases and increases λs, respectively (Doak et al., 2005). In the context of climate drivers, positive covariation arises if all vital rates respond in the same direction to a certain climate driver, whereas negative covariation arises when two (or more) vital rates respond in opposite directions to the same climate driver.
The environmental autocorrelation in climate (component 4), in which the climate at each point in time is correlated to the previous environment, also influences interannual variation in λ and extinction risk. Positive environmental autocorrelation tends to increase extinction risk because populations in decline tend to stay in decline. On the other hand, negative environmental autocorrelation tends to stabilize populations, as declines are followed by increases (Heino & Sabadell, 2003;Pilowsky & Dahlgren, 2020;Schwager et al., 2006). The effect of environmental autocorrelation generally has a small effect on population growth rates when compared to the effects of other components, such as vital rate covariation (e.g. Morris et al., 2011). However, environmental autocorrelation can be important for species that recover slowly from perturbations , such as those with long lifespans and high reproductive output (Salguero-Gómez et al., 2016).
The fifth component, when different vital rates respond to climate drivers in different timeframes (TVRs), has been shown to occur by recent research that considers climate timeframes other than the typical first 12 months prior to vital rate responses (Evers et al., 2021;Tei et al., 2017;Tenhumberg et al., 2018). Here, we show that TVR and climate autocorrelation affect λs via their effect on vital rates covariations. For example, consider a species with two vital rates that respond positively to the same climate driver. In the absence of TVR, the correlation between the vital rates will be strongly positive. However, in the presence of TVR (e.g. survival responds to climate in the current year, and fecundity responds to climate in the previous year), the correlation between the vital rates will depend on the temporal autocorrelation of the climatic driver. Specifically, strong negative autocorrelation will produce a strong negative vital rate covariation, thus increasing λs; vice versa for a strong positive autocorrelation.
We also expect that life history of a species will influence the extent to which TVR will influence populations growth rate. Species with low life expectancies typically have low juvenile survivorship, and species with high iteroparity typically have high adult survivorship. When the means of survivorship are close to zero or to one, high coefficients of variation are not possible (Morris & Doak, 2004), and we expect the effects of TVR to diminish in magnitude. This leads us to the expectation that the relative effects of TVR will vary with life history (e.g. longevity and parity).
Here, we use simulations to investigate how the five components we described above mediate the effect of climate variation on λs. We simulate matrix population models that represent a wide range of life histories. We then run stochastic simulations in which TVR is either present or absent, while modifying the first four components (climate variability, climate effect strength, vital rate covariation and climatic environmental autocorrelation). By doing so, we elucidate how longterm viability responds to TVR across a large range of life histories.

| MATERIAL S AND ME THODS
To investigate the effect of TVR on long-term population growth rate across a broad range of life histories, we used a well-established framework to create Matrix Population Models (MPMs) representing a wide range of life histories (Neubert & Caswell, 2000). We use these MPMs to conduct stochastic simulations of their dynamics under different scenarios of environmental autocorrelation, environmental variance, strength of climatic signal, vital rate covariation and TVR ( Figure 1).

| Simulating temporal sequences
We simulated variation in the vital rates of MPMs starting from normally distributed environmental sequences (V) with standard deviation V . These environmental sequences reflect the response of a vital rate to both a climate driver, C, and unexplained environmental F I G U R E 1 Workflow to simulate the effect of climate autocorrelation (r K ), climate variability (standard deviation of its distribution, σ C ), signal strength (p) and temporally variable response (TVR), on the stochastic natural logarithm of population growth rate (λ s ). (a) Create climate sequences of 10,000 steps, with different levels of autocorrelation and σ c , and combine them with random noise into an environmental sequence (in this example, 50% climate signal, 50% noise, signal strength = 0.5). (b) Using a 2 × 2 MPM, create 147 different life histories, with different values for transition probability from juvenile to adult (γ), juvenile and adult survival (S J and S A ), and with fecundity (φ) set so that matrix A produces a stable population (population growth rate = 1). As an example, one life history (ID) can be seen in the table, with vital rate means and standard deviations (in parentheses, for the three fluctuating vital rates). (c) For each time step, here as example i = 10, calculate the quantile probability of the recent and one-step lagged value from the normally distributed temporal sequence from (a), given a mean of 0 and standard deviation of σ c . Using this quantile probability, calculate the corresponding quantiles on the beta (for vital rates S J and S A ) and gamma (for vital rate φ) distributions given the vital rates' respective mean and standard deviation defined in (b) to populate the A i matrix. (d) Repeat these calculations for all steps in the sequence. (e) Calculate the λ s using the sequence of A matrices. (f) Create the result table with λ s for each of the different life histories, autocorrelation, climate variability, signal strengths and simulation type (TVR) as shown here, or control where φ also responds to recent climate (in green). variation represented as random noise, (Figure 1a). We control the environmental variance (V) explained by climate (C) using signal strength (p). Signal strength varies between 0 and 1, where, for example, 0.5 and 1 imply that climate explains, respectively, 50% and 100% of the environmental variance V 2 . We then converted these normally distributed sequences to the beta and gamma distributions that characterize the survival and fecundity rates, respectively (see Section 2.3 below, Figure 1c).
We simulated the environmental sequences (V) by adding two separate random processes: the climate sequence (C) and the unexplained variation ( ). We first simulated climate sequences, C , using 35 combinations of standard deviation and autocorrelation (see S1.1 for detailed methods). We included five levels of the environmental The final step to produce the environmental sequence, V, was to simulate random noise ( ), partition the variance of C and , and add them together. We included random noise in the temporal sequences to represent other factors that influence population dynamics in the real world, such as anthropogenic disturbances, biotic interactions or other unknown climate drivers. We computed each individual value, i, of this temporal sequence, V, as where i is the ith individual random deviate from a normal distribution with mean 0 and standard deviation σ V , C i is the ith random deviate from climate sequence described above with mean 0, standard deviation σ C , and an autocorrelation level. We multiply each random deviate C and i by parameters p and 1−p to control the proportion (p) of variance in V that is explained by the climate driver C; p can also be seen as the signal strength of C, or the susceptibility of the vital rates to the climate driver C. Because our objective is to produce an environmental sequence V with standard deviation σ V , summing up C and with untransformed variance 2 V would produce a V with standard deviation √ 2 2 V . Multiplying each C i and i random deviate by p and 1−p , respectively, shrinks their standard deviation to produce a V variable with the desired σ V . For example, if the signal strength (p) is 0.5 and we aim to produce a random variable V with a standard deviation (σ V ) of 1, p and 1−p are equal to approximately 0.7071. These values make intuitive sense on the variance scale: they produce two random variables C and E with standard deviation 0.7071 (and therefore variance 0.5), which sum to produce a variable V with standard deviation 1 (and therefore variance 1). We implemented four different p values (0.05, 0.25, 0.5 and 1). As such, we have temporal sequences where hardly any variance is explained by the climate driver, to sequences where the temporal sequence is fully driven by the climate driver. Our sequences thus encompass a range of temporal variance, autocorrelation and of variance explained by the climate driver.

| Population models for a range of life histories
To address how TVR (component 5) affects population demography, we used the Matrix Population Model (MPM) parameterization suggested by Neubert and Caswell (2000; Figure 1b). This MPM has two stages, juvenile and adult, and yearly transitions are described by the following equations: where n t and n t+1 are population size vectors at time t and t + 1, respectively, A t is the transition matrix, γ is the probability of transitioning from juvenile to the adult stage if the individual survives to from t to t + 1, S J,t and S A,t represent the survival probability of juveniles and adults, respectively. Finally, φ t is the number of offspring produced per surviving adult. This MPM can model a large range of life histories depending on the vital rate (γ, S J , S A and ) values (Neubert & Caswell, 2000). For example, changing adult survival so that it approaches 0 (S A → 0) changes the model species from iteroparous to semelparous (Neubert & Caswell, 2000). We recognize While evidence contrary to demographic buffering exists (e.g. Jäkäläniemi et al., 2013;McDonald et al., 2017), we decided to follow Koons et al.'s (2016) method in this as well, as it reflects the more common evidence on demographic buffering. To simulate these standard deviations, we used the elasticities (e) of S J , S A and φ to calculate a proportional measure of buffering: We calculated the standard deviation of survival rate, VR, as VR = i * 0.5 * CV max * VR, where CV max is the maximum coefficient of variation of a probability (Morris & Doak, 2004).
Following Koons et al. (2016), we set CV max to 1 for φ.
comparisons of effect sizes, we transformed ln(life expectancy) and iteroparity into z-scores.

| Environmental sensitive matrix population models
We then created an environmentally sensitive MPM to include the climate effect size (component 2) in our projections (Figure 1c). We kept fixed, but simulated variation in the other vital rates (S J , S A and φ) by mapping the normally distributed temporal sequence v (Equation 1) in beta-distributed (for S t ) and gamma-distributed (for φ t ) values. To do so, we calculated the quantile of each V i value given the distribution of V, and computed the value of that quantile for the beta and gamma distributions ( Figure 1c). For example, if a value of V i was the 98th percentile of its normal distribution, we drew the 98th percentile from the beta distribution for S J,t and S A,t (see S1.2). Note that while the means of these vital rates always remained the same, we scaled standard deviations by a factor V (which ranged from 0.01 to 1).

| Population projections
Using the MPMs and the environmental time series, V, we investigated the effect on populations when some of their vital rates respond to a recent climate driver, and others respond to lagged climate. Using these TVRs to climate drivers, we projected the population dynamics over 10,000 time steps (Equation 2 and 3, Figure 1d) and calculated λs using the popbio package (Stubben et al., 2016) in R (R Core Team, 2021). We obtained a temporal sequence from Equation 1 (hereafter, V recent ) and created the corresponding lagged temporal sequence by offsetting the C sequence of V recent by one step, and a new ε sequence to create V lagged (Figure 1a). In the 'control' simulations, S J , S A and φ responded to the same C sequences, but different sequences. In the 'TVR' simulations, the fecundity vital rate (φ) responded to V recent , but the survival vital rates (S J and S A ) responded to V lagged . In these simulations, all vital rates respond in the same (positive) trend to the V sequences, thus creating positive covariance between the vital rates (component 3). Initial analysis showed that there was no difference in λs when fecundity instead of the survival responds to V lagged ; therefore, we only show the first.
Next, we repeated the simulations described above, but with the assumption of negative covariance (

| Correlation with life-history traits
We investigated how mean ln(life expectancy) and degree of iteroparity correlated with our simulation results. We first calculated the log relative decrease in λ s from σ c = 0.01 to 1 for both the TVR and

| Simulations across life histories
The ( Figure 3a). Second, these effects of TVR are amplified for larger values of c and autocorrelation (Figure 3a). Regarding autocorrelation, its direct effect is minuscule when compared to its interaction with TVR (Figure 3a,b). Third, and importantly, in the case vital rates respond in opposite direction to climate, these two effects of TVR are reversed in sign, resulting in lower λ s values (Figure 3b).
The linear mixed effect model shows that the above patterns

| Correlation with life-history traits
The largest effects of TVR simulations on stochastic population growth rate (λ s ) occur for species with high life expectancy and, to a much lesser degree, species with low degree of iteroparity ( Figure 5 and S2 In this simulation, 50% of the variance was explained by a climate driver and 50% of the variance was random. In the 'control' simulation, all vital rate models respond to recent climate and in the TVR simulation, the survival vital rates (juvenile and adult survival) respond to 1-year lagged climate, whereas the fecundity vital rate responds to recent climate. Simulations were done under (i) 0.6, (ii) 0 and (iii) −0.6 environmental autocorrelation in the climate driver. (b) The density distribution of the interannual difference in λ for the whole 10,000-year sequence under (i) 0.6, (ii) 0 and (iii) −0.6 environmental autocorrelation in the climate driver.

F I G U R E 3
Responding to both lagged and recent climate (i.e. temporally varied responses-TVRs) can either buffer or amplify the negative effect of increasing environmental standard deviation (σ V ). Projected stochastic population growth rates of a life history over a range of environmental variation and climate autocorrelations, using a 2-by-2 matrix population model. In this simulation, 50% of the variance was explained by a climate driver and 50% of the variance was random. We included two types of simulations, the control where all vital rates respond to the same (recent) climate, and the TVR simulations, where the vital rates in the somatic submatrix (survival) respond to climate that is 1 year lagged from that of the reproductive submatrix (fecundity). In (a), all vital rates respond in a positive direction to the climate driver. In (b), the reproductive submatrix responds to the climate driver in the opposite direction of the survival submatrix.

F I G U R E 4
Selected coefficient estimates and 95% confidence interval of the two linear mixed effect models, relating the stochastic log lambda of population dynamic simulations to climate variables under either positive (in black) or negative (in red) correlation between the survival and fecundity vital rates. Climate autocorrelation (r K ) is the autocorrelation in the climate sequence used in the simulations, ranging from −0.6 to 0.6. Climate signal strength (p) is the relative importance of the climate sequence compared to random noise, ranging from 0.01 to 1. Temporally varied response (TVR) simulation type is the difference between the control and TVR simulations. to detrimental (negative correlations, Figure 6).

| DISCUSS ION
There is concern that increased climate variance poses a threat to populations, which has motivated considerable interest in understanding this topic (e.g. Boyce et al., 2006;Vázquez et al., 2017).
Here, we found that when vital rates of populations respond to climate with a mix of more recent and lagged climate driver timing Relative difference in the decrease in stochastic population growth rate (λ s ) between simulations with TVR and control simulations, across a range of iteroparity and life expectancy. Colours show the predicted values of linear mixed effect models relating TVR to stochastic population growth rate (λ s ). Difference in λ s was defined as the change in stochastic population growth rate from climate standard deviation of 0.1 to 1. The relative difference in λ s was calculated by dividing the difference of the control simulations by the TVR simulations. Positive values indicate that TVR is beneficial for the population growth rate (i.e. has a lower decrease in λ s compared to the control simulations), whereas negative values indicate that the responding with all vital rates to the same time window (control) is beneficial for the population growth rate. Each circle represents one of the 147 simulated life histories. In the TVR simulations, survival responds to climate that is 1 year lagged from that of fecundity; in the control simulations, all vital rates respond to the same (recent) climate. The results in the graphs refer to a climate signal strength of 0.5 (i.e. 50% of the vital rate's variance is driven by climate). Columns refer to three levels of autocorrelation (−0.6, 0 and 0.6). Rows refer to positive vital rate correlation (where all vital rates respond positively to the climate driver), or negative vital rate correlation (where the survival and fecundity vital rates respond in different directions to the climate driver).
TVR. Second, the conditions that lead to TVR buffering the effects of environmental variance are likely common in nature. Thus, our results encourage empirical studies to identify TVR, and to include them in population projection models.
Our results are perhaps the first to suggest that environmental autocorrelation might affect λs indirectly, by affecting vital rates covariation via TVRs. Previous studies found that environmental autocorrelation has relatively small direct effects on population growth rates (Eckhart et al., 2011;Paniw et al., 2018). As a result, many researchers currently investigate stochastic population dynamics under the assumption that no autocorrelation is present (e.g. Compagnoni et al., 2016;McDonald et al., 2017). Our simulations show that environmental autocorrelation can have up to 10 times larger an effect on λ s through interaction with TVR, compared to its direct effect. Temporal autocorrelation is expected to increase under climate change (Di Cecco & Gouhier, 2018), which would decrease the indirect buffering effect of TVR. However, significant regional variation in trends are also expected (Di Cecco & Gouhier, 2018), including regional decreases in autocorrelation, which would actually decrease their extinction risk.
We show that the presence of TVR reduces the annual varia- Based on physiological principles, we expect that many natural populations have TVR. For example, many plant species are known to have preformation of leaves and/or inflorescences more than 12 months before emergence (e.g. Diggle, 1997;Inouye, 1986). Thus, the vital rates associated with growth and/or fecundity will respond to climate drivers in this same timeframe as the preformation (e.g. Evers et al., 2021). In combination with F I G U R E 6 Coefficient estimates and 95% confidence interval of the two linear mixed effect models, relating the relative decrease in stochastic population growth rate to different life history and climate variables, and their interactions, under either positive (in black) or negative (in red) correlation between the survival and fecundity vital rates. ln(life expectancy) and degree of iteroparity are scaled variables for comparison. Climate autocorrelation (r K ) is the autocorrelation in the climate sequence used in the simulations, ranging from −0.6 to 0.6. Climate signal strength (p) is the relative importance of the climate sequence compared to random noise, ranging from 0.01 to 1.
possible frost damage that has rapid demographic consequences (e.g. Iler et al., 2019), alpine species could be a prime example of species with TVR. The presence of below-ground rhizomes is another physiological characteristic that has been linked to lagged climate drivers in Heliconia acuminata (Scott et al., 2021).
For this species, more immediate climate responses have been observed as well (Westerband & Horvitz, 2017 (Evers et al., 2021;Scott et al., 2021;Tenhumberg et al., 2018;van de Pol et al., 2016). Thus, we encourage researchers with long-term demographic data to explicitly test for TVR. A large literature on TVR would provide a better understanding of their prevalence and underlying mechanisms. Furthermore, our results are relevant to conservation research aimed at understanding and accurately forecasting the dynamics of populations known to be threatened by climate change (e.g. Compagnoni et al., 2021;Lindell et al., 2022). In the case of species particularly sensitive to climatic variation, explicit modelling of TVR could substantially change forecasts by correctly accounting for the indirect effects of climatic autocorrelation.
We have shown that populations that respond to a mix of temporal climate drivers can be buffered from increasing climate variance.
We have also shown that climatic temporal autocorrelation, often acknowledged but unmodelled, can either increase or dampen the effect of variability when driven by mixed temporal climate drivers.
Thus, explicitly accounting for mixed temporal climatic drivers might be an overlooked avenue to improve our understanding of population responses to future climatic change. Access funding enabled and organized by Projekt DEAL.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors have no conflict of interest to declare.

PE E R R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/api/gatew ay/wos/peer-review/10.1111/ 1365-2745.14155.

DATA AVA I L A B I L I T Y S TAT E M E N T
All code used for the analyses of this paper are archived on Zenodo https://doi.org/10.5281/zenodo.8026981 (Evers et al., 2023).