Modeling Adaptive and Nonadaptive Responses of Populations to Environmental Change

Understanding how the natural world will be impacted by environmental change over the coming decades is one of the most pressing challenges facing humanity. Addressing this challenge is difficult because environmental change can generate both population-level plastic and evolutionary responses, with plastic responses being either adaptive or nonadaptive. We develop an approach that links quantitative genetic theory with data-driven structured models to allow prediction of population responses to environmental change via plasticity and adaptive evolution. After introducing general new theory, we construct a number of example models to demonstrate that evolutionary responses to environmental change over the short-term will be considerably slower than plastic responses and that the rate of adaptive evolution to a new environment depends on whether plastic responses are adaptive or nonadaptive. Parameterization of the models we develop requires information on genetic and phenotypic variation and demography that will not always be available, meaning that simpler models will often be required to predict responses to environmental change. We consequently develop a method to examine whether the full machinery of the evolutionarily explicit models we develop will be needed to predict responses to environmental change or whether simpler nonevolutionary models that are now widely constructed may be sufficient.


Introduction
Ecosystems-from the deep ocean to the high arctic, from deserts to tropical forests-are responding to environmental change. Understanding and predicting these responses is one of the most pressing issues currently facing humanity. For this reason, in the past quarter century, there has been considerable interest in developing ways to understand how the natural world will be affected by environmental change (Ives 1995;Bossdorf et al. 2008;Gilbert and Epel 2009;Wiens et al. 2009;Lavergne et al. 2010;Dawson et al. 2011;Hoffmann and Sgrò 2011). We introduce a new general approach combining insights from structured population modeling and evolutionary genetics that allows us to examine how adaptive evolution and plasticity contribute to the way that populations-and, consequently, the ecosystems in which they are embedded-respond to environmental change.
In order to understand how evolution and plasticity contribute to population responses to environment change, it is necessary to appreciate how different levels of biological organization-alleles, genotypes, phenotypes, populationsare linked, as well as feedbacks between the different levels. First, evolution is defined as a change in allele frequencies (Charlesworth 1994). Allele frequencies change as a direct consequence of changes in the frequencies of the genotypes the alleles occur in, and genotype frequencies can change with a change in the distribution of the phenotypes they code for (Fisher 1930). The dynamics of phenotypic trait distributions are determined by differential birth, death, development, and inheritance rates across phenotypic trait values, where inheritance is defined in the broad sense as the map between parental and offspring phenotypes (Easterling et al. 2000;Rees et al. 2014). Given these links between different levels of biological organization, there can be a cascading dynamic at the level of the phenotype, the genotype, and the allele caused by differences in the demography of individuals with different phenotypic trait values (Lynch and Walsh 1998). Another consequence of this variation is the ecology of the system: population dynamics are an emergent property of who lives, who breeds, and with whom, as are the dynamics of the community and ecosystem the population is embedded within (Caswell 2001).
Although the cascading ecological and evolutionary consequences of variation in demographic rates are relatively straightforward to grasp, the devil is in the details-in particular, how alleles combine to make genotypes, how genotypes influence phenotypes, and how phenotypes influence demographic rates (Coulson et al. 2011). The rate and direction of evolution depends on how these links influence the relative fitness of each allele within the population (Charlesworth 1994). The challenge is that these links are often complicated, particularly for complex phenotypic traits such as body size that are routinely measured by field biologists. The complexity arises not only because large gene networks and multiple cell types can contribute to the phenotype but also because the environment makes a contribution too via plasticity, defined as change in a phenotypic trait distribution that is not caused by genetic change (Baldwin 1896;Gavrilets and Scheiner 1993;Lande 2009).
The environment can be partitioned into biotic and abiotic components (Berryman 2002). The biotic component captures the sizes and structures of the population of the focal species and of all other species with which it interacts. The abiotic environment includes weather, mineral, and water available. The biotic and abiotic environments can influence one another, although the influence of the biotic environment on the abiotic environment typically plays out over geological timescales (one exception being man-made climate change).
The biotic and abiotic environment can influence the map both between genotype and phenotype (Baldwin 1896) and between phenotype and demographic rates (Link et al. 2002). Put another way, demographic rates are a function of phenotype-by-environment interactions, and phenotypic traits are a function of genotype-by-environment interactions. For quantitative phenotypic traits, genotype-by-environment interactions can often usefully be understood by treating the phenotype as consisting of a genetic and an environmental component, with the environmental component determined by aspects of the current and past biotic and abiotic environments (Falconer 1960;Lande 1982;Cheverud et al. 1983). The environmental component of the phenotype can capture phenotypic change caused by individuals altering their physiology, metabolism, behavior, or levels of gene expression. We use the term "epigenetic" to refer to any process that does not involve genetic change that is captured by the dynamics of the environmental component of the phenotype. The biotic and abiotic environment can also influence the generation of new alleles, for example, via retroviral insertions into the germ line of their hosts or via ultraviolet radiation (Salter et al. 1987;Kanjilal et al. 1993). In figure 1a, we depict how different levels of biological organization are linked and feed back to influence one another.
How can this view of biology be used to inform how populations respond to environmental change? Environmental change occurs when the biotic or abiotic environment changes. Biotic changes can result from the arrival of a new species or an extinction within the ecosystem or from evolution. In order to capture such change and to model the links between alleles and demographic rates described above and in figure 1a, it is necessary for models to incorporate (1) the genotype-phenotype map at birth, (2) how the phenotype develops, (3) how the phenotype influences survival at each developmental stage, (4) the population's mating system, and (5) patterns of mate choice based on the phenotype, as well as how these mate choice patterns influence (6) reproductive success, (7) the distribution of genotypes among offspring, and (8) how all these processes result in change in allele frequency and population size from one generation to the next. Processes 1-6 (and consequently also 8) can be influenced by the biotic or abiotic environment. Integral projection models (IPMs) provide a very flexible structured modeling framework that allow each of these processes to be simultaneously modeled (Easterling et al. 2000;Coulson 2012;Merow et al. 2014).
IPMs project the dynamics of phenotype distributions as a function of expected survival and reproduction, the way the phenotype develops, and the distribution of offspring phenotypes (Easterling et al. 2000). Numerous quantities of interest to ecologists and evolutionary biologists describing life-history, population dynamic, and phenotypic traits can be calculated from IPMs (Childs et al. 2003;Ellner and Rees 2006;Coulson et al. 2010Coulson et al. , 2011Rees et al. 2014;Steiner et al. 2014Steiner et al. , 2012Vindenes and Langangen 2015). They consequently offer great potential to study ecological and evolutionary responses to environmental change (Coulson et al. 2011). However, most IPMs to date have been restricted to phenotypic variation in that they do not include genotypephenotype maps (Merow et al. 2014). A small number of evolutionarily explicit IPMs that do include these maps have been developed. For example, Coulson et al. (2011) used IPMs to track the distribution of body size and coat color in wolves, where coat color was determined by genotype at a single biallelic locus. Barfield et al. (2011) and Childs et al. (2016) developed IPMs of quantitative characters determined by a large number of unlinked loci of small effect. However, none of these models incorporates plasticity or different genetic influences on the phenotype at different ages, and these omissions limit their utility in predicting how populations will be influenced by environmental change (Chevin 2015).
The aim of this article is to introduce a general framework to allow prediction of how populations respond to environmental change. We do this by developing IPMs of the bivariate distribution of a phenotype split into its genetic and environmental components. The models incorporate different development and inheritance rules for each component of the phenotype. We develop and illustrate our framework using simple models. Our models reveal new insights into the way that plasticity can influence evolution while also allowing us to retrieve key findings from evolutionary genetics that are already known.

Modeling Approach
We refer to functions as f(...), where the dots inside parentheses define the variables on which the function f operates. Parameters of a function are referenced by the same letter as the function, with subscripts defining the variable they influence. For example, a parameter f Z represents a parameter of function f that operates on variable Z. We reserve I a b c d e f g h Figure 1: a, Linkages and feedbacks in biology. Evolution is defined as the change in allele frequencies but is often inferred from the dynamics of genotypes and phenotypes. Research into links between alleles and genotypes-and particularly between genotypes and phenotypes-often focuses on mechanism (red arrows). Differential survival and reproduction and patterns of mating determine (1) the dynamics of phenotypes, genotypes, and alleles and (2) population, community, and ecosystem dynamics (purple arrows). Ecological dynamics determine the biotic environment that, along with the abiotic environment, can influence the generation of new alleles as well as the maps between genotype and phenotype and between phenotype and demographic rates. b, Integral projection models track the dynamics of the phenotype distribution from t (black line) to t 1 1 (blue line for the intercept of functions and a for age. Age is included only in models for species with overlapping generations. Following standard convention for IPMs (Coulson 2012;Merow et al. 2014;Rees et al. 2014;Ellner et al. 2016), we use primes to indicate the character value of an individual at the end of a time step. For example, this allows us to show how an individual with phenotype Z can develop over a time step to a potentially different phenotype Z 0 or how a parent with genotype G can produce an offspring with genotype G 0 . There are, of course, other notational conventions that could achieve the same objective, and we recognize that primes are used differently in evolutionary genetics; our notation is chosen to make clear how evolutionary processes can be included in the IPM framework. The definitions of all variables and functions are summarized in table 1. Our starting point is a widely used phenotypic modeling approach that many readers will be familiar with (Coulson 2012;Merow et al. 2014;Rees et al. 2014). We then extend this approach by developing dynamic models of the phenotype decomposed into its genetic and environmental components. We start with a two-sex IPM that captures all demographic processes that can contribute to the dynamics of phenotypes-survival, recruitment, development, inheritance, and mating patterns (Coulson et al. 2011;Schindler et al. 2013Schindler et al. , 2015Traill et al. 2014a)-and that iterates for-ward the distribution of the phenotype at time t: N(Z, t) ( fig. 1b).
The model consists of two equations-one for females and one for males-with each equation consisting of two additive components (Schindler et al. 2013). The first component deals with survival and development of individuals already within the population, and the second component deals with reproduction and the generation of phenotypes among newborns entering the population. We assume a prebreeding census such that survival occurs before development and recruitment before inheritance: : : : : Distribution of X at time t. Note that this is an abundance distribution (not a probability distribution): Ð b a N(X )dX is the number of individuals with characters in the interval [a, b], and the integral of N(X , t) over the full range of X gives the total population size at time t. N(A, E, t) Bivariate distribution of the additive genetic and environmental components of the phenotype at time t. Z p z(G, E) Function describing the phenotype as a function of its genetic and environmental components.

S(Z, t)
Survival function. Describes the expected association between Z and survival between t and t 1 1. Only used in agestructured models.

R(Z, t)
Recruitment function. Describes the expected association between Z and the number of offspring produced between t and t 1 1 that survive to recruit into the population at time t 1 1.
Inheritance function. Describes the expected probability of a reproducing individual with character value X at t producing an offspring with character value X 0 at t 1 1 when it recruits to the population. D(E 0 jE, t) Development function. Describes expected probability of a surviving individual with E at t expressing E 0 at t 1 1. Only used in age-structured models. ð1Þ where N f (Z 0 , t 1 1) and N m (Z 0 , t 1 1) are distributions of phenotypes Z 0 in females and males, respectively, at time t 1 1; D f (Z 0 jZ, v, t) and D m (Z 0 jZ, v, t) are the probability of the phenotype developing from Z to Z 0 in females and males, respectively, between t and t 1 1 as a function of environmental drivers v; S f (Z, v, t) and S m (Z, v, t) are survival functions for females and males from t to t 1 1, including effects of phenotype and environmental drivers v; s is the birth sex ratio measured as the proportion of female offspring produced; H f (Z 0 jZ m , Z f , v, t) and H m (Z 0 jZ m , Z f , v, t) describe the probabilities of parents with phenotypes Z m and Z f producing male and female offspring, respectively, with phenotype Z 0 as a function of environmental drivers v at time t; M(Z m , Z f , t) captures the rate of mating between a male with phenotype Z m and a female with phenotype is the expected litter size, given a mating between a male and a female with phenotypes Z m and Z f in environment v at time t; and C N f N m is a normalization constant that is used to specify the mating system. In theory, it could be combined with the mating function, but we follow the notation of Schindler et al. (2013). C N f N m can be used to capture a range of mating systems. For example, if we follow Schindler et al. (2013) and write this adds a minimum size at which females can reproduce Z f (min) . Depending on the mating behavior of the species, C N f N m can be modified in various ways. For example, it can easily be altered such that the number of birth events is determined by the number of the rarer sex, as in monogamous species. Mate choice can be influenced by specifying different functions for M(Z m , Z f , t). Schindler et al. (2013) demonstrate how it can be specified for random mating, assortative mating, disassortative mating, and size-selective mating. In phenotypic IPMs, the phenotypic development functions are usually Gaussian probability functions (Easterling et al. 2000); for example, The functions m D (Z, v, t) and V D (Z, v, t), respectively, describe the expected value of Z 0 given Z and v at time t and the standard deviation around m D (Z, v, t). The Gaussian form can also be used for inheritance functions H(Z 0 jZ, v, t) with functions m H (:::) and V H (:::). The two-sex IPM described above is not evolutionarily explicit because it does not include mechanistic rules for genetic inheritance. We now take this phenotypic model and extend it to be evolutionarily explicit. We do this by writing the phenotype as a function of genetic G and environmental E components Z p z(G, E). We assume that Z is a quantitative phenotype (i.e., measured in integer or real values). The genotypic value G and environmental value E describe the numerical contributions of the genetic and environmental components of the phenotype to an individual's phenotypic trait value. A simple map can consequently be written Z p G 1 E (Falconer 1960).
G is determined by genotype g. When the map between g and G is additive, the dynamics of g and G are identical (Falconer 1960). This means that the dynamics of alleles are identical to the dynamics of genotypes in which they occur. In contrast, when alleles interact-either at a locus (dominance) or across loci (epistasis)-the map between g and G is not additive, and the dynamics of G are not identical to the dynamics of g (Fisher 1930). In classical quantitative genetics, it is assumed that the map between g and G is additive (Falconer 1960). Under these assumptions, it is not necessary to track the dynamics of g, but evolution can be investigated by modeling the dynamics of just G. When the map is additive, we refer to the genetic component of the phenotype as a breeding value and denote it A.
In classical population genetics, when the contribution of dominance and epistasis to evolution are often a key focus, it is necessary to track the dynamics of g and calculate G from each g. The map between G and the phenotype Z is often assumed to be one-to-one (Hartl et al. 1997). In contrast, in quantitative genetics, the environment can influence the map between A and Z by influencing the value of the environmental component of the phenotype, E (Falconer 1960). E can take different values in different individuals and can vary within individuals throughout life. The dynamics of the phenotype may not consequently represent the dynamics of the genotypic value A. Statistical quantitative genetics is concerned with estimating moments of A from Z by correcting for environmental and individual variables that determine E (Kruuk et al. 2008).
The genotype-phenotype map for phenotypic traits measured by biologists in free-living populations is rarely known, and quantitative genetic assumptions are widely adopted (Kruuk et al. 2008). In particular, the infinitesimal model is assumed in which A is determined by a large number of unlinked loci of small additive effect (Fisher 1930). Until we have a better understanding of the genetic architecture of complex traits, this approach is the most powerful available to investigate evolution in the wild (Kruuk et al. 2008). We consequently adopt it here.
We track the joint distribution of the two components of the phenotype N(A, E, t). The utility of this is that we can write expressions to describe the dynamics of each of the components separately, if necessary, before easily combining them to retrieve the dynamics of the phenotype. For Z p A 1 E, we can use a convolution between the two components of the phenotype to construct the phenotype (Barfield et al. 2011).
Phenotypic plasticity and epigenetic inheritance are captured in the dynamics of E. In previous quantitative genetic IPMs, E is a randomly distributed variable that captures developmental noise (Barfield et al. 2011;Childs et al. 2016).
A key contribution of this article is to show how E can be extended to also capture the biotic or abiotic environment as well as signatures of parental A's and E's. E is defined as a function of these factors, and we write E 0 jE, A, v, t to capture the effects of E, A, and the environment v at time t on E 0 .
We now expand terms in our two-sex phenotypic IPM to include the genotype-phenotype map Z p z(A, E). We start with the bivariate distribution of A and E at time t among females that are already within the population at time t: N f (A, E, t). Viability selection now operates on this distribution. Viability selection is a simple multiplicative process describing the expected survival from t to t 1 1 as a function of the phenotype. We can consequently write When it comes to development, A remains fixed throughout life, while E may vary: Recruitment is dealt with in a similar way to survival in that it is a multiplicative process: Note that this is a recruitment-related term of both male and female offspring that is not yet scaled by the normalization factor C N f N m . As with development, inheritance of the genetic and environmental components of the phenotype operates in different ways. For example, once mating pairs have formed and the number of offspring from each mating has been determined, the distribution of offspring genotypes is predictable. We can write the inheritance function for the genetic and environmental components of the phenotype as ::: then The same logic applies to the production of male offspring. We can construct the phenotype from the two components A 0 and E 0 , for example, where Q Z 0 is the set of (A 0 , E 0 ) values satisfying z(A 0 , E 0 ) p Z 0 . For the second integral in equation (8), we have z(A, E 0 ) p Z 0 because the A does not change within individuals and consequently has no prime. The additivity assumption means that models of clonal inheritance can generate very similar predictions to models of two sexes, particularly if both males and females have similar demography. However, clonal models are simpler than two-sex models (Lande 1982). We utilize this consequence of the additivity assumption and initially work with clonal reproduction to examine how the dynamics of A and E influence population and phenotypic trait dynamics and adaptive evolution. We can write a clonal model, The above equations describe the dynamics of a bivariate distribution of the genetic and environmental components of the phenotype. Figure 1b-1g provides graphical examples of how these functions alter the bivariate distribution and, in particular, how development and inheritance rules differ between the environmental and additive genetic components.
To demonstrate these differences, we now focus on developing univariate models of (1) A and (2) E. These models capture limits where all phenotypic variation among individuals is determined by (1) genetic variation and (2) variation in the environmental component of the phenotype. We then combine insights from these univariate model and construct models of the bivariate distribution of A and E. We primarily work with linear functions for three reasons. First, they are easier to interpret and analyze than nonlinear or nonadditive forms. Second, when the environment changes impacting populations, responses-at least in the short term-can be well described with linear or linearized additive models (Cooch et al. 2001). Third, selection-the underpinning of adaptive evolution-is often directional and well described with linear or linearized associations between phenotypic traits and components of fitness (King-ð8Þ solver et al. 2001). Parameters used for all models are provided in section A1 in the appendix (available online), as are expressions to calculate key statistics used to show ecological and evolutionary change from model outputs (sec. A2). Code to produce each figure is available in the Dryad Digital Repository (http://dx.doi.org/10.5061/dryad.4c117; ) as well as GitHub (https://github.com/tncoulson/QG -meets-IPM-figure-code/tree/master).

Adaptive Evolution: The Dynamics of A
Here we start with a simple clonal model of a univariate distribution of A. We go on to show how genetic constraints can be imposed to slow or stop evolution. We then extend this clonal model in two ways: first, we include a multivariate, age-structured distribution of A, and second, we relax the clonality assumption and compare the dynamics of clonal and sexual models. Finally, we introduce a new approximation to describe sexual reproduction and compare its performance with our initial approach.
Genotypes (and hence A) are determined at birth and remain fixed throughout life; neither is influenced by the environment. A consequence of this is that the development function simplifies to a one-to-one map and can be removed from equation (5). We also start by considering clonal reproduction, which means that the inheritance function can also be removed because offspring genotype is identical to parental genotype. The dynamics of A are consequently determined by the survival and reproduction. In these models, as long as there is genetic variation within a population and fitness is a monotonic function of genotype, evolution-defined as E(N(A, t 1 1)) p E(N r (A, t)) ( E(N(A, t)) (where E represents expectations)-will occur.
In our first models we assume nonoverlapping generations, and a linear reproduction function R(A, t) p R I 1 R A A, with expected fitness increasing with the value of A. Over the course of a simulation of 30 generations (model A in section A1.1), the population never achieves an equilibrium structure or growth rate; it grows hyper-exponentially In this model there are two ways to prevent the fitness function from generating change in the location of the distribution. First, the fitness function can take unimodal non-linear forms, such as R(A, t) p R I 1 R A A 1 R A 2 A 2 , with R A 1 0, R A 2 ! 0, and R(A, t) constrained to nonnegative values. This generates stabilizing selection, with the mean breeding value being maintained at the value that maximizes fitness. Eventually, in this model, the breeding value distribution will achieve a trivial equilibrium, a Dirac delta function at this value. Second, continual change in the location of the distribution can be prevented by defining a maximum possible value for A that cannot be exceeded. This captures a genetic constraint in the maximum possible character value; that is, evolution has not evolved a genetic solution to creating a larger breeding value. In our models, this process can be captured by setting the abundance of N(A 1 x, 1) p 0, where x is the maximum possible trait value that evolution can achieve. Selection now pushes the breeding value distribution up to x, again eventually achieving a trivial equilibrium captured by a Dirac delta function where all mass of the distribution is at A p x.
Genetic constraints can also impact the transient dynamics of the breeding value distribution (fig. 2a-2d, red lines). When we impose a genetic constraint (model A in sec. A1.1 with x p 11:5), the genetic variance and skew evolve faster than when no genetic constraint is in place ( fig. 2c, 2d ). These more rapid changes result in a slowing in the evolution of the mean breeding value ( fig. 2b) and of the population growth rate ( fig. 2a).
Genetic covariances between traits can also capture genetic constraints and can also influence the outcome of evolution. We demonstrate this by developing an agestructured model. A now becomes age structured but is still inherited at birth. We construct a multivariate character A describing the breeding values that influence a character at each age (e.g., A1, A2, :::, An for breeding values at ages a p 1, 2, :::, n). If some of the same loci contribute to the genetic components of the character at different ages, there is a genetic covariation across ages. The genetic variances within each age and the covariances between ages can be used to construct a G matrix (Lande 1979). Such age-structured G matrices underpin the character-state approach of quantitative genetics (Lynch and Walsh 1998). In the age-structured model that follows, we define a bivariate normal distribution with a known variance-covariance structure as our starting point and iterate this forward (models B-D in sec. A1.2-A1.4). We consider a simple case: a monocarpic biennial life cycle, where individuals in their first year of life do not reproduce and all age 2 individuals die after reproduction. As with our model for a species with nonoverlapping generations, we assume clonal inheritance, where survival from age 1 to age 2 is specified as with expected survival to age 2 being highest for larger values of A1. Although A2 is not under direct selection, its distribution is modified by its covariance with A1. A2, the genotype at age 2, determines expected reproduction, Although A1 does not directly influence reproduction, there is an association between it and reproduction via its covariance with A2. All age 2 individuals die following reproduction in this model, although it is possible to extend our approach to any arbitrary number of ages. The evolutionary dynamics that particular parameterizations of the fitness functions S(A1, 1, t) and R(A2, 2, t) generate are dependent on (1) the initial covariance between the characters and (2) the fitness functions (models B-D in sec. A1.2-A1.4). Many parameterizations and initial covariances are likely to generate evolutionary dynamics that may be biologically unrealistic. We demonstrate this with three contrasting parameterizations, considering size as our trait ( fig. 2e-2g). In the first example ( fig. 2e; model B in sec. A1.2), the two characters positively covary and experience selection in the same direction. Over the course of the simulation the average developmental trajectory has evolved, with A1 evolving to be 1.76 times larger and A2 evolving to be 1.52 times larger. For a trait such as body size, such a proportional change at different ages may be appropriate. In examples in figure 2f and 2g (models C, D in sec. A1.3, A1.4), the bivariate character evolves in contrasting ways. In model F, A2 evolves much faster than A1, while in model G, A1 evolves to be larger, while A2 evolves to be smaller. These simulations demonstrate that only a constrained set of fit- The colored image plots in e-g represent Gaussian development functions D(Z 0 jZ, t) fitted to the bivariate distributions of A at the beginning and end of the simulation. Evolution of the bivariate character has resulted in different parameterizations of these phenomenological functions. The lighter the shading, the greater the probability of a transition from character value Z at age 1 to Z 0 at age 2.
ness functions and genetic covariances will give biologically realistic evolutionary trajectories for the size-related traits that biologists often study.
We now return to a univariate model and examine the clonality assumption. How can the clonality assumption be relaxed, and what are the consequences? In sexually repro-ducing species, offspring inherit a mix of their parent's genomes. However, genetic segregation means that full siblings do not have the same genotype. When additivity is assumed, the breeding value of offspring is expected to be midway between parental breeding values. However, to obtain the distribution of offspring genotypes, the contribu- Dynamics of inheritance (model E in sec. A1.5 in the appendix, available online). The dynamics of population growth rate R 0 (a) and the mean (b) and variance (c) of A vary between models with clonal inheritance (black line), the convolution in equation (15) (red line), and the Gaussian inheritance function in equation (16) (blue line). Dynamics predicted from the convolution and the Gaussian inheritance function are indistinguishable in this model. d, Temporal dynamics of the clonal model versus the other models. The initial distribution at t p 1 is Gaussian. After 100 generations the character distributions predicted by the clonal and sexual models have diverged only slightly. The infinitesimal model of quantitative genetics assumes that the dynamics of alleles can be inferred from the dynamics of phenotypes. Under this assumption, selection alters genotype and allele frequencies, while inheritance does not (e). In contrast, when dominance variance operates, both selection and inheritance alter genotype frequency, while neither alters allele frequencies ( f ). For a Gaussian distributed character, dominance variance acts as an offset, reducing the intercept of a Gaussian inheritance function ( g). tion of genetic segregation to variation among offspring needs to be taken into account. In two sex models, three steps are required to generate the distribution of offspring genotypes or breeding values, given parental values. First, a distribution of mating pairs needs to be constructed. Second, the distribution of midpoint parental genotypes or breeding values given the distribution of mating pairs needs to be constructed. Third, segregation variance needs to be added to the distribution (Feldman and Cavalli-Sforza 1979;Felsenstein 1981;Turelli and Barton 1994). The mating system and the segregation variance are related: when mating is assortative with respect to genotype, the segregation variance is small, and siblings closely resemble one another and their parents. In contrast, when mating is disassortative with respect to genotype, siblings can differ markedly from one another, and the segregation variance is large.
Expressions have been derived for the segregation variance for the infinitesimal model, where it is assumed that traits are determined by a very large number of unlinked loci of small additive effects and mating is random (Fisher 1930). The infinitesimal model is assumed in most empirical quantitative genetic analyses (Kruuk et al. 2008) and in our initial model. For random mating where both sexes have identical demographies, the distribution of offspring breeding values given parental breeding values is (Barfield et al. 2011) where * represents convolution and f(A, j 2 ) p (1= ffiffiffiffiffiffiffiffiffi ffi 2pj 2 p ) exp(2A 2 =j 2 ) is a Gaussian function with mean 0 and variance j 2 representing the segregation variance.
If males and females have different demographies, then they will have different distributions of genetic values after selection; we represent these as N r M (A, t) and N r F (A, t), respectively. In this case, equation (14) is replaced by where j 2 r ðMÞ (A, t) and j 2 r ðFÞ (A, t) are variances of the postrecruitment selection breeding value of males and females, respectively. We do not superscript the r's with j 2 to avoid a notation making it appear that j is raised to some quantity 2r.
The first two terms on the right-hand side of equation (15) generates the distribution of expected parental midpoint values; it ensures that the mean breeding value among off-spring is midway between the two parental breeding values. However, because the parental distributions are halved, the variance of this distribution is half that of the parental distributions. The third term on the right-hand side of equation (15) adds the segregation variance. For random mating, the variance is assumed to be normally distributed, with a mean of 0 and a variance of half the additive genetic variance among the entire population when the population is at linkage equilibrium (Felsenstein 1981). We approximate this variance as half the additive genetic variance in the parental distribution (Feldman and Cavalli-Sforza 1979). This approach has already been incorporated into IPMs (Barfield et al. 2011;Childs et al. 2016).
We now run two simulations ( fig. 3a-3d) to examine differences in the predictions of clonal and sexual models. The first model assumes clonal inheritance and the second the convolution in equation (15), with both models assuming a linear function R(Z, t) (model E in sec. A1.5). The two models predict slightly divergent dynamics. The reason for this is that equation (15) results in the skew and kurtosis in N R (A, t) is reduced at each time step in the sexual model compared with in the clonal model. If selection is exponential (and the starting distribution proportional to a Gaussian distribution), then there will be no difference between the two approaches. This is because a normal distribution multiplied by an exponential fitness function results in a normal distribution with an unchanged variance (Diaconis and Ylvisaker 1979). These results suggest that insights from clonal models will approximate those from sexual models reasonably well, at least when males and females have similar demography.
In the remainder of this section we explore simple approximations of the models of breeding values described above. Some authors have queried the use of equation (3) as an approximation in IPMs to the inheritance convolution in equation (15) used in models of sexually reproducing species (Chevin et al. 2010;Janeiro et al. 2017). However, being able to construct inheritance functions for A that are of the form of equation (3) would be useful because it would permit methods developed for two-sex phenotypic IPMs to be applied to evolutionarily explicit IPMs (e.g., Schindler et al. 2015). Given that Gaussian approximations frequently perform well in models of evolution (Turelli and Barton 1994), we hypothesize that Gaussian inheritance functions may perform well in evolutionarily explicit IPMs. We consequently constructed a Gaussian inheritance function and compared results with those obtained from the convolution.
Equation (15) results in the mean and variance of the parental and offspring breeding value being the same. We can approximate this by ensuring that the function m H (A, t) passes through the coordinate x p E(N R (A, t)), y p E (N R (A, t)) and that the variance V H (A, t) p j 2 (N R (A, t)). When both sexes have the same demography, we can write where E and j 2 represent expectations and variances, respectively, and h represents the degree of assortative mating. When h p 1, mating is entirely assortative; when v p 0:5, mating is random; and when h p 0, mating is completely disassortative. An equation for the case when males and females have different demographies is provided in section A3. The approximation in equation (16) will increase in accuracy as the distribution of midpoint parental breeding values becomes more Gaussian.
When we compared predictions from equations (15) and (16) with h p 0:5 using the same model used to compare clonal and sexual life histories, results were indistinguishable ( fig. 3a-3d). This reveals that for linear selection, Gaussian inheritance functions for A perform remarkably well.
None of our models to date include any form of mutation. We have not incorporated mutation into our models because we are simulating responses to environmental change over a few tens to hundreds of generations (figs. 1-3), and over that time period mutation is unlikely to play a major role in adaptation. However, for simulations over longer time periods, we can incorporate mutation into our models by slightly increasing the size of the segregation variance (e.g., Lynch and Walsh 1998). This will have the effect of increasing the additive genetic variance, partly countering any loss of genetic variance due to selection.
Our approximation can be used to examine the dynamical contributions of nonadditive genetic processes to population responses to environmental change in a phenomenological manner. Fisher (1930) demonstrated that dominance variance can be treated as an offset, and in our models this would lower the intercept of the function m H (G, t) in equation (16). A consequence of this is that the mean of the offspring genotype is no longer equal to the mean of parental genotype and the dynamics of genotypes no longer exactly match the dynamics of alleles. We demonstrate this with a single-locus two-allele model. When the effects of alleles are additive, the dynamics of the genotype captures the dynamics of alleles ( fig. 3e). In contrast, when the heterozygote has higher fitness, allele frequencies do not change once the equilibrium is achieved. However, selection and inheritance alter genotype frequencies ( fig. 3f ). This effect of dominance variance can be phenomenologically captured within an IPM by setting the intercept of the inheritance function for the genetic component of the phenotype to be less than E R (N R A, t)=2; this imposes an offset that can reverse gains made by selection ( fig. 3g). Because this offset is negative when dominance variance is operating, dominance variance will slow rates of evolutionary change. We could easily phenomenologically explore how a particular value of this offset impacts predicted dynamics; however, further work is required to relate different levels of dominance variance to specific values of the offset in our models. Having shown how IPMs can be formulated to project forward the dynamics of the genetic component of the phenotype, we now turn our attention to the dynamics of the environmental component of the phenotype.

Plasticity: The Dynamics of E
Plasticity is determined by the dynamics of E and in particular in how E is influenced by the ecological environment v. To capture plasticity in IPMs we need to model the probability of transition from E at time t to E 0 at time t 1 1 as a function of the environment v. For most plastic traits we have a poor mechanistic understanding of development and inheritance patterns, and for that reason we use the Gaussian probability density function in equation (3).
In quantitative genetics it is often assumed that the mean of E(E, t) p 0, and any individual departures are purely random (Falconer 1960). In equation (3)  Gaussian transition functions (eq. [3]) can be formulated to predictably modify moments of the distribution of E from time t to time t 1 1. For example, careful choice of intercepts and slopes of m D (E, t), m H (E, t), V D (E, t), and V H (E, t) can be used to predictably grow or shrink the variance of E via either development or inheritance (sec. A4). In addition, specific biological processes can be easily incorporated into the dynamics of E: if the slopes m D E ( 0 or m H E ( 0, then there will be temporal autocorrelation in the value of E among individuals and between parents and their offspring. For example, if m D E 1 0, then individuals with a relatively large value of E at time t will be expected to have a relatively large value of E 0 at time t 1 1. This property of development functions is useful because it allows some memory of E across ages: if an individual has benefited from a particularly good set of circumstances at one age, any phenotypic consequences can persist to older ages. In a similar vein, if m H E 1 0, then a parent with a relatively large E at time t will produce offspring with relatively large E's at time t 1 1, a form of parental environmental effect (Nussey et al. 2007).
Different formulations of m H (:::) and m D (:::) can be used to capture a variety of different forms of plasticity (table 2). When v is incorporated as an additive effect, it acts to shift the intercept of these functions as t changes. This means that the environment influences all values of A in the same manner. If Z p A 1 E, then Z changes as a function of how v influences E if A remains constant. A remains constant as it does not vary within individuals as they age or if A 0 in offspring is the same as A in parents.
Interactions between E, A, and v are listed in table 2. Each form describes a different type of reaction norm (Gavrilets and Scheiner 1993). These forms allow E to develop among individuals (phenotypic plasticity) or be inherited (epigenetic inheritance) as a function of an individual's breeding value A and the environment v as well as the value of E at time t.
Plasticity can be either adaptive or nonadaptive (Ghalambor et al. 2015), and both forms can be captured in our models. Adaptive plasticity enables populations to rapidly respond to an environmental change. For example, if environmental change reduces population size, then adaptive plasticity would result in a change to the mean of the phenotype via either phenotypic plasticity (the development function) or epigenetic inheritance (the inheritance function) that leads to an increase in survival or recruitment rates. In contrast, nonadaptive plasticity does the opposite, potentially exacerbating the detrimental effects of environmental change.
We demonstrate this with an example of a simple IPM of a species with nonoverlapping generations: The model contains no genetic variation, and the phenotype is determined by the density at the time the offspring is born. This means that we can remove A from the model. We assume a linear fitness function and a Gaussian inheritance function, Next, we assume that the phenotypic trait is positively associated with expected recruitment such that R E 1 0. We also assume that the environmental driver is positively associated with expected recruitment such that as v increases in value and fitness increases (R v 1 0). This means that the population growth rate (in a density-independent model) or population size (in a density-dependent model) also increases with v. Now assume that a negative environmental perturbation decreases v such that fitness decreases. For adaptive plasticity to counter this, the effect of the decrease in v on epigenetic inheritance must increase the expected value of E. In our simple model, this can occur only if m H v ! 0. Then, as v declines, m H v v becomes less, and the value of m H I 1 m H v v becomes larger, increasing the mean of E and fitness. In general, in additive linear models such as this, if R E and m H v take opposing signs, then plasticity will be adaptive. We develop three density-dependent models of a phenotype in a species with nonoverlapping generations. In all models we define the fitness function to be R(E, t) p R I 1 R E E 1 R n(t) n(t), where n(t) p Ð N(E, t)dE and R n(t) ! 0. In each model we define m H (E, t) p m H I 1 m H E E 1 m H n(t) n(t). We set in model F m H n(t) p 0; in model G m H n(t) ! 0; and in model H m H n(t) 1 0 (sec. A1.1). The first model (model F) does not include plasticity (m H n(t) p 0), the second (model G) captures adaptive plasticity (m H n(t) ! 0 and R E 1 0), and the third (model H) captures nonadaptive plasticity (m H n(t) 1 0 and R E 1 0). All three models include temporal autocorrelation in the environmental component of the phenotype (sometimes referred to as phenotypic carryover) when m H E 1 0 (table 2). Because the models are not age structured and do not include development, plasticity operates via epigenetic inheritance (e.g., maternal environmental effects). The same logic can be extended to the development function in age-structured populations. In our examples, parameterizations are chosen so that all models converge to the same value of carrying capacity K. Once all three models have converged, we initially impose a one-off perturbation. Model G regains the equilibrium first, followed by model F and then model H ( fig. 4a), showing that adaptive plasticity allows the population to recover from a one-off environmental perturbation much faster than when there is no plasticity or plasticity is nonadaptive. Nonadaptivity plasticity significantly slows the rate at which the population can recover from a perturbation, with the initial population size before perturbation reattained only after 80 generations.
Adaptive and nonadaptive plasticity also impact the way populations respond to permanent environmental change. We demonstrate this by running the same models F, G, and H, except now we impose a constant change in fitness by permanently changing the intercept of the fitness function R I . When we do this, the three models attain different equilibria population sizes ( fig. 4b) and different mean phenotypes ( fig. 4c). Model G achieves a larger population size than the two other models. This buffering of the population against environmental change happens because adaptive phenotypic plasticity results in a change in the mean phenotype ( fig. 2c) that increases the expected recruitment rate and asymptotic population size ( fig. 2b). In contrast, nonadaptive plasticity exacerbates the consequences via a change in the mean phenotype that decreases fitness.
In contrast to our example models in "Adaptive Evolution: The Dynamics of A," the IPMs we have developed in this section-and, indeed, all nongenetic IPMs so far published-achieve an asymptotic population growth rate or equilibrium population size and a stable population structure. These IPMs have monotonically increasing or decreasing fitness functions: an increase in the character results in an increase in expected fitness. A consequence of this is that in these models the recruitment function acts to alter the location of the character distribution and often also alter its shape (Wallace et al. 2013). This is reflected in that the means (and often other moments) differ between the distributions of the phenotype before and after selection. In models at equilibrium with monotonic fitness functions, the inheritance function must reverse the locational and shape changes caused by the fitness function. This is because at equilibrium the moments of the phenotype distribution at times t and t 1 1 must be equal.
In models of species with nonoverlapping generations at equilibrium, such as those above, the inheritance function for E must exactly reverse the changes to the character distribution generated by the fitness function. This requires moments of parental and offspring characters to differ from one another if N R (E, t) 2 N(E, t) ( 0. When there is a correlation between parental and offspring traits in the inheritance function for E, as in our models, the intercept of the inheritance function must take a value such that offspring characters are smaller than their parents' were at the same age (Coulson and Tuljapurkar 2008).
IPMs for species with overlapping generations include development functions D(E 0 jE, a, t). These functions can alter the size (population size) and shape of the distribution of E as individuals age. When generations are overlapping, and at equilibrium, changes to the location of the character distribution via survival, recruitment, and development are all exactly countered by the inheritance functions H(X 0 jX , a, t). Coulson and Tuljapurkar (2008) showed that in red deer, age-specific effects meant that young and old parents were incapable of producing offspring that had the same body weight as they did at birth. This process reversed the effects of viability selection removing small individuals from the population early in life. The same process was observed in marmots (Ozgul et al. 2010) and Soay sheep (Ozgul et al. 2009) and may be general for body size in mammals.
We have now developed IPMs for (1) A where we assumed all individuals had the same constant E and (2) E where we assumed all individuals had the same A. We have shown how IPMs can capture a wide range of biological processes, including adaptive and nonadaptive plasticity and correlated characters, and the circumstances when equilibria are achieved. We now link together these advances into models of the joint dynamics of the bivariate distribution N(A, E, t).

Models for the Phenotype Consisting of Genetic and Environmental Components
Here we construct models where the character can be determined by a mixture of the genetic and environmental components. These models allow us to explore how adaptive evolution is influenced by plasticity. We first develop a dynamic univariate version of the breeder's equation (Falconer 1960) for a species with nonoverlapping generations in a constant environment. In this case, the environmental component of the phenotype is assumed to be a consequence of developmental noise: individuals achieve their genetic potential, plus or minus a departure. At each generation within each breeding value, the distribution of the environmental component of the phenotype is assumed to be Gaussian with a mean of 0 and a constant variance (model I in sec. A1.9).
Our initial conditions are a bivariate Gaussian distribution of A and E, which we iterate forward for 300 time steps. Over time, the mean of the genetic component of the phenotype increases. In contrast, the mean of the environmental component is constant. The population grows hyperexponentially ( fig. 5a), the mean of the phenotype increases in value because of evolution ( fig. 5a, 5d ), and the additive genetic variance is slowly eroded (fig. A2). Because the additive genetic variance is eroded while the phenotypic variance remains constant, the heritability declines over time ( fig. A2).
Our second model (model J in sec. A1.10) has a negative density-dependent term in the fitness function. The phenotype evolves faster in this model than in our densityindependent model ( fig. 5b). Population size grows nearly linearly in this model ( fig. 5d ), although the rate of increase does slow slightly each generation as genetic variation is eroded. The difference between the hyperexponential (density-independent model) and nearly linear increases (density-dependent model) in population size explains the difference in the rates of evolution. This is because the selection differential that determines the rate of evolution (an emergent property from our model [Wallace et al. 2013]) has the population growth rate in its denominator. The population growth rate is smaller in the density-dependent model (just above unity) than in our density-independent one (it increases with time), and this leads to an increase in the strength of selection and the rate of evolution (see also Pelletier and Coulson 2012). A consequence of this is that the additive genetic variation and heritability tend toward 0 faster in the density-dependent model than in the densityindependent one ( fig. A2).
In our third model (model K in sec. A1.11), negative density dependence is included in the inheritance function for the environmental component of the phenotype as well as in the fitness function. This captures adaptive phenotypic plasticity. This results in a negative change in the mean of the environmental component of the phenotype with time ( fig. 5c). This decrease is reflected in a change in the mean of the phenotype itself. Adaptive phenotypic plasticity leads to a decline in the population growth rate, which results in a slight increase in the rate of evolution compared with the density-dependent model with no plasticity. However, the effect is not large and is only just distinguishable when comparing figure 5b, 5c.
In our final models (models L-N in sec. A1.12-A1.14) we examine how a one-off perturbation influences the mean of the phenotype, its components, and the population growth rate ( fig. 5g-5l ) when there is no plasticity, adaptive plasticity, and nonadaptive plasticity. We set the variance in the genetic and environmental component of the phenotype to be equal, giving an initial heritability of h 2 p 0:5. In each model we allow the population to achieve the same equilibrium population size in the absence of selection (R Z p 0). We then impose a one-off mortality event when 99% of individuals above the mean of the phenotype are killed off. At this point we also impose selection (R Z p 0:1). In all three models the mortality event results in a small change in the mean value of the phenotype (for an explanation, see sec. A5; fig. 5g-5i, red lines) but a halving of population size ( fig. 5j-5l). Adaptive plasticity results in the environmental component of the phenotype returning to its preperturbation value very quickly ( fig. 5g-5i, blue lines). In contrast, although the perturbation causes a modest change in the mean of the genetic component of the phenotype, it takes 110 generations for evolution to reverse the change ( fig. 5g-5i, black lines). This demonstrates that a strong selective effect can leave a large population dynamic impact but leave only a small initial signature in the phenotype, even when the trait is highly heritable.
Over the longer term, the dynamics of the components of the phenotype, the phenotype itself, and the population dynamics all depend on whether plasticity is adaptive or nonadaptive. Adaptive plasticity allows the population size to initially recover from the perturbation more quickly than when plasticity is absent or nonadaptive ( fig. 5j-5l). However, over a longer time period, nonadaptive plasticity results in the population achieving a larger size than when plasticity is absent or adaptive. These differences in population growth rate impact rates of evolution: immediately following the perturbation, the rate of evolution is greatest when plasticity is nonadaptive. However, the rate of evo-lution then increases when plasticity is adaptive (figs. A2, A3). As with our previous models, the effects of adaptive and nonadaptive plasticity on rates of evolution are relatively small, but our results demonstrate how the two processes can interact.

Signatures of Evolution in Models That Are Not Evolutionarily Explicit
The models in the previous section are quite complex. Do we always need to construct such evolutionarily explicit IPMs The dynamics of the phenotype (red lines) and its genetic (black lines) and environmental (blue lines) components (a-c, g-i) and the dynamics of the population (d-f, j-l). In the first model (a, d), both fitness and inheritance of the environmental component of the phenotype are independent of density (model I in sec. A1.9 in the appendix, available online). In the second model (b, e), fitness is negatively density dependent, and inheritance of the environmental component of the phenotype is density independent (model J in sec. A1.10). In the third model, both fitness and inheritance of the environmental component of the phenotype are negatively density dependent (model K in sec. A1.11). The treatment of plasticity can dramatically influence the dynamics of the phenotype and population size (models L-N in sec. A1.12-A1.14). Adaptive phenotypic plasticity (h, k) leads to the population size and phenotype recovering from a perturbation much faster than nonadaptive plasticity (i-l). The absence of a plastic response ( g, j) results in the population recovering from a perturbation at an intermediate rate between cases where adaptive and nonadaptive plasticity are operating.
to predict population responses to environmental change, or can we rely on simpler, phenotypic IPMs? There are two reasons why it may be preferable to not construct evolutionarily explicit models. First, evolutionarily explicit IPMs are more complicated to construct than those that do not include genotypes or breeding values. Second, when data are unavailable to explicitly include breeding values into models (Traill et al. 2014b), the effects of evolution on predicted dynamics can still be explored by examining the consequences of perturbing parameter values (Traill et al. 2014a). When evolution occurs within a system we would expect parameters in phenomenological inheritance and development functions that are fitted to data to change with time.
We can see this in figure 2e-2g. In these age-structured evolutionarily explicit models, the bivariate breeding value distribution (black contours) changes location as evolution occurs. We have fitted Gaussian development functions to these bivariate distributions at the beginning of each simulation and at the end (colored image plots). The parameters that determine these developments functions have clearly changed as the location of the functions have changed. A similar process occurs for inheritance functions (not shown).
Numerous authors have previously noted this phenomenon in models of evolution. For example, in population genetic (Charlesworth 1994) and ecoevolutionary models (Yoshida et al. 2003;Coulson et al. 2011) when genotype frequencies change with time, macroscopic, populationlevel quantities such as mean survival and recruitment also change; in adaptive dynamic models, as one strategy invades another, population-level parameters inevitably change with strategy frequency over time (Metz et al. 1996); in quantitative genetic predator-prey models, population-level parameters of both predators and prey vary over time, leading to persistence of the interaction (Doebeli 1997); and in evolutionarily explicit IPMs, parameters in inheritance functions have been shown to change with time as evolution progresses . These insights are useful because if evolution is occurring within a system, then temporal trends in statistical estimates of model parameters would be expected; in other words, the effect of time-either additively or in an interaction with other parameters-would be expected in m H (Z, t), m H (Z, a, t), or m D (Z, t). If marked temporal trends are observed in parameters in development and inheritance functions that cannot be attributed to a changing environmental driver, then evolutionarily explicit IPMs may be required.
What about parameters in fitness functions S(Z, t) and R(Z, t)? Can any inferences from temporal trends in these parameters be made? In our approach, evolution of a focal trait would not be expected to alter statistical estimates of the fitness functions. In our models, evolution simply moves the location and shape of the phenotype distribution, but not its association with survival or recruitment.
We have identified one circumstance where evolution will leave a signature in the dynamics of fitness function parameters. Parameters in these functions can evolve in the presence of a genetically unmeasured correlated character that is also evolving. To demonstrate this we construct a model of a bivariate character and examine the dynamics it predicts before exploring the consequences of failing to measure one of the characters.
We assume clonal inheritance such that dynamics of the characters are solely determined by a bivariate fitness function, The dynamics this model predicts depend on the initial covariance between the two characters in a similar way to our age-structured model (eq. [11]). In our first example the two characters negatively covary, while in the second they positively covary (for model parameterizations, see sec. A1.1). The initial negative covariation allows rapid evolution, with population growth ( fig. 6a), the mean of the characters ( fig. 6b), their variances ( fig. 6c), and the covariance between them ( fig. 6d) evolving relatively quickly. In contrast, when the two characters positively covary, evolution is much slower, with the character means, variances, and covariance changing much more slowly, even though the fitness functions are identical in each model ( fig. 6e-6h). We now construct a fitness function for A1 when A2 is not measured. We start by defining mean fitness, an observable, as E(R:t) p E(R (A, t)). The slopeR A1,t is given byR The intercept can be calculated in the usual manner by estimating the means of fitness and A1, giving Equation (20) is what would be estimated from data if A2 were not measured and included in analyses (Söderström and Stoica 2002;Kendall 2015). It will correctly describe the consequences of selection on A1 even though A2 could be correlated with it. This is because the unmeasured correlated character impacts fitness whether it is measured or not and consequently impacts the association between the focal character and fitness in its absence (Lande and Arnold 1983). However, the fitness function cannot provide accurate predictions over multiple generations when it is assumed that the fitness function is constant.
Over multiple generations the existence of unmeasured correlated characters will alter parameters in the fitness function in equation (20) if selection alters genetic variances and covariances of measured and unmeasured correlated characters ( fig. 6i, 6j). This is becauseR I,t andR A1,t are both functions of the covariance between the two characters (eqq. [18]- [20]). If selection alters this covariance, parametersR I,t and R A1,t will evolve with time. It is also why we use the subscript t forR I,t andR A1,t . Evidence of correlated characters under selection can consequently be inferred if parameters in fitness functions are observed to change with time in a system in the absence of a changing environmental driver. Note that a nonstationary unmeasured environmental driver could also generate trends in parameter values in fitness functions in phenomenological IPMs.

Discussion
In this article, we develop an approach that allows prediction of how populations respond to environmental change Character 2 Figure 6: Dynamics of bivariate characters and evolution of fitness functions in the presence of an unmeasured, genetically correlated character (models P, Q in sec. A1.15 in the appendix, available online). We construct a simple model with clonal inheritance of two correlated characters that both influence fitness. We explore two initial starting conditions that differ only in their genetic covariance (models P, Q in sec. A1.15). In a-d, the covariance accelerates the rate of evolution compared with e-h. i, j, Dynamics of the fitness function for each character when the other character is not measured. Regardless of the covariance between characters, the fitness functions evolve during the course of 120 time step simulations.
via adaptive evolution and plasticity. We do this by incorporating insights from evolutionary genetics into data-driven structured population models. Our approach is to split the phenotype into its genetic and environmental components and to model the dynamics of the genetic component with functions based on understanding of the mechanisms of inheritance. In contrast, the dynamics of the environmental component of the phenotype are modeled with phenomenological functions that can be identified from the analysis of data. Our approach is appropriate for sexually reproducing or clonal species with either overlapping or nonoverlapping generations.

Evolutionarily Explicit Structured Models
IPMs are now a widely used tool in ecology and evolution because of their versatility and the ease with which they can be parameterized (Merow et al. 2014). All key statistics routinely estimated in population ecology, quantitative genetics, population genetics, and life history describe some aspect of a character distribution or its dynamics . IPMs are so versatile because they describe the dynamics of these distributions. Characterization of the determinants of these statistics gained via sensitivity or elasticity analysis of models has provided insight into how ecological and evolutionary quantities that interest biologists are linked (Coulson et al. 2011). Although this logic was developed several years ago, there has recently been criticism that IPMs cannot be used to track the dynamics of multivariate breeding values expressed at different ages (Chevin 2015;Janeiro et al. 2017). Our article addresses this criticism head-on; we show how IPMs can be formulated to capture a mechanistic understanding of inheritance and development. In demonstrating this we develop a general modeling approach to capture population responses to environmental change. Having done this, we are now in a position to construct IPMs of quantitative characters and examine how perturbing the environment will influence not only the dynamics of the phenotype and its genetic and environmental components but also the life history (Steiner et al. 2012(Steiner et al. , 2014 and population dynamics (Easterling et al. 2000). The work we present here adds to a growing literature that explicitly incorporates evolution into structured models and IPMs in particular. Within the population genetics paradigm, Charlesworth (1994) developed structured models with a one-to-one map between genotype and phenotype in age-structured populations. Building on this work, Coulson et al. (2011) showed how simple genetic architectures can be incorporated into IPMs, developing a model to explore how evolution at a single locus would occur simultaneously with phenotypic change of discrete and continuous characters, life history, and population dynamics.
Working in the quantitative genetic paradigm, Lande (1982) derived age-structured models that tracked the dynamics of the mean of the additive genetic component of the phenotype (E(A) in our notation) and the mean of the phenotype itself (E(Z)). He assumed a constant genetic variance-covariance matrix and consequently weak selection and normally distributed character values-assumptions we relax. Barfield et al. (2011) extended Lande (1982's approach to track the dynamics of the entire character distribution and to stage-structured populations. In doing so, they developed a general, flexible approach to track the entire distributions of A and Z. Childs et al. (2016) extended this approach to two sexes. Because A is inherited with mechanistic rules that are not impacted by the environment, while inheritance and development of E are plastic and can be impacted by the ecological environment (Falconer 1960), it is difficult to incorporate the effects of the environment on the dynamics of the phenotype by focusing on A and Z as Lande (1982), Barfield et al. (2011) andChilds et al. (2016) have done. In contrast, our approach (which otherwise has a similar logic to Barfield et al. [2011] and Childs et al. [2016]) tracks the dynamics of E and A (or G-the full genotypic value, including nonadditive componentsif desired), making incorporation of environmental drivers that influence inheritance and development of E more straightforward. We show that it is possible to have selection operating on the phenotype while incorporating modern understanding of genetic inheritance into the dynamics of the genetic component of the phenotype and phenomenological insight into the role of the ecological environment on the dynamics of the environmental component of the phenotype. By doing this, we show how population responses to environmental change via adaptive evolution, phenotypic plasticity, and epigenetic inheritance can be simultaneously explored. This opens up the way to provide novel insights into the circumstances when each process is expected to contribute to population responses to environmental change.

Population Responses to Environmental Change
Unlike previous evolutionarily explicit IPMs (Barfield et al. 2011;Childs et al. 2016;Rees and Ellner 2016), our approach requires explicit consideration of the inheritance and development of E, the environmental component of the phenotype. This allows our models to capture a range of plastic responses to environmental change along with adaptive ones. What do our findings say about the contributions of plasticity, evolution, and their interaction to population responses to environmental change?
Detrimental environmental change often causes a decline in population size. When there is an association between a phenotypic trait and survival and recruitment rates, phenotypic change can lead to increased survival and re-cruitment rates (Ozgul et al. 2010) and consequently an increase in population growth rate and size. Two processes can lead to phenotypic change: plasticity and adaptive evolution. There has been considerable discussion about the relative roles of each in allowing populations to respond to change (e.g., Chevin et al. 2010;Bonduriansky et al. 2012).
Genotypes and breeding values remain fixed within individuals throughout life, which means that differential survival and recruitment rates are the processes that alter these distributions and underpin evolution. The strength of differential survival and recruitment can be impacted by environmental variation generating fluctuating selection (Lande 2007). Environmental variation does not influence genetic inheritance: once mating pairs are formed, inheritance of breeding values, A, does not alter the mean or variance of breeding value distributions (Fisher 1930). In contrast, distributions of the environmental component of the phenotype can be altered via survival, recruitment, development, and inheritance, with each process potentially impacted by environmental variation (Reed et al. 2010). Given these differences between the dynamics of A and E, plasticity can lead to more rapid change than evolution in our models (e.g., fig. 5). This is because more biological processes can directly alter the distribution of plastic characters than can impact distributions of breeding values. These results are consistent with those of other authors, including Lande (2009) and Chevin et al. (2010), who also concluded that plastic change should be faster than evolutionary change. But how quickly will evolution alter phenotypic trait distributions?
Our results on the speed of evolution suggest that claims of detectable rapid evolution in quantitative phenotypes is likely to take a few tens of generations. For example, environmental change increases mortality, leading to a decline in population size, but for mortality selection to lead to evolutionary change over the course of a generation, a large proportion of the population needs to be selectively removed, and the phenotype needs to be highly heritable. This is seen in our model results ( fig. 5g-5i) and with a simple numerical example: when all individuals above the mean of a normally distributed character are removed from the population and the trait has a heritable of h 2 p 0:5, population size halves in a single time step, but the mean of the character will only shift from the 50th percentile to the 37.5th percentile. For a standard normal distribution with a mean of 0 and a standard deviation of unity, this means that the mean would shift by only 0.319, that is, less than one-third of a standard deviation. In reality, mortality selection resulting from environmental change will likely result in a change to the mean of the distribution that is only a fraction of a standard deviation compared with our example. Given this, reports of rapid evolution due to environmental change increasing mortality selection over a small number of generations (e.g., Coltman et al. 2003) should be treated with caution. It is much more likely that change is a consequence of phenotypic plasticity. Over multiple generations, recruitment selection can also contribute to evolutionary change, and our approach allows the role of this to be investigated. However, unless reproduction is restricted to individuals with extreme phenotypic trait values in both sexes, it seems unlikely that evolution can generate statistically demonstrable evolutionary change over a small number of generations ). This is not to say that evolution is not important over longer timescales. Over tens of generations evolution can shift phenotypic trait means to a greater extent than phenotypic plasticity ( fig. 5g-5i, blue vs. black lines).
In order for plasticity to allow populations to rapidly respond to environmental change, a large proportion of individuals within the population must exhibit the same plastic response. A good example of such a dynamic is for sizerelated traits that are determined by resource availability, particularly when scramble competition is operating. When resources become limiting, all individuals will be unable to develop as rapidly as when resources are more common. A consequence of this is that individuals that developed in cohorts when resource were sparse will exhibit smaller body sizes compared with individuals in those cohorts that developed when resources were more abundant. We can capture this form of plasticity in our framework with an additive effect of density in the inheritance or development function for E (e.g., fig. 4). In contrast, when contest competition operates, larger individuals would acquire more resources than those that are smaller and would develop faster. We can capture this in our models with interactions between density, E, and A in either the inheritance or development functions for E.
This discussion demonstrates how our approach can be used to capture different forms of plasticity. However, for plasticity to help populations respond to environmental change, it must be adaptive: plasticity must change the mean trait value in a way that increases fitness (Ghalambor et al. 2007). We demonstrate that for additive, linear models, adaptive and nonadaptive plasticity can be specified by altering the sign for the effect of the environment in the function specifying the mean dynamics of the inheritance or development functions ( fig. 4). When interactions are included in these functions, specifying general rules for whether plasticity is adaptive or nonadaptive will likely be more challenging. However, our approach provides a way in which to investigate when plasticity is adaptive or nonadaptive and how different types of plasticity will influence population responses to environmental change.
Our results also show how plasticity can influence evolutionary rates. Plasticity-operating via development and inheritance functions for the environmental component of the phenotype-alters the distribution of the phenotype, and this can alter the strength of selection, which can then influence the dynamics of the genetic component of the phenotype (evolution). The effects of plasticity on selection and evolution can be surprisingly complex. We examined only the evolutionary consequences of plasticity following an environmental shock that influenced all individuals in the same way, but even in this simple case we found that adaptive plasticity initially slowed the rate of evolution compared with nonadaptive plasticity before increasing it ( fig. 5; appendix). In general, in order to understand how plasticity will influence selection, it is necessary to understand how it influences both the numerator and the denominator of the selection differential that underpins evolution (Pelletier and Coulson 2012). The numerator is the covariance between the phenotype and absolute fitness (Falconer 1960), and the denominator is mean fitness. In our models of species with nonoverlapping generations, this is mean recruitment-the population growth rate (Fisher 1930). Selection is linear in our models, where plasticity influences all individuals in the same way via an additive effect of density on inheritance of the environmental component of the phenotype (fig. 5), and this means that plasticity influences the population growth rate rather than the numerator of the selection differential. A consequence of this is that it is differences in the population growth rate that generates the differences in evolutionary rates between models when plasticity is adaptive and nonadaptive. In more complex cases when plasticity influences the covariance between the phenotype and fitness via genotype-environment interactions within a generation, to understand how selection influences evolution it is necessary to understand not only how plasticity influences mean fitness but also how it generates differences between the covariance between the genetic component of the phenotype and fitness and the covariance between the phenotype itself and fitness. Because the components of the selection differential can be calculated from IPMs Wallace et al. 2013), the approach we develop here provides a flexible way to examine how different types of plasticity can influence evolution following environmental change.
We have not considered bet hedging in this article. Bet hedging is another form of plasticity that can influence the way populations respond to environmental change, and it can be incorporated into IPMs . Deterministic IPMs incorporate probabilistic transitions when V H (E 0 jE, A, t) 1 0 and V D (E 0 jE, A, t) 1 0. These probabilities do not vary from one time step to the next. In stochastic models these functions can include terms for an environmental driver v, such that the variation in trajectories changes with the environment. In evolutionarily explicit models, the variance in transition rates among different values of E can be made to depend on v, A, and their interaction (if desired). This means that individuals with spe-cific values of A can produce offspring with more variable values of E (and consequently Z) in particular environments than individuals with other values of A. In this article we focused on the incorporation of v into m H (E 0 jE, A, v, t) and m D (E 0 jE, A, v, t), but responses to environmental change could also be incorporated into functions for the standard deviation that we use to construct our kernels. In order to explore how the various forms of plasticity influence rates of evolution for real systems, it will be necessary to parameterize our models with data.

Parameterizing and Analyzing Evolutionarily Explicit IPMs
A large literature exists on how to statistically parameterize IPMs (Easterling et al. 2000;Merow et al. 2014;Rees et al. 2014). The vast majority of IPMs have been constructed phenomenologically, using statistical descriptions of observational data. Several authors have shown how fixed and random effects incorporated into these statistical functions can be formulated within IPMs (Childs et al. 2003;Rees and Ellner 2009;Coulson 2012), but additional statistical estimation is required to parameterize the evolutionarily explicit IPMs we have developed. Fitness functions in evolutionarily explicit IPMs can be parameterized using standard general, generalized, and additive regression methods that are routinely used to parameterize phenomenological IPMs (Rees and Ellner 2009). If relatedness information is available and the infinitesimal model is assumed, genetic and phenotypic variances and covariances can be estimated using the animal model (Lynch and Walsh 1998). These quantities can be used to construct the initial distributions of the genetic and environmental components of the phenotype. Parameter estimates of ecological drivers fitted as fixed or random effects in the animal model can be used to parameterize inheritance and development functions for the environmental component of the phenotype. It is consequently possible to parameterize models using our approach with existing methods.
There is also a large literature on how to analyze IPMs (Ellner and Rees 2006;Steiner et al. 2012Steiner et al. , 2014. The majority of these tools-including sensitivity and elasticity analysis of model predictions to transition rates and function parameters (Ellner and Rees 2006;Coulson et al. 2010Coulson et al. , 2011Steiner et al. 2012Steiner et al. , 2014-are likely sufficiently general to be applicable to evolutionarily explicit IPMs. In future work we plan to parameterize models for bird, mammal, and fish species with overlapping generations and to analyze them with existing methods. Once evolutionarily explicit IPMs have been parameterized and analyzed, we will be able to explore how populations, phenotypic characters, and life histories are predicted to respond to a range of environmental changes via plasticity and adaptation.

When Should Evolutionarily Explicit IPMs Be Used to Predict Population Responses to Environmental Change?
Chevin (2015) and Janeiro et al. (2017) speculated that published IPMs that did not include explicit evolutionary processes could provide spurious insight. Three strands of evidence suggest that this speculation may often be unwarranted.
First, the signature of evolutionary change in model predictions is a function of the heritability of the trait: when the phenotypic variance is dominated by the environmental component of the phenotype, then the dynamics of that component will dominate model predictions. Most IPMs to date have been constructed for body weight (Merow et al. 2014), a trait that often has a heritability of !0.2 in vertebrates (e.g., blue tits; Garnett 1981) and often around 0.1 (e.g., bighorn sheep; Wilson et al. 2005). This means that model predictions will be dominated by the dynamics of the environmental component of the phenotype and that a phenomenological statistical approach to parameterizing these models has the potential to capture observed dynamics well.
Second, even when phenotypic traits are heritable, they rarely evolve in the wild as predicted: evolutionary stasis of heritable phenotypic traits in the presence of directional selection is frequently observed in nature (Merilä et al. 2001). When fitness functions are monotonic in the phenotypic value and selection is directional (which is typical for body size [Kingsolver et al. 2001]), then in order to maintain an equilibrium trait distribution the inheritance function must reverse the phenotypic changes caused by selection. Coulson and Tuljapurkar (2008) showed this for the mean phenotypic trait. However, when the genotypephenotype map is additive and there is additive genetic variance for the trait, directional selection is expected to result in evolutionary change, and the inheritance function for the genetic component of the phenotype cannot reverse genetic changes attributable to selection. Unmeasured genetically correlated characters can prevent evolutionary change in these circumstances, although the cases when this is likely to prevent evolution are restrictive, and evidence for such characters playing a major role in limiting evolution in the wild is lacking (Agrawal and Stinchcombe 2009). Assuming that selection on the phenotype has been measured appropriately and is directional, this suggests that the assumption of an additive genotype-phenotype map may be violated, and the mean of the parental and offspring breeding value distributions may not be equal. A mechanism such as overdominance can achieve this (Fisher 1930). Our approach allows the effects of relaxing assumptions of quantitative genetics on evolutionary change to be approximated through the use of phenomenological inheritance functions for the genetic component of the phenotype.
Third, because evolutionary change is rarely observed in the wild when it is predicted, observed phenotype change in natural populations is usually attributable to plasticity (e.g., Ozgul et al. 2009Ozgul et al. , 2010. In these cases, standard nonevolutionarily explicit IPMs have accurately captured observed dynamics (Childs et al. 2003;Ozgul et al. 2010;Merow et al. 2014).
These three strands of evidence suggest that evolutionarily explicit IPMs may frequently not be required to gain useful insight into population responses to environmental change. If there is no statistical evidence of temporal trends in inheritance, development, or fitness function parameters once variation in the ecological environment has been corrected for, then the use of evolutionarily explicit IPMs may result in the construction of unnecessarily complex models. There is often a temptation to include ever more complexity into models, but this comes at the cost of analytical tractability: as more mechanisms or processes are incorporated into models, understanding why a model produces the predictions it does becomes increasingly challenging. However, when evolutionary change is convincingly documented (e.g., Reznick et al. 1997) or is proposed to be a possible mechanism generating rapid phenotypic change (Coltman et al. 2003), the construction of evolutionarily explicit IPMs is advised because the models allow separation of the roles of adaptive and plastic responses to environmental change.
We have shown how evolutionarily explicit IPMs can be constructed, invalidating the criticisms of Chevin (2015) and Janeiro et al. (2017) that IPMs have not been developed to incorporate the character-state approach of quantitative genetics. IPMs that are not evolutionarily explicit have been used to address many questions in ecology, and their application has proven insightful (Merow et al. 2014). They are likely to remain widely used, and we expect this use to result in important new insights. However, we have extended their utility to cases where evolutionary processes are known or proposed to be drivers of phenotypic change.

Conclusions
In this article we have developed a theoretical modeling approach that links demography and quantitative genetics to explore how populations will respond to environmental change. The approach is general, providing formal links between ecology and evolution. Our work builds on a growing literature of developing evolutionarily explicit structured population models. This body of literature shows how flexible IPMs are. They provide a powerful tool with the potential to unify ecology and evolution.