Evolutionary responses to a changing climate: Implications for reindeer population viability

Abstract If we want to understand how climate change affects long‐lived organisms, we must know how individuals allocate resources between current reproduction and survival. This trade‐off is affected by expected environmental conditions, but the extent to which density independent (DI) and density dependent (DD) processes interact in shaping individual life histories is less clear. Female reindeer (or caribou: Rangifer tarandus) are a monotocous large herbivore with a circumpolar distribution. Individuals that experience unpredictable and potentially harsh winters typically adopt risk averse strategies where they allocate more resources to building own body reserves during summer and less to reproduction. Such a strategy implies that the females do not reproduce or that they produce fewer or smaller offspring. A risk averse strategy thus results in females with large autumn body reserves, which is known to increase their survival probabilities if the coming winter is harsh. In contrast, females experiencing predictable winters may adopt a more risk prone strategy in which they allocate more resources to reproduction as they do not need as many resources to buffer potentially adverse winter conditions. This study uses a seasonal state‐dependent model showing that DD and DI processes interact to affect the evolution of reproductive strategies and population dynamics for reindeer. The model was run across a wide range of different winter climatic scenarios: One set of simulations where the average and variability of the environment was manipulated and one set where the frequency of good and poor winters increased. Both reproductive allocation and population dynamics of reindeer were affected by a combination of DI and DD processes even though they were confounded (harsh climates resulted in lowered density). Individual strategies responded, in line with a risk sensitive reproductive allocation, to climatic conditions and in a similar fashion across the two climatic manipulations.

.1. Simulated environmental conditions ( ) for 250 years for environmental conditions generated from a: (a) standard normal distribution [ 0,1 ]: (b) skew-normal distribution with a shape parameter of -2.5 [ , , 2.5 ]; and (c) skew-normal distribution with a shape parameter of 2.5. Blue lines in the right panel show the probability density for the skewnormal distributions, whereas the red lines show the probability density for the standard normal distribution. The SKEW simulations are in line with the predictions from the literature: future global climate change will most likely result in a shift towards more frequent extreme precipitation events (e.g. Wilby & Wigley, 2002, Semmler & Jacob, 2004, Tebaldi et al., 2006, Benestad, 2007, Sun et al., 2007, a trend that is already empirically evident on several continents (Sun et al., 2007 and references therein). Hanssen-Bauer et al. (2005), for example, review several studies predicting how climate will change in Fennoscandia in the future (see also Benestad, 2011): 1) increased warming rates with distance to the coast; 2) higher warming rates in winter compared to summer; and 3) increased precipitation especially during winter. The shifts between warm and cold periods during winter coupled with an year-round increased intensity of precipitation (Hanssen-Bauer et al., 2005), will lead to an increased frequency of wet weather, deep snow and ice crust formation that has negative consequences for large herbivores (e.g. Solberg et al., 2001). Rising temperatures and changing precipitation patterns have already been suggested to lead to population declines for Rangifer (Vors & Boyce, 2009). In the European Arctic and Sub-Arctic, increased occurrences of rain-on-snow and freeze-thaw events have had negative impacts on reindeer demography and population growth (Hansen et al., 2011, Solberg et al., 2001, and these are candidates for causing more frequent population collapses (see also Pape & Löffler, 2012). Yet, not all the predicted changes have to be negative (Pape & Löffler, 2012: Table 3). For semidomestic reindeer in Europe, pasture quality (earlier spring, longer growing season and higher biomass of edible plants) combined with climate change are believed to have negative effects in Scandinavia (Sweden and Norway), to be neutral in Finland and to have positive effects in Russia (Rees et al., 2008). It is thus, uncertain if the overall effect of climate change will be positive or negative, or, alternatively neutral.

NORMAL DISTRIBUTION (NORMAL simulations)
For these simulations, climatic conditions were simulated as white noise characterized by a given average ( ) and a standard deviation ( ): [ , ]. The three averages in the middle (i.e. 0.15, 0 and 0.15; see main text for details) represented the same environments as in Bårdsen et al. (2011), whereas two more extreme environments has been added in the present study (i.e. 0.30, and 0.30). The standard normally distributed environment, i.e. 0,1 , may be viewed as a reference as this resents the distribution currently observed for the AO (see Appendix A1 in Bårdsen et al., 2011 for details). For this distribution there exists three continuous climatic gradients going from (Fig. S1.2a): 1) good to poor (defined as increased ); 2) predictable to unpredictable (defined as increased ); and 3) benign to harsh (defined as the simultaneous increase in and ).

SKEW-NORMALLY DISTRIBUTION (SKEW simulations)
For these simulations, climatic conditions were simulated as a skew-normal distribution [ , , ], where , and represent the location, scale and the shape parameters (Azzalini,

Appendix S2: Model Initiation
BACKGROUND I initiated all simulations with the same population, i.e. using the exact same individuals, in order to keep initial conditions constant across the simulations. The model was initiated at by creating 200 individuals ( ), which resembles a population density of 0.40 individuals km -2 , with a constant age of 2 years: (1) Spring body mass was a random number drawn from a normal distribution with a given average ( ) of 60.7 kg and a standard deviation ( ) of 5.0 kg [i.e. , ].
(2) Similarly, the threshold spring body mass for reproduction was generated as 43.2,3.0 . (3) Then, the population was divided into two equally sized subsets ( and ). (4) and , which substitutes the range in the deterministic values for or , respectively, were defined as follows: ∈ , … , added to ; and ∈ , … , added to (the steps between the minimum and maximum was equally spaced). (5) Finally, and , which substitutes the realized values for and at , was created by adding a normally distributed error to these deterministic values: For , (defined as 0, ; and for (defined as 0, . Figure S2.1. (a) Between-individual variability in strategies, where each line represents one female, present in the population when the simulation was initiated (i.e. at ). The thick lines in blue and red shows the two most winning strategies in the previous study (Bårdsen et al., 2011). (b) Normal quantile-quantile (QQ) plots for initial spring body mass. The point and line shows the given average and standard deviation used for generating the body masses. At first glance, this might be viewed as unnecessary complex, but the point of this exercise was to ensure a high degree of between-individual variability during model-initialization with respect to the parameters defining the strategies (Fig. S1.1a) and body mass ( Fig. S1.1b). This was assured by combining different values for the three genotypic traits (the deterministic part of the initiation; point 3-4 above) while inducing some randomization to this process (point 5 above; Table S2.1 provides information about the parameter used).   (Bårdsen et al., 2011) Maximum (max) possible age 16 (year) (Bårdsen et al, 2011) Initial spring body mass was generated from a normal distribution with a given mean ( ) and standard deviation ( ) [i.e. , ] 60.70 (kg) (Bårdsen et al., 2011) 5.00 (kg)

REFERENCES
Reproductive spring body mass threshold was generated from a normal distribution with a given mean ( ) and standard deviation 43.20 (kg) (Bårdsen et al., 2011) 3.00 ( Winter climatic conditions I ran the model with two sets of winter climatic conditions ( ). I modelled winter climate on a relative scale where 'less is better' in the sense that large positive values represent poor conditions (details provided in Appendix S1).

Climatic scenario I: normally distributed environments (NORMAL)
First, similar to Bårdsen et al. (2011), was modelled as white noise characterized by a given average ( ) and standard deviation ( ): [ , ]. Five different average climatic regimes ( 0.30, 0.15,0,0.15,0.30 ) was then applied (the two most extreme values were not used in Bårdsen et al., 2011). Different levels of climatic variability was added to each average condition ( 0.04,0.08,0.12, … , 2.00 ), and this resulted in a total of 250 possible environments ( × ). Nonetheless, due to extinction before convergence in the most extreme environments, i.e. those with high values for and/or , the total number of simulations was only 147 (i.e. extinction occurred in ~40 % of the simulations). For this distribution there exists three continuous climatic gradients going from (see Appendix S1): (1) good to poor (defined as increasing ); (2) predictable to unpredictable (defined as increasing ); and (3) benign to harsh (defined as a simultaneous increase in and ).

Climatic scenario II: skew-normally distributed environments (SKEW)
was modelled as a skew-normal distribution in order to simulate an increased frequency of both poor and good environments, through varying the shape parameter (hereafter referred to as environmental skewness: 2.5, 2.46, 2.42 … ,2.5 , see Appendix S1 for details), which resulted in a total of 126 potential simulations. For this distribution there exists only one continuous climatic gradient going from beneficial (high frequency of good conditions) to detrimental [high frequency of poor conditions; defined as increased skewness values ( ; Appendix S1)]. Extinction before convergence happened in the most extreme environments (the most extreme positive with live population was 0.58). Only 72 simulations contained live populations in the end so similar to the NORMAL simulations extinction occurred in ~40 % of the runs. When 0 this resembles the standard normal distribution, i.e. 0,1 , which was the only simulation directly comparable to any of the NORMAL-scenarios. When 0 and 0 the distribution of the simulated environments were skewed towards negative and positive values, respectively. This, however, also changed the mean and the variability of the distribution (Appendix S1).

Convergence
In the previous model convergence was reached when only one strategy was left in the population. Because strategies were allowed to evolve continuously in the present model, assessing convergence became more difficult. Natural selection is a process that leads to traits becoming more or less common in a population, which means that natural selection affect the distribution of traits within a population over time (e.g. Coulson et al., 2006). Consequently, convergence was achieved when evolution had stabilized the distribution of the three genotypic traits defining the reproductive strategies. I thus assumed convergence to be achieved when the variability of all genotypic traits showed no temporal trends even though they were allowed to vary both between-and within-years. Convergence was assessed within a subjectively defined time window ( ) set to 75 years, and convergence was achieved when no significant temporal trends in any of the traits occurred within this time window. Convergence was thus assessed iteratively for each time step ( ) in a simulation. As convergence cannot be achieved when , convergence for genotypic trait (where can be substituted by , and ; see eqn. 1-2 in the main text) was assessed and updated as follows when : (1) the year-specific coefficient of variation, i.e. the estimated standard deviation divided by the estimated average, was calculated for each trait [ ]; (2) for a given time step ( ′) in a simulation, the values for from to was extracted; (3) the residuals from a linear regression model fitting only the intercept was fitted using as the response [in the software R: lm( ~1)]; (4) any significant, positive or negative, correlations (using the cor.test function in R) between the residuals and the specified time period revealed that the model had not converged with respect to genotypic trait (when convergence was achieved, i.e. when this correlation was non-significant, no further testing was performed); and (5) the simulation converged when no such significant correlation was apparent for any of the three genotypic traits.
Strictly speaking, I assessed convergence as a lack of any temporal trends in the variance of each genotypic trait, but a visual inspection of the annual averages and medians also indicates a lack of any temporal trends for these statistics. To show details about this I checked the details for a selection of environments: 1) the standard normal environment (Fig S3.1); 2) the most benign environment (Fig S3.2); and 3) the harshest environment ( Fig. S3.3).

ANALYZING AND INTERPRETING RESULTS Responses and predictors
At the point of convergence the simulations ran for 200 more years (the data collection period) reaching terminal time ( ) when the simulation was ended. The following annual population-level data was recorded: (1) across-year average density ( ); and (2) winter climate conditions; acrossyear environmental average ( ) and its standard deviation [sd ] for the NORMAL simulations and skewness ( ) for the SKEW simulations. Moreover, the average across all individuals and years of the following individual-level data was recorded across the whole data collection period: (1) the three genotypic traits ( , and ) as well as the resulting phenotypic expression ( ); (2) female reproductive success ( ; offspring per female on log e -scale); (3) spring and autumn body mass ( and ) for females and offspring, respectively; and (4) female age ( ). The utilization of values across years and individuals differ from the previous study (Bårdsen et al., 2011:249).

PSEUDO-EMPIRICAL ANALYSES Evolution of strategies
The analyses of the dynamic optimization assessed how reproductive allocation, and the three genotypic traits, was dependent on environmental conditions and density. In addition, I also analyzed more commonly used empirically used measures of individual life histories.

Allocation of resources to reproduction and survival
Generalized Additive Models (GAM), using the mgcv library (Wood, 2012), were fitted using estimated averages responses (see above) from a given simulation. The resulting dataset contained one data point for each response/predictor(s) from a given simulation. Thin plate regression splines were used to model potential non-linear effects of average population density ( ) and the interaction between both environmental characteristics [i.e. the estimated climatic variability [sd ] and average climate ( )] (Wood, 2006). The following code was applied in the R software: 1) the Density Independent (DI) model specified as 'gam(response ~ s(sd , , bs=tp, k=4)' and 'gam(response ~ s( , bs=tp, k=4)' for the NORMAL and SKEW simulation, respectively; and 2) the Density Dependent (DD) model specified as 'gam(response ~ s( , bs=tp, k=4)'. Plotting of results from the DI model in the NORMAL simulations was performed using the vis.gam function (as this model contained two continuous predictors), whereas plotting the DI and DD model in the remaining analyses were performed using the predict.gam function (Wood, 2012).

Population dynamics: time series analyses using autoregressive models
As in the previous study, I assessed to what extent climatically induced changes in life histories affect population dynamics, which were analyzed using time series analysis. The population growth rate, i.e.
, where refers to populating abundance in a given year ( ), was modelled using second-order autoregressive models [AR(2), fitting an ARIMA(p = 2, d = 0, q = 0) model using the arima-function where was included as a covariate in the xreg argument]. The linear predictor of the models included direct ( ) and delayed ( 1) density dependence (regulation) and the direct effect of on . I, thus, estimated the first-order AR coefficient (1 ), the second-order AR coefficient ( ) and the direct effect of winter climate conditions ( ). Statistical analyses, using the resulting coefficients from the time series analyses as responses, and plotting was performed similarly as in the previous paragraph.

BACKGROUND
In this appendix I show more detailed results compared to what is shown in the main text.

Normally distributed environments (NORMAL)
Evolution of strategies Reproductive allocation ( � ), which is the strategies' phenotypic expression, was positively related to both environmental average ( � ) and environmental variability [sd( ) � ], and negative related to population density ( � ; Table S4.1a & Fig. S4.1). This lead to the surprising result that female reindeer allocated more resources to reproduction when environmental conditions were harsh (i.e. in poor and unpredictable) environments. The range in the predicted values of � was, however, small while at the same time as the deviance explained was >20% higher for the density dependent (DD) model compared to the density independent (DI) model ( Fig. S4.1). Additionally, population density was highest in benign (i.e. in good and predictable) environments and lowest in harsh environments ( Fig. S4.2). The fact that climate was such a powerful predictor of density ( Fig. S4.2) implies that density and climate was confounded -and hence their effects may be difficult, if not impossible, to disentangle in simulations where harvest or predation does not occur (see main text for details). The relative difference in the deviance explained for the DD model, however, indicates that density might be a more important predictor of � than climate, which showed the pattern that I a priori expected. The intercept (�) and the body mass threshold (�; the lowest spring body mass the females' needed for engaging in reproduction), i.e. two of the strategies' genotypic traits, showed, similar to � , positive relationships with � and sd( ) � , and both were subject to negative density dependence (Table S4.2b,d & Fig. S4.1b,d). Judging from the percentage deviance explained, � was more affected by climate than density, whereas this was the opposite for body mass threshold. The body mass threshold was, however, poorly explained by both the DD and DI models. The slope for spring body mass ( � ), which was the only genotypic trait being optimized in the previous study, showed negative relationships with both � and sd( ) � , and a positive relationship with density (Table S4.1c & Fig. S4.1c).

Reproductive allocation: reproductive success and offspring body mass
By using measures of reproduction commonly used in empirically studies, I attempted to relate model output to the empirical world. Reproductive success ( � : the number of offspring per female on loge-scale) was highest in benign and lowest in harsh environments [i.e. negatively related to both � and sd( ) � ]. This was expected, but it was the opposite what I found for � . Females also produced more offspring in high-vs. low-density environments (Fig. S4.3a & Table S4.2a), but such a positive density dependence was unexpected. Again, this could be because density and climate was confounded: for � , the DI model the deviance explained was more than three-fold compared to the DD model. This might indicate that the relationship between � and climate, which supported my a priori expectations, was stronger than the relationship between � and density ( Fig. S4.3a). It is, however, important to note that the range in the predicted values for � was low. Offspring spring and autumn body mass was highest in harsh and lowest in benign environments, but showed evidence of negative density dependence (Fig. S4.3b,c & Table S4.2b,c).

Somatic allocation: female age and body mass
Female age was highest in benign and lowest in harsh environments, whereas females became older as density increased (Fig. S4.4a & Table S4.3a). In the autumn, females were largest in harsh environments, and showed clear evidence of negative density dependence (Fig. S4.4b & Table  S4.3b). In spring, however, females were smallest in poor and predictable environments and largest in good and unpredictable environments, but showed no evidence for density dependence (Fig.  S4.4c & Table S4.3c).

Population dynamics: time series analyses
The effect of direct regulation (1 − 1 � ) was most pronounced, i.e. the most negative, in harsh environments and negatively related to density (Fig. S4.5a & Table S4.4a). This means that direct negative density dependence increased, i.e. became more negative, along the benign-harsh gradient and as density increased. Delayed regulation ( 2 � ) was not significantly related to climate, but were negatively related to density (even though the explanatory power of both these model were poor: Fig. S4.5b & Table S4.4b). The direct effect of climate (ω 1 � ) on population growth was highest in poor and relatively stable climatic conditions, and unaffected by density (Fig. S4.5c & Table S4.4c).

Skew-normally distributed environments (SKEW)
Evolution of strategies Unexpectedly, reproductive allocation ( � ) was positively related to the environmental skewness ( ), which means that females allocated more to reproduction when the frequency of poor winter conditions increased (i.e. under detrimental conditions: Appendix S1), and negatively related to density ( � ): Table S4.5a & Fig. S4.6a]. The fact that � increased along the beneficial-detrimental gradient could, similar to the findings for the normally distributed environment, be due to a confounding between climate and density (Fig. 3b in the main text). Again, the DD model, which supported my a priori expectations, had a higher explanatory power compared to the DI model ( Fig. S4.6a), which gave unexpected results. Both the intercept (�) and the body mass threshold (�) were positively related to �, and both were subject to negative density dependence (Table  S4.5b,d & Fig. S4.6b,d). The slope for spring body mass ( � ) showed negative relationships with environmental skewness, and a positive relationship with density (Table S4.5c & Fig. S4.6c). Judging from the percentage deviance explained all four responses were more affected by density compared to climatic conditions.

Reproductive allocation: reproductive success and offspring body mass
Even though the positive relationship between � and � was unexpected, � was negatively related to � but positively related to density (Fig. S4.7a & Table S4.6a). Offspring spring and autumn body mass was, however, positively related to �, but showed evidence of negative density dependence ( Fig. S4.7b,c & Table S4.6b,c). As in the normally distributed environments climate seemed to be a more important predictor of these responses than population density (judging from the deviance explained by the DD and DI models, respectively). This means that the relative importance of the DD and DI model was the opposite and in favour of climate for all the 'empirical' measures of reproductive allocation than it was for the 'theoretical' measures presented above.

Somatic allocation: female age body mass
Female age was negatively related to �, indicating that age was highest when the frequency of good winters was high, whereas female age increased as density increased (Fig. S4.8a & Table S4.7a). In the autumn, female body mass was positively related to �, and showed clear evidence of negative density dependence (Fig. S4.8b & Table S4.7b). Neither � nor density explained spring body mass (Table S4.7c).

Population dynamics: time series analyses
The effect of direct regulation (1 − 1 � ) showed a curved (concave-down) relationship with �, being lowest for intermediate values of environmental skewness, and a linear and negative relationship with density ( Fig. S4.9a & Table S4.8a). Delayed regulation ( 2 � ) was not statistical significantly related to either climate or density (Table S4.8b). The direct effect of climate (ω 1 � ) on population growth was highest, i.e. most negative, at high frequency of poor environments, and the estimated effect of climate in the AR-models increased as density increased (Fig. S4.9c & Table S4.8c). Once again, density was to a high degree explained by climatic conditions (Fig. S4.10).  Table S4.1 provides detailed GAM output.  Table S4.2 provides detailed GAM output.  Table S4.3 provides detailed GAM output.  ( 2 ) and (c) direct effects of climate ( 1 ), was a function of the interaction between environmental variability and average (left panel) and population density (right panel). Table S4.5 provides detailed GAM output.  Table S4.5 provides detailed GAM output. Please note that zero values for � represents the standard normal distribution, and represents a baseline for comparison, whereas a negative � gives a skew towards the left (negative values for the outcome) and a positive � gives a skew towards the right.  Table S4.7 provides detailed GAM output.  Table S4.8 provides detailed GAM output.  Table S4.8 provides detailed GAM output.  i.e. the interaction between environmental average, � , and standard deviation, sd( ) � ] and population density (ii). Estimated degrees of freedom (edf) provide an estimate of the degree of complexity in the relationship, whereas deviance explained (D) provides an estimate of how well the model explains the variance in the response (Fig. S4.1 provides a   i.e. the interaction between environmental average, � , and standard deviation, sd( ) � ] and population density (ii). Estimated degrees of freedom (edf) provide an estimate of the degree of complexity in the relationship, whereas deviance explained (D) provides an estimate of how well the model explains the variance in the response (Fig. S4.3 ( 2 ) and (c) direct effects of climate ( 1 ), was a function of environmental skewness [i.e. the shape parameter ( �) in the skew normal distribution; left panel] and population density (ii). Estimated degrees of freedom (edf) provide an estimate of the degree of complexity in the relationship, whereas deviance explained (D) provides an estimate of how well the model explains the variance in the response (Fig. S4.9 provides a graphical view of these results).