The important choice of reference environment in microevolutionary climate response predictions

Abstract 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.

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 u t , as, for example, temperature in °C, a suitable zero point (reference environment) u ref 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 u t − u ref and y i,t are the environmental cue and the individual phenotypic value, respectively, as functions of time t measured in generations. Here, a i,t and b i,t are the additive genetic components of the intercept and plasticity slope, respectively, while v i,t and i,t 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 a i,t + v i,t , and the individual plasticity slope b i,t + i,t , as two quantitative traits in their own right. Microevolution thus involves changes in the mean trait values a t and b t from generation to generation.
The generations are here assumed to be non-overlapping.
From Equation (1) follows the mean trait reaction norm model, y t = a t + b t u t − u ref , and from this simple equation follows the basic questions discussed in this article. How can u ref be estimated, and how can the evolution of a t and b t be predicted, provided that u t and y i,t are known? And how will the predictions be affected by errors in the estimated or assumed value û ref ? It turns out that in order to answer these questions we also need information on individual fitness values W i,t .
The reference environment u ref 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 u ref 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 u t − u ref changes over time, the mean trait values a t and b t as follow from Equation (1) may evolve due to selection, and as a result also the mean phenotypic value y t will evolve (Lande, 2009). Without changes due to selection, that is, if the mean trait values a t and b t are constant, the value of y t may still change when u t − u ref 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 a t and b t , and thus also y t , evolve as functions of the environmental cue u t − u ref and the phenotypic selection gradient y,t . Third, it is shown how u ref and y,t , as well as initial mean trait values and the parameter values in the G matrix, can be estimated by means of a prediction error minimization method (Ljung, 2002), using data from known time series of u t and y i,t , as well as of individ- It must be underlined that the theory in Section 2 assumes that the phenotypic trait y i,t in Equation (1)  (1) 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.

| 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), Here, Equation (2b) is the multivariate breeder's equation (Lande, 1979), where W i,t is found from any given fitness function. It is assumed that the phenotypic trait y i,t 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: Definition 1 Assuming a single environmental variable u t , and given a reaction norm model, the reference environment is as follows: where u 0 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, f(reaction norm parameter values) is a correction term that may be 0.
Definition 2 With u ref according to Definition 1, the environmental cue is given by u − u ref .
For the reaction norm model (1) we find, for example, (using which by setting cov(y, b) = 0 and u � = u 0 − u ref gives the reference environment For G ab = 0, 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.

| Mean trait prediction equations
A fundamental equation for mean trait predictions follows from Δy t = y t+1 − y t are changes per generation. From this follows that the value of u ref has nothing to say in special cases with constant phenotypic plasticity slopes, that is, when Δb t = 0. In such cases, we simply F I G U R E 1 Reaction norms for 100 individuals in a population according to Equation (1), with G aa + 2 v = 0.05, G bb + 2 = 0.02, and G ab = 0. The reference environment is u ref = 10°C, which since G ab = 0 also is the temperature u 0 to which the population is fully adapted. The mean trait values are a t = 0 and b t = − 0.5 . Solid lines indicate the range of data used for parameter estimation and mean trait predictions in simulations. Note that u ref = u 0 = 10°C is not within that range As shown in Appendix S1, Equation (5)

| Prediction error minimization method
From predicted mean intercept and plasticity slope values found by Equations (6a, 6b) follow predicted values of y t 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-6c), we can then set G aa to any value, and estimate G ab , G bb , 2 v , and 2 relative to that value. For small values of G bb , that is, when G bb → 0 and G ab → 0, it fol-

| Effects of errors in the reference environment
. In this case, an error in û ref has very little effect on the change in â t per generation, as also follows from Equation (5).
For larger values of G bb , the predicted change per generation Δâ t will be affected by an error in û ref , and with G ab = 0 good predictions ŷ t ≈ y t for t = 1 to T can then only be obtained by parameter tuning such that b t ≈ b t over all generations. That is possible because u ref appears in both nominator and denominator of Equation (6b).
According to Equation (7) we then find â t ≈ a t +b t (û ref − u ref ), which as shown in Section 3 may result in large errors in predicted changes in a t over time.
Equation (7) also shows why it may be difficult to estimate u ref from data based on small populations. The reason is that minimiza- over all generations, also if there is a large error in û ref .

| 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, c i,t + i,t is a perception trait, as discussed in Ergon and Ergon (2017).

| Description of toy example
In the toy example, the environmental input (u t ) is a noisy positive trend in spring temperature, resulting in a noisy negative trend in mean clutch initiation date (y t ) 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 (W i,t ) are integers from 0 to 10, with (6a) F I G U R E 2 Block diagram of microevolutionary PEM, with dynamical tuning model based on an intercept-slope reaction norm model with mean traits a t and b t . Here, u t and y t are the known environmental input and the known mean phenotypic value at time t , respectively. A t is the additive genetic relationship matrix, which here is assumed to be A t = I, while y t and w t are vectors of individual phenotypic and relative fitness values, respectively. The Ĝ and P matrices include the system parameters, while â init , b init , and û ref are the initial mean trait values and the reference environment, respectively. Assuming data over T generations, all these model parameters are tuned until ∑ T with y t =ŷ 1 = 0 and â 1 = −b 1 u 1 −û ref number of fledglings as unit. Generations are assumed to be nonoverlapping, and the population size is assumed to be constant.
Data for u t , y i,t , and W i,t are generated over 60 generations, where the positive temperature trend begins at generation 10. The population is assumed to be fully adapted to the mean spring temperature 10°C before generation 10, which is thus the reference environ-

| True model, fitness function, and environmental input signals
Assume that what we consider to be true mean responses, y t , a t , and b t , are generated by the state-space model (2a, 2b). Here, G ab = 0 in the true system but left as a free parameter in the tuning model in Also assume a stationary or slowly varying mean U,t of a stochastic environment, with added iid zero mean normal random variations u n,t with variance 2 U n , that is, u t = U,t + u n,t , and that the population is fully adapted to a stationary stochastic environment with U,t = u ref = u 0 = 10 • C (as in Figure 1). In a corresponding way, assume that t = Θ,t + n,t , where n,t is iid zero mean normal with variance 2 Θ n , and where u n,t and n,t are correlated with covariance Θ n U n . Following Lande (2009), we may assume that juvenile birds of generation t are exposed to the environment u t− 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 t = − 2 u t − 10 , which implies a linear relationship Θ,t = − 2 U,t − 10 , variances

| 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 t = 31 to 60 with population size N = 100 are given in  (9), there are in all six parameter values to be estimated (while â 31 follows from Equation (2a) with ŷ 31 set to 0). In the optimizations, the initial values of Ĝ bb , Ĝ ab , ̂ 2 v , ̂ 2 , and b 31 were set to 0, while the initial value of û ref was set to 10 (when û ref was a free variable). The true value Ĝ aa = 0.025 was used, such that estimates of G bb , G ab , 2 v , and 2 are found relative to G aa = 0.025. Table 1 presents  The results for Case 1 and Case 2 are very much improved with population size N = 10, 000, while the standard errors in Case 3 were only marginally improved by an increased population size (Appendix S3). The results for û ref in Case 2 were, for example, improved from 9.92 ± 0.66 to 9.95 ± 0.19. Table 1, a large error in the assumed reference environment û ref results in large errors in predicted changes in a t over 30 generations (Case 3). Table 2 shows prediction results for more moderate errors in û ref , as well as for û ref = 11.75 (expected mean value in optimization data), and they are all in accordance with Equation (7). Figure 5 shows predicted mean values ŷ t , â t , and b t , as compared to true mean values y t , a t , and b t , for Case 1 and Case 3 in Table 1. Note that y 31 =ŷ 31 is set to 0, and the zero point is thus not the same as in Figure 4, Panel a.

| DISCUSS ION
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 u ref is defined, and that this should not be done arbitrarily. Instead, it should be based on the environment u 0 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 as- In the toy example, the mean plasticity slope b t is predicted quite well also when there is a large error in û ref , and one reason for this is that an error in Ĝ ab to some extent compensates for the error in û ref (Table 1, Case 3). This does not, however, prevent a large error in the predicted mean intercept â t (Equation (7)). In theoretical studies, it is often assumed that the environmental variable is scaled such that u ref = u 0 = 0 (Chevin & Lande, 2015;Lande, 2009). This can be done also in databased applications, provided that u 0 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 100 to 10, 000 (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 a t over time ( -10 −5 (23 ± 10) 10 −5 (23 ± 10) 10 −5 (29 ± 11) 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 inputoutput data, and that must inevitably result in prediction errors. An

ACK N OWLED G EM ENTS
The derivation of prediction equations in Appendix S1 was inspired by a review comment by Michael Morrissey, although on a quite different manuscript with a derivation directly based on the Price equation. I thank University of South-Eastern Norway for support and funding.

CO N FLI C T O F I NTE R E S T
The author declares no conflict of interest.  Table 1. True y t values are shown by solid blue lines.

AUTH O R CO NTR I B UTI
All parameter values except û ref and Ĝ aa = 0.025 were initially set to 0, which gave predictions ŷ t,start = cov W i,t , y i,t ∕W t as shown by dashed blue lines. Final predictions ŷ t are shown by blue dots.