RESEARCH ARTICLE EVOLUTION A method to predict the response to directional selection using a Kalman filter Lisandro Miloccoa,1 and Isaac Salazar-Ciudada,b,c,1 Edited by Nils Stenseth, Universitetet i Oslo, Oslo, Norway; received September 29, 2021; accepted May 2, 2022 Predicting evolution remains challenging. The field of quantitative genetics provides predictions for the response to directional selection through the breeder’s equation, Significance but these predictions can have errors. The sources of these errors include omission of traits under selection, inaccurate estimates of genetic variance, and nonlinearities in Evolution is a historical process the relationship between genetic and phenotypic variation. Previous research showed with contingency, where out- that the expected value of these prediction errors is often not zero, so predictions are comes are sensitive to past events. systematically biased. Here, we propose that this bias, rather than being a nuisance, Due to this, predicting evolution can be used to improve the predictions. We use this to develop a method to predict remains challenging. In this paper, evolution, which is built on three key innovations. First, the method predicts change as we propose a method to predict the breeder’s equation plus a bias term. Second, the method combines information from the response to selection that the breeder’s equation and from the record of past changes in the mean to predict change incorporates history. The method using a Kalman filter. Third, the parameters of the filter are fitted in each generation uses tools from quantitative using a learning algorithm on the record of past changes. We compare the method to genetics and combines them with the breeder’s equation in two artificial selection experiments, one using the wing of the fruit fly and another using simulations that include a complex mapping of genotypes information from past evolution to phenotypes. The proposed method outperforms the breeder’s equation, particularly that is extracted through a when traits under selection are omitted from the analysis, when data are noisy, and when combination of signal processing additive genetic variance is estimated inaccurately or not estimated at all. The proposed and learning techniques. This method is easy to apply, requiring only the trait means over past generations. information is underexploited by existing methods quantitative genetics | evolutionary prediction | Kalman filter | breeder’s equation | G matrix to predict evolution but is of great value since it reflects singularities Evolutionary prediction is an important and active field within evolutionary biology of the evolutionary system. (1–4). Aside from its theoretical value, predicting evolution has important applications, We show that this combination such as developing strategies for the persistence of populations amid rapid environmental of information coming from change (5) and designing interventions to control the spread of a disease (6). the time series and quantitative Quantitative genetics is a widely used approach to study and predict short-term evolution of continuous traits (7, 8). The backbone of this theory is the breeder’s equation genetics methods outperforms (9, 10). In its multivariate form, it provides predictions of the change in the mean of a set classical methods in predicting of traits, from one generation to the next, in response to directional selection: the response to selection. Δz̄ i = Gi Pi−1 s i , [1] where Δz̄ i = z̄ i+1 − z̄ i is the vector of change in trait means from generation i to Author affiliations: a Institute of Biotechnology, Helsinki Institute of Life Science, University of Helsinki, 00014 i + 1; Gi and Pi are additive genetic and phenotypic variance–covariance matrices Helsinki, Finland; b Centre de Recerca Matemàtica, Uni- between traits in generation i , respectively; and s i is the vector of selection differentials in versitat Autònoma de Barcelona, 08193 Barcelona, Spain; generation i . and c Genomics, Bioinformatics and Evolution, Departa- ment de Genètica i Microbiologia, Universitat Autònoma A major appeal of the equation is that its elements can be estimated without detailed de Barcelona, Barcelona, Spain knowledge of the genetic architecture and development underlying the focal traits. Indeed, estimates of Gi and Pi can be obtained using only phenotypic data and known genetic relatedness among individuals in a population (11, 12), while estimates of s i need Author contributions: L.M. conceptualized and developed knowledge of individual fitness (8, 13). The simplicity of the equation, however, is the method, with input from I.S.-C.; L.M. performed the achieved at the cost of some assumptions. research and wrote the original draft; and L.M. and I.S.-C. designed the experiments and edited the original draft. The breeder’s equation assumes an infinitesimal model of genetic effects (i.e., a large The authors declare no competing interest. number of loci, each of small effect) or at least a linear parent–offspring regression (and This article is a PNAS Direct Submission. a few additional assumptions) (8, 14). It further assumes linkage equilibrium, Hardy– Copyright © 2022 the Author(s). Published by PNAS. Weinberg proportions, random associations of environmental effects, and that the traits This article is distributed under Creative Commons being analyzed are uncorrelated with other unmeasured traits under selection (2, 8, 13). Attribution-NonCommercial-NoDerivatives License 4.0 Moreover, the equation is local, meaning that the accuracy of the predictions can only be (CC BY-NC-ND). 1 ensured for a single generation (8). To whom correspondence may be addressed. Email: milocco.lisandro@gmail.com or isaac.salazar@uab.cat. When applied to real systems, the assumptions of the breeder’s equation are violated to This article contains supporting information online at some extent. A common violation is the so-called missing character problem, where the https://www.pnas.org/lookup/suppl/doi:10.1073/pnas. particular traits chosen for study do not account for all selection (2, 15). Another violation 2117916119/-/DCSupplemental. arises when the equation is used to predict the response for several generations under Published July 6, 2022. PNAS 2022 Vol. 119 No. 28 e2117916119 https://doi.org/10.1073/pnas.2117916119 1 of 11 the assumption that the statistical parameters remain constant variables of interest using a model of how the variables are expected over these generations. However, the constancy of the G matrix to change and a series of measurements observed over time. Here, is a debated issue (16, 17) since it can be affected by mutation, we adapt it to be used in the prediction of the response to selection. drift, and selection itself (18), and work on nonlinear genotype– In Results, we first develop the KF method in three parts. Part phenotype maps (19) and gene–environment interactions I: The Breeder’s Equation Plus a Bias Term is the introduction of (20, 21) shows that the G matrix can change rapidly even in the extended equation that consists of the breeder’s equation plus a few generations. This issue is aggravated by the fact that the a bias term. Part II: The Kalman Filter is the development of the breeder’s equation assumes that the real statistical parameters Kalman filter for this application. Part III : Learning the Parameters Gi , Pi , and si are known without error; in reality, these are of the Kalman Filter is the explanation of the learning algorithm inaccessible parameters that must be replaced by estimates. While to learn the parameters of the filter at each generation. Later in this is not strictly a violation of the assumptions, this introduces Results, the KF method is used to predict the response to selection uncertainty and possibly biases to the predictions, particularly in two artificial selection experiments. The datasets are used to when Gi is estimated and used in different environments (22) or explore common situations where the assumptions of the breeder’s when relevant effects are not controlled for during the estimation equation are violated to some extent, as explained above. of Gi [e.g., maternal effects (7, 8, 15)]. Indeed, when applied to real systems, violations of the as- Results sumptions of the breeder’s equation can lead to prediction errors (2, 8, 14, 15, 23–25). A notable example is the problem of Part I: The Breeder’s Equation Plus a Bias Term. We want to stasis (2, 26), where no response to selection is observed in a predict the change in the mean of a set of traits between gener- population that both has ample additive genetic variance and ations, Δz̄ i . We propose the following equation consisting of the is under strong directional selection [e.g., body size in several breeder’s equation plus a bias term b i , a vector of length equal to species, including Soay sheep and red deer (26)]. Prediction errors the number of traits: have also been reported in artificial selection experiments when the parent–offspring regression is nonlinear (23, 27) and when Δz̄ i = Gi Pi−1 s i + b i . [2] selection is applied in the direction opposite to the sign of the genetic correlation between the two traits under selection (7). The bias term can be understood as the part of the response to An important feature of the prediction errors when using the selection that is not captured by the breeder’s equation and that breeder’s equation is that their mean over time can be nonzero arises from violations of assumptions. As such, the systematic bias (14, 25), indicating the presence of a systematic bias. For example, is structured, and we expect the bias of generation i − 1 to be if a trait under selection is missing from the analysis, the prediction similar to the bias at generation i (14, 25). using the breeder’s equation can be biased because there is an Here, we propose to estimate the bias term by using measure- indirect effect of selection that is systematically omitted in the ments of the system up to generation i . In principle, one could prediction (26). Moreover, if the G matrix has changed from its estimate the b i−1 as the difference between the realized change in original estimate, predictions will also be biased because the G the mean, Δz̄ i−1 , and the prediction from the breeder’s equation, −1 used for predictions is incorrect. In this way, the total prediction Gi−1 Pi−1 s i−1 . Assuming that the bias changes slowly, one could error (i.e., the difference between the prediction of the breeder’s further assume that b i ≈ b i−1 and obtain an estimation for the equation and the true response to selection) is composed of bias at generation i . The problem with this approach is that two parts. One part of the error is stochastic due to drift and both the breeder’s prediction and the change in the mean for measurement noise. The other part is deterministic, a systematic the trait are measured with noise, which typically is very large. bias. Due to this systematic bias, the error at a given generation This random component of the prediction error is due to several i is informative of the error at generation i + 1. This indicates factors, including drift and sampling, and can be so large as to that there is potential to improve predictions by incorporating obscure the bias estimated using only the past generation (30). this bias if one could retain the information of past generations Furthermore, the stochastic component of the error is indepen- as a “memory.” dent in each generation, so it contains no useful information that Here, we propose a method to predict the response to direc- can be exploited to improve predictions. This component of the tional selection that yields better predictions when some of the error, therefore, has to be separated from the deterministic part assumptions of the breeder’s equation do not hold. The method, (i.e., the bias), which can be used to improve predictions. which we refer to as the KF method, uses the record of the To minimize the stochastic component of the prediction error, means of the traits in past generations to improve predictions. we propose here to use a Kalman filter to estimate Δz̄ i and b i There are three key innovations in the method. First, it uses in each generation. The filter is explained in the next section. a model for the change in the mean of the traits that is the Most importantly, the filter works by efficiently separating the breeder’s equation plus a bias term, which is the term with stochastic component of the error from the bias. To simplify the memory. Second, the method predicts the change in the traits bookkeeping and notation, we will develop the equations of the and the bias in each generation using a Kalman filter (28). filter for each trait separately. We then rewrite Eq. 2 for each The filter integrates the information of the breeder’s equation trait as and the record of past means of the traits, and it efficiently Δi = ΔB i + bi , [3] deals with the stochastic component of the prediction error. where Δi is the change in the mean of a given trait in generation Third, the method incorporates a learning scheme to learn the i , ΔB i is the prediction using the breeder’s equation, and bi is the parameters of the filter required to provide predictions in each bias. In the next section, we show how we use a Kalman filter to generation. estimate the state variables Δ and b in each generation. The Kalman filter is a hallmark of signal processing theory (28, 29) and has a wide variety of technological applications from Part II: The Kalman Filter.The Kalman filter is a general al- navigation of aircrafts to econometrics. The filter is a general gorithm to estimate the value of a set of state variables from algorithm that allows the estimation of the value of a set of noisy measurements (29). To achieve this, the filter integrates 2 of 11 https://doi.org/10.1073/pnas.2117916119 pnas.org two sources of information. First, it uses a model of how we weight to the first term of the equation. The “trust” is quantified expect the state vector to change from one generation to the by the associated error covariance matrices. This together with the next. This makes the algorithm recursive since the estimate of calculation of the gain Ki is explained in Appendix A. the state vector at time i − 1 is used to make an estimate of the As mentioned above, the algorithm is recursive; the estimates state vector at i . This estimate is combined with a second source obtained in generation i − 1 using Eq. 4 are the starting point for of information to make the estimate of the state vector at time the prediction in generation i . We then require initial estimates at i . This second source of information is a set of measurements time i = 1 to begin the recursion. For our state variables, b̂1 = 0, from the system taken at time i that are related to the state and Δ̂0 is the prediction using the breeder’s equation. vector. The filter combines these two sources of information by a weighted average. How the average is obtained is the central part Part III: Learning the Parameters of the Kalman Filter. The of the filter, and it is achieved by calculating a weight matrix that matrices Qi and Ri have to be provided by the user to im- minimizes the error in the estimates (29). Note that both sources plement the filter explained in Part II: The Kalman Filter. Qi of information described above have associated noise summarized is the covariance of the vector η i , and Ri is the covariance of by the covariance matrices Ri and Qi (explained below). These vector v i , which describes measurement noise. These matrices matrices are the parameters of the filter that have to be provided are hard to calculate analytically. For example, the variance in by the user (Part III: Learning the Parameters of the Kalman Filter). the measurement noise for Δ̃i−1 is affected by drift, selection, For this particular application of the Kalman filter, the state measurement, and sampling (8). An added difficulty is that Ri vector at time i is composed of Δi−1 and bi . Note that with the and Qi can change in time. above definitions, estimating Δi gives us a prediction for z̄i+1 Instead of calculating the matrices Ri and Qi analytically, since z̄i+1 = z̄i + Δi . In developing the algorithm below, we use we learn them using the time series of the trait means. Several the symbol ˆ to refer to estimates of the variables (e.g., Δ̂i is methods exist to identify these matrices using data from the time the estimate of the state variable Δi ). We make the assumption series (31–33). Here, we propose a simple method to learn the that the response to directional selection does not show abrupt matrices at generation i using a window of the last L recorded changes from one generation to the next (8). Additionally, in this changes in the mean {Δ̃i−L , . . . , Δ̃i−1 } to learn the values of application we will assume that the bias changes slowly in time. Ri and Qi (similar to the method presented in ref. 33). This In this way, we define Δi−1 = Δi−2 + ηi and bi = bi−1 + ηib , is done by running the filter inside the window with several where η i = (ηi , ηib ) is a vector of small changes that we assume combinations of Ri and Qi . We then calculate the prediction to be normally distributed with mean zero and covariance matrix error of the method in the window for each combination of Ri Qi . Note that the assumption here is not that the state variables are and Qi . The error is calculated as the difference between the constant in time but rather, that they do not show abrupt changes prediction Δ̂i and the true Δi . If the trait mean in each generation from one generation to the next. is measured with error, the true Δi can be estimated by first There are two measurements at time i that we can use to making a linear regression of the means against generations in improve our estimates. We use the symbol ∼ to indicate that the window and then, calculating the per-generation change (11). the variable has been measured with noise. The measurements The combination of Ri and Qi that results in the smallest error are Δ̃B B B i = Δi + vi and Δ̃i−1 = Δi−1 + vi , where we assume is then used to make the actual prediction of interest at time i . B that v i = (vi , vi ) is a vector of Gaussian measurement error with Note that this process is done in every generation i for each time mean zero and covariance matrix Ri . series separately. In this way, the method learns the best Ri and The Kalman filter combines the estimates of the state variables Qi possible for the specific system at time i . in i − 1 (i.e., Δ̂i−2 and b̂i−1 ) and the new measurements (i.e., To learn the matrices, we assume that Ri and Qi are roughly Δ̃B constant inside the window. This sets a limit to how large the i and Δ̃i−1 ) to provide the best possible estimates of the window can be since if the window is too large, then the matrices state variables in generation i (i.e., Δ̂i−1 and b̂i ). Given the may change substantially inside the window. Then, the size of the relationships described above, this is done using the following window should be kept relatively small, making it hard to learn all formula: the elements of the 2 × 2 matrices Ri and Qi (i.e., more elements      B    Δ̂i−1 Δ̂i−2 Δ̃i Δ̂i−2 − b̂i−1 to learn require a larger dataset). To reduce the number of elements = + Ki − . to learn, we make additional simplifications about Ri and Qi , as b̂i b̂i−1 Δ̃i−1 Δ̂i−2 [4] detailed in Appendix B. These simplifications allow us to reduce The first term of the right-hand side is the state vector estimates the identification to a single parameter to be learned from the in step i − 1. The second term is the correction, which is the data in the moving time window. We call this parameter ρi , and product of the matrix Ki and the error. The error is formed by the it summarizes the relationship between the Qi and Ri matrices (Appendix B). The value of ρi will adjust in each generation to difference between the measurements Δ̃B i and Δ̃i−1 and their reflect the amount of stochastic noise in the data inside the time expected values using the estimates at step i − 1 (Appendix A). window, which is of key importance to determine the value of The estimate of the bias is finally used to predict the change at K as explained above. For example, if measurement noise is large, generation i using Eq. 3, as Δ̂i = Δ̃B i + b̂i . then ρi adjusts accordingly by becoming small so that the gain Ki Ki is a 2 × 2 matrix called the Kalman gain, which assigns assigns little weight to the second term of Eq. 4. For the analyses weights to the correction. The calculation of Ki is the key of the in this paper, we use a window of size L for i > L and size L = i filter, and it is done for each i . Ki is a trade-off between the for 2 < i  L (i.e., we use the available generations in the record). confidence we have on the estimate of the state vector at i − 1 Apart from using the window to learn the parameter ρi of the and the confidence we have on our measurements at generation Kalman filter, we also use it to approximate the uncertainty in the i , and it is calculated to minimize the error covariance matrix of predictions using the KF method. To do this, we calculate the SD the estimates (28, 29). If the measurements are to be trusted, then of the residuals of the predictions against the true change inside the gain will give more weight to the second term of Eq. 4. If the the window. We use this as the uncertainty for the predictions of estimates at i − 1 are to be trusted, then the gain will assign more the method. PNAS 2022 Vol. 119 No. 28 e2117916119 https://doi.org/10.1073/pnas.2117916119 3 of 11 Testing the Method. We compare the performance of different to the true changes. As explained in Materials and M ethods, the methods in predicting the response to directional selection using true change for the teeth experiments is directly the measured two artificial selection experiments, one with teeth simulations change. For the fly experiments, the true change is obtained from and one with the wing of the fruit fly (Fig. 1). We compare a quadratic regression of the means (Fig. 1G). The RMSE is a three prediction methods: the multivariate breeder’s equation, general measurement of the goodness of prediction for the whole the method we introduce, and a univariate method based on time series, which is affected by both the accuracy (i.e., bias) and realized heritability (h 2 ). This last method uses the time window the precision (i.e., variance) of the predictors. It is, therefore, the of past generations to estimate the realized h 2 as the slope of main criteria to compare methods used in this work. Note that the regression of the cumulative response to selection against the the RMSE reported here removes the error from the two first cumulative selection differential (11). This realized h 2 is then used generations since at least two generations are required to estimate to predict change at generation i . the first realized h 2 . The performance of the prediction methods is assessed by calculating the error for single-generation predictions obtained Predicting the Response to Selection in Teeth Simulations. The as the relative rms error (RMSE) between the multivariate series teeth artificial selection experiments are in silico simulations of of predictions and the multivariate series of true changes. The evolution in a population. A key feature of these simulations is that relative RMSE is calculated as the squared root of the sum of the the mapping between genetic and phenotypic variation is done us- squared differences between true and predicted changes, relative ing a model of development that produces realistic morphological A E B F C G 6.0 6.5 4 Trait 1 Trait 3 2 3 5.5 1 0 5 10 15 20 25 30 5 10 15 20 Generations Generations D H 0.6 0.08 0.2 0.06 0.04 -0.2 0 5 10 15 20 25 30 5 10 15 Generations Generations Fig. 1. Summary of the artificial selection experiments. A–D correspond to the teeth experiments, and E–H correspond to the fly wing experiments. A shows the tooth morphology and the three landmarks used for the experiments. The five traits are the x and y coordinates for anterior and posterior landmarks and the y coordinate for central landmark. B shows the directions of selection for 4 of the 32 evolutionary simulations as examples. C is the mean of trait 1 in time for simulation 2, and D shows the change in the trait mean. We do not make a regression in C because there is little measurement noise. E shows the morphology and the five landmarks on the wing. There are six phenotypic traits that are obtained after aligning the 10 coordinates of these landmarks using Procrustes superimposition. F shows the direction of selection. G shows the mean of trait 3 for line R1 together with a quadratic fit to the data and its 90% CI (remaining traits and replicate lines are shown in SI Appendix, Fig. S1). H shows the corresponding change in the mean of trait 3. The line with its 90% CI is the true change estimated from the quadratic regression shown in G, while the points are the measured changes with noise. 4 of 11 https://doi.org/10.1073/pnas.2117916119 pnas.org variation (34). Importantly, the genotype–phenotype map of this Fig. 2A summarizes the total error for the different prediction model is known to be complex and lead to biases in the estimation methods in the teeth simulations. Different scenarios are explored, of the response to selection (26). There are a total of 32 simulations characterized by the combinations of m, σ, and L as described with different optima, each of 30 generations (Fig. 1C). Fig. 1 above. The figure shows that the KF method outperforms the shows the tooth morphology and the three landmarks used. The alternatives on average. The KF method is particularly better than x and y coordinates of these landmarks are the five measured the breeder’s equation when Gi and Pi are not updated in every traits. Fig. 1 also shows the response to selection for one trait in generation (i.e., m = 30, regardless of σ and L) since it is able an example simulation. to correct the prediction bias generated by outdated estimates of Because the data are simulated, all conditions are controlled. variance components. In the idealized scenario with m = 1 and This allows us to study the performance of the prediction methods σ = 0, the breeder’s equation is the least biased possible. In this in different scenarios, where specific sources of prediction error are case, the KF method reduces the spread of the prediction error isolated and when sources are interacting with each other. To study toward lower values but cannot reduce the median error with how using incorrect estimates of Gi and Pi affects predictions, respect to the breeder’s equation. This is because the uncertainty we introduce the parameter m, which determines how often the of the method is sometimes larger than the value of the bias. estimates of Pi and Gi are updated (i.e., m = 1 means that Pi When compared with the realized h 2 method, the KF method and Gi are updated in every generation). In the simulations, the is particularly better when selection is measured with noise selection differential si is known in each generation since all indi- (σ = 0.3) and when the size of the time window is larger (L = 15). viduals are measured with no error. However, in many scenarios, This occurs because the KF method can manage noise more there can be large uncertainty in the measurement of selection. efficiently through the Kalman filter and because it is less sensitive This is the case when selection is not measured in every generation to the size of the time window L. The latter results from the or when only a sample of the population is used to estimate fact that the realized h 2 method assumes that the h 2 is constant selection, leading to SEs that can be as large as the measurement inside the time window and therefore, has problems if it changes of selection itself (8, 35). We investigate the sensitivity of the rapidly. The KF method, on the other hand, only assumes that the different predictors to noise in the measurement of selection by parameter ρi of the Kalman filter is constant inside the window, introducing the parameter σ, which defines the distribution of which is a more robust assumption. Further, the KF method noise that we add to si in each generation. For this, we multiply greatly benefits from frequent updates of the variance components the selection differential si in each generation by 1 + u, where u (m = 1), while the realized h 2 method is unable to exploit this is a random number from a normal distribution with mean zero information. An additional limitation of the realized h 2 method and SD σ. Finally, we also try different values of L, the number of is that it is intrinsically univariate, so it cannot use information of past generations in the moving window used by the KF method (to selection on other traits. learn the parameter ρi ) and the realized h 2 method (to estimate Fig. 2B shows the predicted and true changes for some traits in the realized h 2 ). By trying different combinations of the above example simulations for the scenario with m = 30, σ = 0.3, and parameters, we can explore different scenarios and evaluate how L = 5 (the remaining traits are shown in SI Appendix, Fig. S2). well the methods predict. For example, for m = 1 and σ = 0.3, we Note that the difference between the prediction using the breeder’s can evaluate how the methods deal with noise when Gi and Pi are equation and the KF method in each generation i is the bias bi , estimated in each generation, and for m = 30 and σ = 0.3, we can which is shown explicitly for each example in Fig. 2C. further study how the methods perform when the matrices Gi and We further use the teeth artificial selection experiments to study Pi are assumed constant in addition to the stochastic noise in si . the situation where traits that are under selection are omitted from KF method A σ σ 0.3 KF method True B 0.02 0.00 0.05 0.10 0.15 0.00 -0.02 1 2 4 -0.02 30 -0.06 -0.04 0 10 20 30 0 10 20 30 0 10 20 30 Generations Generations Generations C Bias trait 4 Bias trait 1 Bias trait 2 -0.004 0.00 0.01 30 -0.06 -0.01 -0.012 0 10 20 30 0 10 20 30 0 10 20 30 Generations Generations Generations Fig. 2. Comparison of the prediction methods for the teeth experiments. A shows box plots of the prediction error for the 32 simulations using the tooth model for the three prediction methods compared (i.e., the KF method, the breeder’s equation, and the realized h2 method). Different scenarios were explored by modifying the parameters L (the size of the time window), m (the update time of variance component estimates), and σ (the stochastic error in the estimate of the selection differential). B shows examples of time series for the predictions and true changes with m = 30, σ = 0.3, and L = 5. The time series correspond to traits 1, 2, and 4 in simulations 6, 11, and 15, respectively (remaining traits are shown in SI Appendix, Fig. S2). The insets show a scheme of the tooth, with the direction of selection in each case represented with arrows and the trait plotted shown with a bar. C shows the dynamics of the bias for each of the examples given in B. PNAS 2022 Vol. 119 No. 28 e2117916119 https://doi.org/10.1073/pnas.2117916119 5 of 11 the prediction, what is known as the missing character problem. different pedigree depths ranging from 2 to 15, which correspond Specifically, we predict the change in traits 1, 2, and 3 without data to 400 to 3,000 individuals from the control. We also include a from traits 4 and 5. As before, we explore different scenarios by pedigree depth of one, which means assuming that G1 = P1 (i.e., varying the values of m, σ, and L. This allows us to further study that all phenotypic variation is additive genetic). how omitting traits interacts with other sources of bias, such as Fig. 4 shows the predictions for the change in the traits us- using outdated estimates of variance components (i.e., m = 30). ing the KF method, the realized h 2 method, and the breeder’s Fig. 3A shows a summary of the prediction errors of the equation for the first three traits in selection line R2 (remaining different methods in the different scenarios studied. We find that traits and lines are shown in SI Appendix, Fig. S2). A pedigree omission of traits can lead to biases and that the KF method depth of three was used here, and window sizes of 5 and 15 is able to correct the errors to a large extent. Fig. 3B includes were explored. The time series of the predictions also includes the the time series of the same trait in an example simulation, with measured changes in each generation as well as the true change predictions being done for different combinations of m, σ, and that is obtained from the quadratic regression of the response L. The corresponding bias is shown in Fig. 3C. Note that the bias (Materials and M ethods). It can be seen that the KF method becomes larger when in addition to the omission of traits in the yields predictions that are closer to the true change than the other prediction, it is assumed that variance components are constant methods on average. (compare m = 1 and m = 30 in the time series). Further, note We compare the predictions against the measured changes that the realized h 2 method becomes heavily biased for larger without the regression by calculating the cumulative error in gen- window size, while the KF method is robust to these changes eration i as the sum of the differences between the prediction and (compare L = 5 and L = 15). the measured change (i.e., without regression) from generation 2 to i . Cumulative errors are shown in Fig. 4 and confirm that Predicting the Response to Selection in the Wing. The artificial the KF method has cumulative errors closer to zero. Further, note selection experiments in the wing show the full complexity of that the KF method is more robust to changes in the window size the problem of predicting the response to selection in a real L than the realized h 2 method as seen also for the teeth data. population. This is the most common scenario in which the KF Fig. 5A shows the prediction error for the KF method, the method can be applied. There are three replicates with selection realized h 2 method, and the breeder’s equation for different and one control line, all coming from the same base population pedigree depths and with a window size L = 5 (L = 15 gives very and running for 20 generations (Fig. 1, Materials and M ethods, similar results as shown in SI Appendix, Fig. S4). For low pedigree and SI Appendix, Fig. S1). In each generation, 100 males and 100 depths (i.e., inaccurate estimates of G1 and P1 ), the realized h 2 females are measured. Selection is applied on five landmarks of method outperforms the breeder’s equation, while the opposite the wing as shown in Fig. 1E by selecting the 50% of measured is true for high pedigree depths. This is not surprising since the individuals in the direction shown in Fig. 1F. performance of the breeder’s equation is largely affected by how For this experiment, we only calculate the G matrix at the well we can estimate the variance components, while the realized beginning (i.e., G1 ). Since the control line and the selection lines h 2 method is completely independent of the estimates of G1 and all start from the same base population, we use the pedigree and P1 and relies exclusively on past time series data. Note that our phenotypic data of the initial generations of the control to estimate method exploits both sources of information, namely the time G1 . We call the pedigree depth the number of generations of the series and the estimates of the variance components. control line used to estimate G1 . The larger the pedigree depth, Notably, the KF method using a G matrix with a small pedigree the more precise the estimate of G1 . We test the predicting ability depth outperforms the breeder’s equation using a G matrix with a of the different prediction methods using estimates of G1 for large pedigree depth. This is important because experimentally, it A KF method KF method B C Fig. 3. Comparison of the prediction methods for the teeth experiments when traits are missing. A shows the prediction error for the 32 simulations using the tooth model when predictions for traits 1, 2, and 3 are made without information of traits 4 and 5. Errors are shown for different scenarios characterized by different combinations of parameters L (the size of the time window), m (the update time of variance component estimates), and σ (the stochastic error in the estimate of the selection differential). B shows the time series of predictions and true changes for trait 2 in example simulation 32. The insets show a scheme of the tooth, with the direction of selection in each case represented with arrows and the trait plotted shown with a bar. Different combinations of m, L, and σ were used in each plot. The corresponding dynamics of the prediction bias are shown in C. 6 of 11 https://doi.org/10.1073/pnas.2117916119 pnas.org KF method Fig. 4. Comparison of the prediction methods for the first three traits in selection line R2 of the fly wing experiments. The figure shows the time series for the different methods together with the corresponding cumulative errors for window sizes (L) of 5 and 15 generations in Left and Right, respectively. The time series of the predictions are shown on the left of each panel together with the true change (i.e., obtained from the regression of trait means) (Materials and Methods) and the measured change as points. The predicted and true changes are shown with their approximated uncertainties (90% CI for the true change, one SD for the predictions). The cumulative error for each method is calculated in each generation i as the sum of the difference between the prediction and the measured change from generation 2 to i. A pedigree depth of three was used to estimate the variance components. The rest of the traits for line R2 as well as the six traits for lines R1 and R3 are shown in SI Appendix, Fig. S3. is much more expensive to increase the accuracy of the estimate of and show that it outperforms the multivariate breeder’s equa- G1 than to apply the KF method. The latter only requires record- tion and a univariate method based on realized heritability on ing the trait means in past generations, while the former requires average. The method is general and can be applied to virtually phenotypic and relatedness data in particular breeding designs. any evolving system that is under sustained directional selection. We found that there are two sources of prediction bias for the Most importantly, the KF method only requires the record of breeder’s equation in these experiments. First, there is bias asso- means of the trait for past generations, which is relatively easy to ciated with using wrong estimates of G1 and P1 . This prediction collect, at least compared with alternatives like obtaining better error is most evident when using estimates of G1 with low pedigree estimates of Gi . In this way, the method can be applied to a depth and even more when assuming G1 = P1 . Increasing the wide variety of scenarios, especially when the assumptions of the pedigree depth of the estimates can correct much of these error. breeder’s equation are not met, like in the later stages of long- The second source of bias is the fact that Gi changes during the term selection studies when the full set of traits under selection experiment. This leads to possible errors at later stages of the is not known and when Gi cannot be accurately estimated. The experiment if the G matrix estimated for the base population is more the assumptions are violated, the more the KF method will used, even if the estimate of G1 is obtained with high accuracy outperform the breeder’s equation, as shown in Figs. 2, 3, and 5. (i.e., a deep pedigree in this case). We discuss three key elements of the KF method. The first key Both of the errors described above can be seen in Fig. 5B element is the introduction of the bias term in Eq. 2. This is pro- for trait 6 of line R1 using different pedigree depths, namely 1 posed on the grounds of previous theoretical and empirical work (G = P ), 4, and 10. The corresponding time series of the bias that shows that the expected value of the prediction error using is shown in Fig. 5C. A big part of the error is reduced when the breeder’s equation may not be zero (8, 14, 15, 25). The bias is increasing the pedigree depth. However, even when using a precise introduced here as a single term, b, which can be understood as estimate of G1 (pedigree depth of 10), the breeder’s predictions the quantitative effect of violating the assumptions of the breeder’s remains biased toward the end of the experiment. equation, and its value reflects the complex mixture of sources that is specific to each system. This simple way of modeling the bias Discussion allows for improved predictions but makes it hard to disentangle how different sources of prediction error quantitatively contribute We developed a method to predict the response to directional to the total value of the bias. If this separation of the bias is selection by combining the breeder’s equation with data from desirable, however, it could be performed after the experiment the time series. We tested this method, which we refer to as is over and the complete time series of the bias is available. For the KF method, with two sets of artificial selection experiments example, it is possible to perform analyses of the type proposed PNAS 2022 Vol. 119 No. 28 e2117916119 https://doi.org/10.1073/pnas.2117916119 7 of 11 A The method introduced here can be classified as recursive KF method because it forecasts the variables of interest in i using the estimate in i − 1. There has been recent interest in recursive models to make predictions of future evolution (1, 36–39). For example, Nosil et al. (1) fitted an autoregressive model using several years of data of frequency changes of coloration and pattern in a population of stick insects. They examined whether data from early time points in the series could predict data in later time points of the series (similar to what we do using the window of past generations). They were able to successfully predict changes B in frequency for a trait under clear frequency-dependent selection but failed to predict change for a trait under a more complex, unknown form of selection. They conclude that predictability was limited by the understanding of selection. The authors suggest that knowledge of selection could be determinant in improving predictions when using recursive models. The method we propose in this paper does exactly this; it combines a recursive model with knowledge of selection given by the breeder’s equation. Used like C this, the breeder’s prediction contributes the type of information that purely recursive models are lacking. At the same time, purely theoretical models, like the breeder’s equation, are based on sim- plifying assumptions that may miss some of the complexity of the system. The efficient combination of the recursive model, which is data driven, and the breeder’s equation, which is theoretical, is what results in the method proposed here to outperform each Fig. 5. Comparison of prediction methods for different pedigree depths in approach when used separately. the artificial selection experiments of the wing. A shows the relative RMSE of the prediction using the breeder’s equation, the method based on realized A notable finding is that the KF method provides good pre- h2 , and the KF method for the three replicated selection lines (squares dictions even when G is not estimated at all. This is assuming correspond to line R1, circles correspond to line R2, and triangles correspond that G = P . This is shown in Fig. 5 for the fly experiments and to line R3). The G and P matrices were estimated in the first generation with varying pedigree depths. A pedigree depth of one means that it is assumed SI Appendix, Fig. S5 for the teeth experiments. This is important that G = P. B shows the predicted and true time series for trait 6 of R1 for because the P matrix has been used as a proxy of the G matrix different pedigree depths. The corresponding bias is shown in C. for morphological traits, a simplification suggested due to the difficulty in estimating the latter. The simplification that G is by Le Rouzic et al. (36) to analyze time series data, which are proportional to P is known as Cheverud’s conjecture (40–42) based on proposing autoregressive models for the parameters and has been used, for example, in ref. 43 to infer the evolution and keeping the model that maximizes the likelihood of the of patterns of genetic covariation under directional selection. data. When used in this method, assuming G = P still provides good The second key element of the method is the use of a Kalman predictions because the resulting deviations are corrected by the filter. This was necessary to minimize the effect of stochastic noise bias term. Note that even in this case, the information of selection on the estimation of the bias, which can be large (11, 30). Figs. is still exploited, as it enters the predictions through the selection 2 and 3 show that the Kalman filter is very robust to stochastic differential, s i . An important note is that both the breeder’s errors, particularly when compared with the alternative methods. equation and the KF method perform significantly better when The third key element of the KF method is that we use a Gi is estimated than when it is assumed that Gi = Pi (compare window of data to learn the parameter ρi of the filter in each precision 1 and 2 in Fig. 5A). This means that Gi contains useful generation using a learning algorithm. Apart from enabling the information, even when estimated with relatively low precision. method to be used in real time (i.e., during the experiment), it has This improvement does not occur for the realized h 2 method, the important quality that it allows the parameters to change in which is unable to incorporate Information from variance com- time. Moreover, it exploits the dynamical properties of the time ponents and only relies on time series data. series (36), which are specific to the population of interest and The method proposed here was developed under the assump- its singularities. Note that the moving time window is used to tion of directional selection sustained for several generations. The learn the parameter ρi , which is then used by the filter to estimate method is, therefore, limited to this type of selection. However, the bias in each generation. In this way, the time window is not the formalism of the Kalman filter has the potential to be used used to estimate the value of the bias directly as the average of for evolutionary predictions in several other scenarios, particularly the prediction errors using the breeder’s equation in the window. when stochastic noise is a serious issue. For example, the Kalman This is a valid way of estimating the bias, but it provides worse filter could be coupled with an autoregressive model to predict estimates than the method based on the Kalman filter (over 23% evolution under frequency-dependent selection based on previous worse for the fly wing experiments) (SI Appendix, Table S1). Note evidence that such models are adequate for this type of selection that here we propose to reduce the covariance matrices Ri and (1). For the case of fluctuating selection, prediction is severely Qi to the single parameter ρi , allowing for a straightforward limited by the difficulty in obtaining information of how selection learning algorithm that relies on exploring different values of ρ is acting in each generation. Recent efforts (39) have tried to inside the sliding time window. This simple method to estimate map environmental fluctuations to fluctuations in selection since the covariance matrices gives very good predictions in the datasets certain environmental cues, such as temperature, are much easier we explore, but other more complex methods can be used for this to measure than selection itself. By modifying the equations that purpose (31–33). relate the states with the measurements, one could include the 8 of 11 https://doi.org/10.1073/pnas.2117916119 pnas.org information of these other environmental cues to improve the landmarks shown in Fig. 1B were obtained using a semiautomatic landmarking predictions under the Kalman filter formalism similar to the one software (47). In the control line, 50 males and 50 females were chosen randomly introduced here. as parents for the next generation. In the selected lines, the 50 males and Data-driven methods are only becoming more popular in the 50 females with wings with the shortest distance to the optimum morphology future. This change from more classical, theoretical methods is were selected as parents. The distance of each individual to the optimum was fueled by the rapid accumulation of data. The method we propose calculated as the Euclidean distance between the values of the traits in the here is line with this change by combining theory and data. As individuals and the optimal values of the traits. The optimum morphology is suggested by other authors (38), this is a promising future for shown in Fig. 1D, and it is the same for the three lines with selection. The process developing better predictions in evolutionary biology. of image processing and selection was repeated in each generation. Sibling mating was avoided to reduce inbreeding. The process was repeated for a total of 20 generations, equivalent to 4,000 flies per line (16,000 in total). If some of the Materials and Methods formed couples did not produce offspring for the next generation either because one of the parents died or due to infertility, we measured more offspring from Experiment 1: Teeth. We used data of in silico artificial selection experiments other couples to complete the 200 individuals per generation. We also formed on teeth. Details of the simulations are given in previous work, and the data are three extra couples in each generation to provide extra individuals in case some publicly available (19, 26). Briefly, each evolutionary simulation has a popula- of the original 50 couples failed to produce offspring. tion of genotypes. Each genotype is mapped to a tooth morphology through a As mentioned above, we measured the x and y coordinates of five landmarks, deterministic model of tooth development (34, 44). The tooth model recapitu- resulting in 10 traits. The data were aligned by generalized Procrustes least lates the process of development for a tooth, starting from a flat epithelium to squares superimposition. Four degrees of freedom are lost in this process: one a complex three-dimensional morphology. The dynamics of development are to estimate wing size and three to standardize the orientation of wing shapes. determined by the value of a set of parameters that are determined by the Therefore, there are only six independent traits in the data. For these traits to genotype. Traits were measured on each tooth. These were the x and y coordinates be comparable between lines and through the generations, we use the six first of three landmarks located in the three tallest cusps of the tooth (Fig. 1A). In components of a principal component analysis of generation 1 of the control as each generation, once the genotypes of all individuals had been mapped to a reference and project all the data to that space. The resulting six phenotypic their corresponding phenotypes using the tooth development model, selection traits are a linear combination of the original 10 traits that conserves all relevant was applied by choosing 50% of the individuals with morphology closest to the variation in all lines. In this paper, we refer to these six traits as the phenotypic optimum. Each simulation had an optimum shape defined at the beginning, traits. The means of these traits against generations for all four experimental lines which determined the direction of selection (Fig. 1C). Selected parents were are shown in SI Appendix, Fig. S1. paired randomly and produced the next generation of genotypes. Each couple Estimation of variance components and true change. All lines start from the produced four offspring, resulting in a constant population size. Recombination same founding population. We estimate G1 and P1 for this founding population and mutation were included in each generation, and the process was iterated to by pooling the first n generations of the control. For these n generations, we have simulate evolution. There are in total 32 simulations, each with a different selec- the pedigree and phenotypic data. We call n the depth of the pedigree. Here, we tion optimum. Each simulation was run for 30 generations using a population of explore values of n from 2 to 15. REML estimates of G1 and P1 were obtained 300 males and 300 females. using the software WOMBAT (45), and sampling variation was estimated using Estimation of Variance Components and True Change. In each genera- the REML–MVN method (46). The linear, mixed effect model (i.e. animal model) tion, the elements of the breeder’s equation were estimated (i.e., Gi , Pi , and used included sex, generation, and identification of the person measuring as s i ). Variance components were estimated from a half-sibling breeding design fixed effects. using individuals at generation i as the base population (details are in ref. The estimation of the means in each generation inevitably has noise. Noise 25). The animal model used was the simplest possible (i.e., with only addi- arises from the imaging and landmarking process, finite sampling of the popu- tive genetic merit fitted to each individual). Restricted maximum likelihood lation, and drift. Because we focus on directional selection, this noise has to be (REML) estimates of Gi and Pi were obtained using the software WOMBAT (45). removed. We perform a quadratic regression to the 20-generation time series of Sampling variation in the estimation of Gi was accounted for using the REML– the means, which is the standard type of regression used for long-term artificial multivariate normal (REML-MVN) method (46). For each generation, we resam- selection data (8, 49–51). The fitted values are used as Δi , which we call the true pled 100 G and P matrices from this distribution and used them to calculate change. This is compared with the change predicted by the different methods to 100 predicted changes using the selection differential and the breeder’s equa- calculate the prediction error (Fig. 1 F and H). We also compare predictions with tion. We plot the mean and one SD of these predictions. Note that the tooth the measured change in trait means (i.e., without regression) by calculating the development model is deterministic and that there is no measurement error. cumulative prediction error. For generation i, the cumulative error is the sum of Moreover, we have a large sample size. This allows for very precise estimates the differences between the measured change and the predicted change from of Gi and Pi . Due to the fact that there is little measurement noise for the generation 1 to i. population mean in the simulations, the true change was obtained directly as Δi = z̄i+1 − z̄i . This is the amount that we look to predict at generation i Appendix A (Fig. 1 E and G). Here, we will derive the equations of the Kalman filter and show Experiment 2: Fruit Fly Wing. We performed artificial selection experiments how the gain matrix Ki is calculated in each generation. For this, on the wing of the fruit fly Drosophila melanogaster. The starting population we first express the model in matrix notation. The state equation was founded from 250 isofemale lines derived from flies captured during the summer of 2017 in Groningen, the Netherlands by the Billeter laboratory. From is given by each line, 25 females and males were collected and merged to make a large x i = x i−1 + η i , [5] outbred population that was maintained in laboratory conditions. For the initial  T T generation of the experiments, 100 virgin males and 100 virgin females from the where x i = (Δi−1 , bi ) and η i = ηi , ηib . The measure- large population were randomly assigned to one of four lines. Three of these lines ments used in the model are the prediction of the breeder’s were subjected to selection (R1, R2, and R3), with the remaining being a control equation, Δ̃B i , and the most recent measured change in trait without selection (C1). Lines were kept at 25 ◦ C with alternating 12-h light and dark cycles during the experiment. means, Δ̃i−1 . ΔB i is related to the states via the following In each generation, 100 males and 100 females were collected as virgins. relationship: The left wing of each collected, anesthetized fly was taken by the automatic system known as the WingMachine (47, 48). The x and y coordinates of the five i = Δi − bi = Δi−1 + ηi+1 − bi . ΔB [6] PNAS 2022 Vol. 119 No. 28 e2117916119 https://doi.org/10.1073/pnas.2117916119 9 of 11 The measurement equation is then Appendix B y i = C x i + B η i+1 + v i [7] To implement the Kalman filter given in Appendix A, we need the matrices Qi and Ri∗ . These parameters are, however, gen-  T  T erally unknown and have to be identified from the data. Here, with y i = Δ̃B i , Δ̃i−1 , v i = viB , vi , and we use a simple method based on exploring different values of the parameters and keeping the ones that result in the smallest     prediction error in a moving time window. For this, we reduce the 1 −1 1 0 C= ;B= ; [8] matrices to a single parameter ρi by assuming that the matrices 1 0 0 0 are diagonal with equal elements in the diagonal: that is, Qi = the state-space model is determined by [5] and [7]. The objective qi I and Ri∗ = ri∗ I , where I is the 2 × 2 identity matrix. The of the Kalman filter is to provide estimates of the states at reduction is justified by the fact that the noise in the measure- generation i (x̂ i ) using the previous estimate (x̂ i−1 ) and new ments, as well as in the states, is in the same units and of the measurements (y i ). Particularly, the Kalman filter provides the x̂ i same order of magnitude. We further assume that these values are that minimizes the variance of the error e i = x i − x̂ i denoted by constant inside the window: that is, qi ≈ qi+1 and ri∗ ≈ ri+1∗ . Φi = E[e i e T If we use these definitions, we can rewrite the equations from i ]. For this, it uses the expression Appendix A as x̂ i = x̂ i−1 + Ki (y i − C x̂ i−1 ), [9] Ki = (Φ∗i−1 + ρi I )C T (C (Φ∗i−1 + ρi I )C T + I )−1 , [16] which depends on the matrix Ki . To find the Ki that minimizes x̂ i = x̂ i−1 + Ki (y i − C x̂ i−1 ), [17] Φi , we first need to obtain an expression for the error. Using [5], Φ∗i = (I − Ki C )(Φ∗i−1 + ρi I ), [18] [7], and [9], we get where we define Φ∗i = Φi /ri∗ and ρi = qi /ri∗ . For each gener- e i = x i − x̂ i ation i and for a window size of L, the Kalman filter is run     for generations k ∈ {i − L, . . . , i − 1} with different values of I ηi = (I −Ki ) e i−1 + , ρi and using Eqs. 16–18. Then, ρi is chosen as the value that C C η i + B η i+1 + v i minimizes the mean square error between the predictions Δ̂k and [10] the true value Δk . For the teeth simulations, Δk is directly the where I is the identity matrix. Considering ηi , vi and ei−1 measured change in trait mean (i.e., Δ̃k ) since there is very little independent, the expected value of the cross-products between measurement noise for the trait means. For the wing data, a better e i−1 and both η i and v i vanishes. Then, we get estimate of Δk is obtained by first making a linear regression of the means of the traits inside the window and then calculating    T the change in trait mean as the slope of the regression [11]. For I I Φi = (I −Ki ) Φi−1 the first two iterations of the algorithm, the window is too small C C to calculate ρi , so we set it to zero. The initial conditions are     x̂ 1 = (Δ̃B T ∗ Qi Qi C T I 1 , 0) and Φ1 = 0. + , [11] CQi CQi C T + Ri∗ −KiT Data Availability. The data from in-silico experiments using the tooth model is available in Dryad at https://doi.org/10.5061/dryad.9cnp5hqdr (52). The R script where Ri∗ = BQi+1 B T + Ri , Qi = E[η i η T i ], and Ri = to apply the method and the data from the D. melanogaster experiments are E[v i v T ]. We want to find K i such that Φ i is minimized. stored in GitHub at https://github.com/millisan/Learning-from-mistakes (53). i This is a convex quadratic minimization problem with a unique solution that can be obtained, for example, by using the method ACKNOWLEDGMENTS. We thank Ruben H. Milocco for advice on applying the Kalman filter, David Houle and his team for instructions on how to perform artifi- of completing the squares (ref. 29, pp. 430–433). The solution is cial selection experiments with the fruit fly using the WingMachine, Pınar Güler given by Kohlmeier from the Billeter team for providing the flies from the wild population, and Jukka Jernvall for useful comments. We also thank Niko Björkell, Shuja Ki = (Φi−1 + Qi )C T (C (Φi−1 + Qi )C T + Ri∗ )−1 , [12] Wahid, Ardit Fejzullahi, and Sara Goncalves for help carrying out the experiments. x̂ i = x̂ i−1 + Ki (y i − C x̂ i−1 ), [13] The research was funded by Academy of Finland Grant 315740, the Integrative Φi = (I − Ki C )(Φi−1 + Qi ). [14] Life Sciences Graduate Program of the University of Helsinki, Finnish Cultural Foundation Grant 00210730, Spanish Ministry of Science and Innovation Grant PGC2018-096802-B-I00, and the Spanish State Research Agency through the From this recursion, we get estimates of the states Δ̂i−1 and b̂i . Severo Ochoa and Marı́a de Maeztu Program for Centers and Units of Excel- Finally, the prediction at generation i is given by lence. We also thank Centres de Recerca de Catalunya Programme/Generalitat de Catalunya for institutional support. Drosophila work was supported by the Hi-Fly Δ̂i = Δ̃B i + b̂i . [15] core facility funded by the Helsinki Institute of Life Science and Biocenter Finland. 1. P. Nosil et al., Natural selection and the predictability of evolution in Timema stick insects. Science 5. R. Gomulkiewicz, R. G. Shaw, Evolutionary rescue beyond the models. Philos. Trans. R. Soc. Lond. B 359, 765–770 (2018). Biol. Sci. 368, 20120093 (2013). 2. R. G. Shaw, From the past to the future: Considering the value and limits of evolutionary prediction. 6. S. Cobey, Modeling infectious disease dynamics. Science 368, 713–714 (2020). Am. Nat. 193, 1–10 (2019). 7. D. A. Roff, A centennial celebration for quantitative genetics. Evolution 61, 1017–1032 (2007). 3. A. Le Rouzic et al., Unidirectional response to bidirectional selection on body size II. Quantitative 8. B. Walsh, M. Lynch, Evolution and Selection of Quantitative Traits (Oxford University Press, 2018). genetics. Ecol. Evol. 10, 11453–11466 (2020). 9. J. L. Lush, Animal Breeding Plans. (Iowa State College Press, ed. 2, 1937). 4. M. Wortel et al., The why, what and how of predicting evolution across biology: From disease to 10. R. Lande, Quantitative genetic analysis of multivariate evolution, applied to brain: Body size biotechnology to biodiversity. EcoEvoRxiv [Preprint] (2021). https://ecoevorxiv.org/4u3mg/ (Accessed allometry. Evolution 33, 402–416 (1979). 28 September 2021). 11. M. Lynch, B. Walsh, Genetics and Analysis of Quantitative Traits (Oxford University Press, 1998). 10 of 11 https://doi.org/10.1073/pnas.2117916119 pnas.org 12. L. E. B. Kruuk, Estimating genetic parameters in natural populations using the “animal model.” Philos. 33. Y. Huang, F. Zhu, G. Jia, Y. Zhang, A slide window variational adaptive Kalman filter. IEEE Trans. Circuits Trans. R. Soc. Lond. B Biol. Sci. 359, 873–890 (2004). Syst. II Express Briefs 67, 3552–3556 (2020). 13. R. Lande, S. J. Arnold, The measurement of selection on correlated characters. Evolution 37, 34. I. Salazar-Ciudad, J. Jernvall, A computational model of teeth and the developmental origins of 1210–1226 (1983). morphological variation. Nature 464, 583–586 (2010). 14. S. H. Rice, Evolutionary Theory: Mathematical and Conceptual Foundations (Sinauer Associates, 2004). 35. M. B. Morrissey, In search of the best methods for multivariate selection analysis. Methods Ecol. Evol. 15. B. Pujol et al., The missing response to selection in the wild. Trends Ecol. Evol. 33, 337–346 (2018). 5, 1095–1109 (2014). 16. S. J. Steppan, P. C. Phillips, D. Houle, Comparative quantitative genetics: Evolution of the G matrix. 36. A. Le Rouzic, D. Houle, T. F. Hansen, A modelling framework for the analysis of artificial-selection time Trends Ecol. Evol. 17, 320–327 (2002). series. Genet. Res. 93, 155–173 (2011). 17. J. D. Aguirre, E. Hine, K. McGuigan, M. W. Blows, Comparing G: Multivariate analysis of genetic 37. M. Rescan, D. Grulois, E. Ortega-Aboud, L. M. Chevin, Phenotypic memory drives population growth variation in multiple populations. Heredity 112, 21–29 (2013). and extinction risk in a noisy environment. Nat. Ecol. Evol. 4, 193–201 (2020). 18. A. G. Jones, S. J. Arnold, R. Bürger, Stability of the G-matrix in a population experiencing pleiotropic 38. P. Nosil et al., Ecology shapes epistasis in a genotype-phenotype-fitness map for stick insect colour. mutation, stabilizing selection, and genetic drift. Evolution 57, 1747–1760 (2003). Nat. Ecol. Evol. 4, 1673–1684 (2020). 19. L. Milocco, I. Salazar-Ciudad, The evolution of the G-matrix under nonlinear genotype-phenotype 39. M. Rescan, D. Grulois, E. O. Aboud, P. de Villemereuil, L. M. Chevin, Predicting population genetic maps. Am. Nat. 199, 420–435 (2022). change in an autocorrelated random environment: Insights from a large automated experiment. PLoS 20. C. M. Sgrò, A. A. Hoffmann, Genetic correlations, tradeoffs and environmental variation. Heredity 93, Genet. 17, e1009611 (2021). 241–248 (2004). 40. J. M. Cheverud, A comparison of genetic and phenotypic correlations. Evolution 42, 958–968 (1988). 21. C. W. Wood, E. D. Brodie III, Environmental effects on the structure of the G-matrix. Evolution 69, 41. S. M. Sodini, K. E. Kemper, N. R. Wray, M. Trzaskowski, Comparison of genotypic and phenotypic 2927–2940 (2015). correlations: Cheverud’s conjecture in humans. Genetics 209, 941–948 (2018). 42. A. C. Love et al., Evolvability in the fossil record. Paleobiology, 1-24 (9 November 2021). 22. M. Pigliucci, Genetic variance–covariance matrices: A critique of the evolutionary quantitative genetics https://doi.org/10.1017/pab.2021.36. research program. Biol. Philos. 21, 1–23 (2006). 43. A. P. A. Assis, J. L. Patton, A. Hubbe, G. Marroig, Directional selection effects on patterns of phenotypic 23. A. Gimelfarb, J. H. Willis, Linearity versus nonlinearity of offspring-parent regression: An experimental (co)variation in wild populations. Proc. R. Soc. B Biol. Sci. 283, 20161615 (2016). study of Drosophila melanogaster. Genetics 138, 343–352 (1994). 44. E. Harjunmaa et al., Replaying evolutionary transitions from the dental fossil record. Nature 512, 24. M. B. Morrissey, L. E. B. Kruuk, A. J. Wilson, The danger of applying the breeder’s equation in 44–48 (2014). observational studies of natural populations. J. Evol. Biol. 23, 2277–2288 (2010). 45. K. Meyer, WOMBAT: A tool for mixed model analyses in quantitative genetics by restricted maximum 25. L. Milocco, I. Salazar-Ciudad, Is evolution predictable? Quantitative genetics under complex likelihood (REML). J. Zhejiang Univ. Sci. B 8, 815–821 (2007). genotype-phenotype maps. Evolution 74, 230–244 (2020). 46. D. Houle, K. Meyer, Estimating sampling error of evolutionary statistics based on genetic covariance 26. J. Merilä, B. C. Sheldon, L. E. B. Kruuk, Explaining stasis: Microevolutionary studies in natural matrices using maximum likelihood. J. Evol. Biol. 28, 1542–1549 (2015). populations. Genetica 112-113, 199–222 (2001). 47. D. Houle, J. Mezey, P. Galpern, A. Carter, Automated measurement of Drosophila wings. BMC Evol. 27. J. S. Heywood, An exact form of the breeder’s equation for the evolution of a quantitative trait under Biol. 3, 25 (2003). natural selection. Evolution 59, 2287–2298 (2005). 48. J. G. Mezey, D. Houle, The dimensionality of genetic variation for wing shape in Drosophila 28. R. E. Kalman, A new approach to linear filtering and prediction problems. J. Fluids Eng. 82, 35–45 melanogaster. Evolution 59, 1027–1038 (2005). (1960). 49. E. J. Eisen, Long-term selection response for 12-day litter weight in mice. Genetics 72, 129–142 29. K. Åström, B. Wittenmark, Computer-Controlled Systems: Theory and Design (Prentice-Hall, (1972). Englewood Cliffs, NJ, ed. 3, 1997). 50. J. J. Rutledge, E. J. Eisen, J. E. Legates, An experimental evaluation of genetic correlation. Genetics 75, 30. C. Pélabon et al., Quantitative assessment of observed versus predicted responses to selection. 709–726 (1973). Evolution 75, 2217–2236 (2021). 51. P. Grassini, K. M. Eskridge, K. G. Cassman, Distinguishing between yield advances and yield plateaus 31. M. Laine, “Introduction to dynamic linear models for time series analysis” in Geodetic Time Series in historical crop production trends. Nat. Commun. 4, 2918 (2013). Analysis in Earth Sciences, J.-P. Montillet, M.S. Bos, Eds. (Springer, Cham, Switzerland, 2020), pp. 52. L. Milocco, Data from “Is evolution predictable? Quantitative genetics under complex 139–156. genotype-phenotype maps. Dryad. https://doi.org/10.5061/dryad.9cnp5hqdr. Deposited 18 32. J. Dunı́k, O. Straka, O. Kost, J. Havlı́k, Noise covariance matrices in state-space models: A survey and December 2019. comparison of estimation methods—part I. Int. J. Adapt. Control Signal Process. 31, 1505–1543 53. L. Milocco, millisan/Learning from mistakes. GitHub. https://github.com/millisan/Learning- (2017). from-mistakes. Deposited 7 February 2022. PNAS 2022 Vol. 119 No. 28 e2117916119 https://doi.org/10.1073/pnas.2117916119 11 of 11