The important choice of reference environment in microevolutionary climate response predictions.

It is well documented that individuals of wild populations can adjust to climate change by means of phenotypic plasticity, but few reports on adaptation by means of genetically based microevolution caused by selection. Disentanglement of these separate effects requires that the reference environment (the environmental zero point) is defined, and this should not be done arbitrarily. The problem is that an error in the reference environment may lead to large errors in predicted microevolution. Together with parameter values and initial mean trait values, the reference environment can be estimated from environmental, phenotypic and fitness data. A prediction error method for this purpose is described, with the feasibility shown by simulations. As shown in a toy example, an estimated reference environment may have large errors, especially for small populations. This may still be a better choice than use of an initial environmental value in a recorded time series, or the mean value, which is often used. Another alternative may be to use the mean value of a past and stationary stochastic environment, which the population is judged to have been fully adapted to, in the sense that the expected geometric mean fitness was at a global maximum. Exceptions are cases with constant phenotypic plasticity, where the microevolutionary changes per generation follow directly from phenotypic and environmental data, independent of the chosen reference environment.

INTRODUCTION

Wild populations respond to changing environments by means of phenotypic plasticity and microevolution, and especially climate change responses have been extensively studied. The aim is then to disentangle phenotypic changes owing to genetically based microevolution caused by natural selection, and changes due to individual phenotypic plasticity. Relying on 11 review articles, including reviews of altogether 66 field studies, Merilä and Hendry (2014) arrived at the conclusion that evidence for genetic adaptation to climate change has been found in some systems, but that such evidence is relatively scarce. They also concluded that more studies were needed, and that these must employ better inferential methods. The aim of the present article is to give a contribution in that last respect.

It is obvious that for all evolutionary systems with interval‐scaled environmental variables , as, for example, temperature in °C, a suitable zero point (reference environment) must be chosen, and as argued in Section 2, this should not be done arbitrarily. A zero point is in general defined as “the point on a scale that denotes zero and from which positive and negative readings can be made” (Collins English Dictionary). The problem is that an error in the reference environment may lead to large errors in predicted microevolution. In most cases where the environmental variable is, for example, a temperature, the reference environment should not, for example, be set to 0°C (or 0°F). Neither should it, without further consideration, be set to the initial or mean environmental value of a specific time series. It appears that the need for a proper reference environment definition, and thus also an environmental cue definition, has been largely ignored in the reviewed studies referred to in Merilä and Hendry (2014).

The present article is an attempt to clarify some important questions relating to reference environments, and for that purpose a method for model‐based predictions of microevolutionary changes is also proposed. This method is based on parameter estimation by means of prediction error minimization, including estimation of reference environment and initial mean values of quantitative traits.

For a discussion on the general microevolution versus plasticity disentanglement problem, we may for simplicity assume the intercept‐slope individual reaction norm model: where and are the environmental cue and the individual phenotypic value, respectively, as functions of time measured in generations. Here, and are the additive genetic components of the intercept and plasticity slope, respectively, while and are independent iid zero mean normal non‐additive effects. As done in Lande (2009) and Ergon and Ergon (2017), we may consider the individual reaction norm intercept , and the individual plasticity slope , as two quantitative traits in their own right. Microevolution thus involves changes in the mean trait values and from generation to generation. The generations are here assumed to be non‐overlapping.

From Equation (1) follows the mean trait reaction norm model, , and from this simple equation follows the basic questions discussed in this article. How can be estimated, and how can the evolution of and be predicted, provided that and are known? And how will the predictions be affected by errors in the estimated or assumed value ? It turns out that in order to answer these questions we also need information on individual fitness values .

The reference environment is determined by the environment at which the phenotypic variance has its minimum, as defined in more detail in Section 2, and as discussed in Ergon and Ergon (2017) and Ergon (2018). In theoretical work, it is often assumed that the population has fully adapted to a stationary stochastic environment with a given mean value, such that the expected geometric mean fitness is at a global maximum, and the reference environment is then set to zero (Chevin & Lande, 2015; Lande, 2009). Although there is nothing wrong with this theoretical approach, it disguises the underlying problem discussed here, and is therefore included in Equation (1). This formulation also makes it possible to distinguish between the environment as such and the environmental cue. In some cases, it may be possible to determine the reference environment experimentally, see, for example, Fossen et al. (2018), but that may obviously be difficult for wild populations.

When the environmental cue changes over time, the mean trait values and as follow from Equation (1) may evolve due to selection, and as a result also the mean phenotypic value will evolve (Lande, 2009). Without changes due to selection, that is, if the mean trait values and are constant, the value of may still change when changes, as also follows from Equation (1).

Section 2 discusses several aspects of the general microevolution versus plasticity disentanglement problem. First, a definition of the reference environment is given. Second, it is shown how the mean trait values and , and thus also , evolve as functions of the environmental cue and the phenotypic selection gradient . Third, it is shown how and , as well as initial mean trait values and the parameter values in the matrix, can be estimated by means of a prediction error minimization method (Ljung, 2002), using data from known time series of and , as well as of individual fitness values . Fourth, it is discussed why it may be difficult to estimate , as revealed by simulations, and which consequences errors in estimated values of will have. Exceptions are cases with constant phenotypic plasticity, where the microevolutionary changes per generation follow directly from phenotypic and environmental data, independent of the chosen reference environment.

It must be underlined that the theory in Section 2 assumes that the phenotypic trait in Equation (1) is not correlated with other phenotypic traits having causal effects on fitness, see Morrissey et al. (2010) for a discussion. Also note that the need for a proper reference environment is not specific for the simple case according to Equation (1).

Simulations in Section 3 make use of a toy example, utilizing the intercept‐slope reaction norm model in Equation (1). The environmental input is here a noisy positive trend in spring temperature, while the individual phenotypic values are the clutch initiation dates for a certain bird species. The toy example also assumes that the individual (mid‐parent) fitness values are the numbers of offspring. The essential questions are how microevolutionary changes in mean intercept and plasticity slope can be predicted, and how these predictions are affected by errors in the reference environment in Equation (1). The simulations show that errors in the estimated or assumed value of may cause large mean trait prediction errors. They also show the feasibility of the proposed parameter estimation method.

Finally, follows a discussion in Section 4. Derivations of prediction equations, simulation results with modeling error and increased population size, and a short comparison with BLUP/REML parameter estimation are given in Appendices S1–S4.

THEORY AND METHODS

Example system

For a study of the general reference environment problem, and for a test of the proposed parameter estimation method, we may consider a true evolutionary system based on Equation (1), with the additive genetic covariance matrix  , and the phenotypic covariance matrix  . Here, Equation (2b) is the multivariate breeder's equation (Lande, 1979), where is found from any given fitness function. It is assumed that the phenotypic trait in Equation (1) is not correlated with other phenotypic traits having causal effects on fitness, and that generations are non‐overlapping.

Reference environment and environmental cue definitions

As discussed in the Introduction, there is a need for reference environment and environmental cue definitions:

where is the environment at which the phenotypic variance is at a minimum, and where the covariance between the plastic phenotypic value and reaction norm slope is equal to 0. Here, is a correction term that may be 0.

Assuming a single environmental variable , and given a reaction norm model, the reference environment is as follows:

With according to Definition 1, the environmental cue is given by .

For the reaction norm model (1) we find, for example, (using ) which by setting and gives the reference environment

For the reference environment is thus the environment where the phenotypic variance is minimized (see Figure 1 for illustration). This is also the environment where the expected geometric mean fitness has a global maximum, and thus the environment the population is fully adapted to. In this environment the environmental cue will be 0.

FIGURE 1

Reaction norms for 100 individuals in a population according to Equation (1), with , and . The reference environment is °C, which since also is the temperature to which the population is fully adapted. The mean trait values are and . Solid lines indicate the range of data used for parameter estimation and mean trait predictions in simulations. Note that °C is not within that range

Mean trait prediction equations

A fundamental equation for mean trait predictions follows from Equation (2a) as where , , and are changes per generation. From this follows that the value of has nothing to say in special cases with constant phenotypic plasticity slopes, that is, when . In such cases, we simply have , where is constant, or only if does not evolve.

As shown in Appendix S1, Equation (5) leads to equations for and as functions of the phenotypic selection gradient , and where

In addition to time series of and , we thus need parameter values for , , , , , and , and a time series of individual fitness values . For mean trait predictions, we also need initial values. Note that these equations are valid only when the genetic relationship matrix is a unity matrix (Ch. 26, Lynch & Walsh, 1998).

Prediction error minimization method

From predicted mean intercept and plasticity slope values found by Equations (6a, 6b) follow predicted values of from Equation (2a). The prediction equations can thus be used for parameter estimation in a prediction error minimization method (PEM), as shown in Figure 2. As follows from Equations ((6a), (6b), (6c)), we can then set to any value, and estimate , , , and relative to that value.

FIGURE 2

Block diagram of microevolutionary PEM, with dynamical tuning model based on an intercept‐slope reaction norm model with mean traits and . Here, and are the known environmental input and the known mean phenotypic value at time , respectively. is the additive genetic relationship matrix, which here is assumed to be , while and are vectors of individual phenotypic and relative fitness values, respectively. The and matrices include the system parameters, while , , and are the initial mean trait values and the reference environment, respectively. Assuming data over generations, all these model parameters are tuned until is minimized, with and

Effects of errors in the reference environment

With a reference environment instead of , predictions based on Equation (2a) can be written where and are found from Equations (6a, 6b) with use of estimated parameter values, assuming initial values known.

For small values of , that is, when and , it follows from Equations ((6a), (6b), (6c)), that is independent of , and that is constant. This results in , such that only must be tuned in order to minimize . In this case, an error in has very little effect on the change in per generation, as also follows from Equation (5).

For larger values of , the predicted change per generation will be affected by an error in , and with good predictions for to can then only be obtained by parameter tuning such that over all generations. That is possible because appears in both nominator and denominator of Equation (6b). According to Equation (7) we then find , which as shown in Section 3 may result in large errors in predicted changes in over time.

Equation (7) also shows why it may be difficult to estimate from data based on small populations. The reason is that minimization of results in and over all generations, also if there is a large error in .

Effects of modeling errors

Modeling errors will obviously affect predictions of the mean traits. As an example, simulations with the true individual model are included in Appendix S2. Here, is a perception trait, as discussed in Ergon and Ergon (2017).

SIMULATION RESULTS

Description of toy example

In the toy example, the environmental input () is a noisy positive trend in spring temperature, resulting in a noisy negative trend in mean clutch initiation date () for a certain bird species, approximately as in figure 2 in Bowers et al. (2016). The individual phenotypic values are discrete, with days as unit. The individual (mid‐parent) fitness values () are integers from to , with number of fledglings as unit. Generations are assumed to be non‐overlapping, and the population size is assumed to be constant. Data for , and are generated over generations, where the positive temperature trend begins at generation . The population is assumed to be fully adapted to the mean spring temperature 10°C before generation 10, which is thus the reference environment , but only data from generations to are used for parameter estimation and mean trait predictions. Note that 10°C may not be within the range of input data used for parameter estimation (depending on realization). The essential questions are how may be estimated, how well microevolutionary changes in mean intercept and plasticity slope over generations 31 to 60 can be predicted by means of the PEM method in Figure 2, and how errors in the estimated or assumed value of will affect the predictions.

True model, fitness function, and environmental input signals

Assume that what we consider to be true mean responses, , and , are generated by the state‐space model (2a, 2b). Here, in the true system but left as a free parameter in the tuning model in Figure 2. The individual effects , , and are at each generation drawn from populations with normal distributions around , and , respectively.

The individual fitness function is assumed to be rounded values of where is the phenotypic value that maximizes fitness, while . The discrete values of (number of fledglings) are thus integers from 0 to 10.

Also assume a stationary or slowly varying mean of a stochastic environment, with added iid zero mean normal random variations with variance , that is, , and that the population is fully adapted to a stationary stochastic environment with (as in Figure 1). In a corresponding way, assume that , where is iid zero mean normal with variance , and where and are correlated with covariance . Following Lande (2009), we may assume that juvenile birds of generation are exposed to the environment during a critical period of development a fraction of a generation before the adult phenotype is expressed and subjected to natural selection. We will define , which implies a linear relationship , variances , and covariance , where is the autocorrelation of background environmental fluctuations. We will assume and . The optimal value of the mean plasticity slope in a stationary stochastic environment is then (as in Figure 1) (Ergon & Ergon, 2017).

Further assume that and are noisy ramp functions as shown in Figure 3, with the ramp in starting from 10°C at generations. The choice of a negative trend in , and thus in , results in earlier clutch initiation dates as a result of the positive temperature trend.

FIGURE 3

Noisy ramp functions (panel a) and (panel b), with , , , , and The solid parts of the curves indicate data used for parameter estimation and mean trait predictions (compare with Figure 1, where reaction norms for are indicated by solid lines)

Figure 4 shows typical individual phenotypic (clutch initiation date) and fitness (number of fledglings) values for the true model with population size at generation 45 in Figure 3. The figure shows that the most negative (earliest) dates give the highest number of offspring, and the population is thus under directional selection toward earlier clutch initiation dates. The zero‐point date is the mean clutch initiation date before the positive temperature trend sets in at generation 10 in Figure 3, when the population is assumed to be under stabilizing selection and fully adapted to the stationary stochastic temperature.

FIGURE 4

Individual clutch initiation dates , with a range from to  days (phenotypic values, Panel a), and number of fledglings (fitness, Panel b), at a generation where the population is under directional selection toward earlier clutch initiation dates. As indicated by dots, the number of offspring are integers, while the phenotypic values are integers divided by 10 (i.e., days)

Parameter estimation and mean trait prediction results

Parameter estimation and mean trait prediction results were found by use of the MATLAB function fmincon in the PEM method in Figure 2. Results with use of input–output data from to with population size are given in Table 1. The relative errors in total change in predictions over  generations are included, computed as etc., where and . The final values of are also included, as they indicate the degree of optimization success. Results are presented as mean values and standard errors, Mean ± SE, based on 100 repeated simulations with different realizations of random inputs.

TABLE 1

Estimation and prediction results with true system responses generated by means of Equations (2a, 2b) and (9). Results are for cases with population size and perfect observations of and , and they are based on 100 simulations with different realizations of all random input variables

Parameter etc.	True value	Results Case 1	Results Case 2	Results Case 3
G^_bb	0.01	0.0103 ± 0.0023	0.0099 ± 0.0054	0.0065 ± 0.0023
G^_ab	0	0.0033 ± 0.0141	0.0028 ± 0.0152	0.0067 ± 0.0022
σ^v2	0.025	0.0341 ± 0.0141	0.0301 ± 0.0348	0.0264 ± 0.0040
σ^η2	0.01	0.0118 ± 0.0091	0.0123 ± 0.0115	0.0157 ± 0.0037
b¯^_31	–	−0.4965 ± 0.0086	−0.4969 ± 0.0083	−0.4985 ± 0.0080
u^_ref	10	10	9.9210 ± 0.6553	11
∑εt,final2	–	10−5 (23 ± 10)	10−5 (23 ± 10)	10−5 (29 ± 11)
Δ30errora¯^_t%	–	1 ± 4	−5 ± 42	68 ± 8
Δ30errorb¯^_t%	–	−1 ± 3	0 ± 3	−4 ± 4
Δ30,correrrora¯^_t%	–	1 ± 4	0 ± 5	−4 ± 5

Case 1: (the true value). Case 2: as a free variable. Case 3: (expected initial value in optimization data). Here, 9% of the simulations were discarded because .

Given the model in Equations (2a, 2b) and (9), there are in all six parameter values to be estimated (while follows from Equation (2a) with set to 0). In the optimizations, the initial values of , , , and were set to 0, while the initial value of was set to (when was a free variable). The true value was used, such that estimates of , , and are found relative to . Table 1 presents results for three cases, first for (Case 1), second for as free variable (Case 2), and third for (Case 3), which is the expected initial value in the time series used. The estimates of and have in all cases rather large standard errors, and in some cases also large bias errors. What is more interesting is that the prediction errors and are small in Case 1 (with the correct value . In Case 2 (with as a free variable), has a large standard error. With (Case 3), has a large bias error. In this case also the estimates of , and are clearly biased. In Case 2 and Case 3, the results are close to the results with the correct value (Case 1), as explained in Subsection 2.5.

Table 1 includes theoretical prediction error results based on Equation (7), . These are in all cases close to the results with (Case 1).

The results for Case 1 and Case 2 are very much improved with population size , while the standard errors in Case 3 were only marginally improved by an increased population size (Appendix S3). The results for in Case 2 were, for example, improved from to .

As shown in Table 1, a large error in the assumed reference environment results in large errors in predicted changes in over 30 generations (Case 3). Table 2 shows prediction results for more moderate errors in , as well as for (expected mean value in optimization data), and they are all in accordance with Equation (7).

TABLE 2

Errors in predicted total relative change in and over 30 generations, as functions of the reference environment used in the optimization procedure

u^_ref	Δ30errora¯^_t%	Δ30errorb¯^_t%
9.75	−16 ± 4	0 ± 3
10	1 ± 4	−1 ± 3
10.25	17 ± 4	0 ± 3
10.5	34 ± 6	−2 ± 3
11	68 ± 8	−4 ± 4
11.75	115 ± 14	−8 ± 6

Figure 5 shows predicted mean values , and , as compared to true mean values , and , for Case 1 and Case 3 in Table 1. Note that is set to 0, and the zero point is thus not the same as in Figure 4, Panel a.

FIGURE 5

Typical responses for Case 1 and Case 3 in Table 1. True values are shown by solid blue lines. All parameter values except and were initially set to 0, which gave predictions as shown by dashed blue lines. Final predictions are shown by blue dots. True and responses are shown by green lines, while predictions and are shown by magenta circles. Panels a and b show results for Case 1 with (true value), while panels c and d show results for Case 3 with . Here, the theoretical predictions are included as black dashed line. Note that is set to 0, such that and , where is not quite the same from realization to realization

DISCUSSION

It is well documented that populations adjust to climate change by means of individual plasticity, but few reports on adaptation by means of genetically based microevolution caused by phenotypic selection (Merilä & Hendry, 2014). The main point in this article is that disentanglement of these separate effects requires that the reference environment is defined, and that this should not be done arbitrarily. Instead, it should be based on the environment where the phenotypic variance is at a minimum (Definition 1 and Figure 1). This definition can be extended to multivariate cases.

As shown in a toy example, large errors in the estimated or assumed value may lead to large errors in predicted microevolutionary changes over time (Table 1 and Figure 5). Such large errors in may occur when the range of environmental data used for predictions is far from the mean value of the stochastic environment the population is adapted to.

In the toy example, the mean plasticity slope is predicted quite well also when there is a large error in , and one reason for this is that an error in to some extent compensates for the error in (Table 1, Case 3). This does not, however, prevent a large error in the predicted mean intercept (Equation (7)).

Although the plastic response to climate change is a result of individual plasticity, it should be noted that individual traits do not determine the environmental value in Definition 1, but instead the phenotypic variance (Figure 1). Similarly, individual traits do not enter into the prediction Equations (6a, 6b), but instead the individual phenotypic and fitness values (Figure 4).

In theoretical studies, it is often assumed that the environmental variable is scaled such that (Chevin & Lande, 2015; Lande, 2009). This can be done also in databased applications, provided that is known, and that the correction term in Definition 1 is 0.

The toy example used in the simulations is a simplification of reality because changing spring temperatures affect fitness in a complex way (Bowers et al., 2016). It still suffices to show that the reference environment together with initial mean trait and parameter values can be estimated from environmental, phenotypic, and fitness data, by use of the prediction error minimization method in Figure 2. The simulations make use of an environmental trend, as in noisy temperature trends caused by climate change (Figure 3), and a correct reference environment then results in quite good predictions of changes in mean traits over time (Table 1, Case 1). Although these predictions are based on parameter estimation, all the separate parameter estimates as such are not especially good, although they were considerably improved when the population size was increased from to (Appendix S3). The reference environment can also be estimated, but with large standard errors, especially for small population sizes, and this results in a correspondingly large standard error in predicted change in mean intercept over time (Table 1, Case 2). An estimated reference environment may still be a better choice than use of an initial environmental value in a recorded time series, or the mean value, which may give large errors in predicted changes in mean traits (Table 1, Case 3). Another alternative may be to use the mean value of a past stationary stochastic environment, which the population is judged to have been fully adapted to.

It is here assumed that the genetic relationship matrix is an identity matrix, and the simulation results are obtained by use of a very simple system and a prediction error minimization method. However, the fact that errors in the reference environment may cause large errors in predictions of microevolution, as discussed in Subsection 2.5, is a generic problem. Independent of prediction method and the complexity of the model, an error in the reference environment implies that an erroneous model is fitted to the input–output data, and that must inevitably result in prediction errors. An alternative view, when is constant (Case 3), is that the tuning model in Figure 2 still uses the correct value of , but then with an error term added to the environmental input. In order to minimize , this input error must as good as possible be compensated by errors in estimated parameter values, resulting in prediction errors. Note that this argument is independent of the specific parametrizations used in the microevolutionary system and tuning model in Figure 2. It is in any case no reason to believe that prediction errors caused by will disappear in cases where the genetic relationship matrix is not a unity matrix, and when other parameter estimation and mean trait prediction methods are used. A more specific argument regarding BLUP/REML parameter estimation is given in Appendix S4. It must thus be expected that predictions of microevolutionary changes over time depend on the chosen reference environment, and such predictions cannot therefore be trusted unless the chosen reference environment can be trusted. Exceptions are here cases with a constant mean plasticity slope, where the change in mean reaction norm intercept per generation according to Equation (5) is independent of the reference environment. This implies that a nearly constant mean plasticity slope must be expected to result in small errors in the predicted changes in the mean intercept, also if there is an error in the reference environment.

CONFLICT OF INTEREST

The author declares no conflict of interest.

AUTHOR CONTRIBUTION

Rolf Ergon: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal).

Supplementary Material

Click here for additional data file.