The ﬁ tness costs of adaptation via phenotypic plasticity and maternal effects

Summary 1. Phenotypes are often environmentally dependent, which requires organisms to track environmental change. The challenge for organisms is to construct phenotypes using the most accurate environmental cue. 2. Here, we use a quantitative genetic model of adaptation by additive genetic variance, within- and transgenerational plasticity via linear reaction norms and indirect genetic effects respectively. 3. We show how the relative influence on the eventual phenotype of these components depends on the predictability of environmental change (fast or slow, sinusoidal or stochastic) and the developmental lag s between when the environment is perceived and when selection acts. 4. We then decompose expected mean fitness into three components (variance load, adaptation and fluctuation load) to study the fitness costs of within- and transgenerational plasticity. A strongly negative maternal effect coefficient m minimizes the variance load, but a strongly positive m minimises the fluctuation load. The adaptation term is maximized closer to zero, with positive or negative m preferred under different environmental scenarios. 5. Phenotypic plasticity is higher when s is shorter and when the environment changes frequently between seasonal extremes. Expected mean population fitness is highest away from highest observed levels of phenotypic plasticity. 6. Within- and transgenerational plasticity act in concert to deliver well-adapted phenotypes, which emphasizes the need to study both simultaneously when investigating phenotypic evolution.


Introduction
Phenotypes are complex, built from genetic and non-genetic components which are modified by physiological, developmental and environmental cues (West-Eberhard 2003). This modification is often provoked by external stimuli, whose effects are filtered by alleles, genotypes and phenotypes to affect fitness (Lewontin 1974;Ricklefs & Wikelski 2002;Coulson et al. 2006;Kokko & Lopez-Sepulcre 2007). Each path has dynamical consequences for the rate of phenotypic evolution (Day & Bonduriansky 2011) because of the lags between an environmental stimulus and the evolutionary response it provokes (Ricklefs & Wikelski 2002). Here, we focus on the consequences of genetic evolution and two forms of phenotypic flexibility: transgenerational parental effects (Falconer 1965;Mousseau & Fox 1998;R€ as€ anen & Kruuk 2007) and within-generation phenotypic plasticity (Scheiner 1993;Lande 2009).
Phenotypic plasticity is the ability of a genotype to alter its resultant phenotype in response to environmental change (Scheiner 1993;Lande 2009). This ability can be heritable (Scheiner 2002) and enables, for example, individual great tit Parus major mothers to match the timing of breeding with maximal food abundance (Charmantier et al. 2008), orange-tip butterflies Anthocharis cardamines to adjust to large-scale temperature gradients (Phillimore et al. 2012) and bird populations to considerably reduce their extinction risk (Vedder, Bouwhuis & Sheldon 2013). In theoretical work, phenotypic plasticity has been predicted to play a major role in accelerating adaptation to novel (Lande 2009) or sinusoidal environments (Tufto 2000) as well as to marginal habitats (Chevin & Lande 2011).
Maternal effects (Falconer 1965;Mousseau & Fox 1998) are the influence of the mother's genotype or phenotype on her offspring (Wolf & Wade 2009). Maternal effects are obviously transgenerational, using some aspect of maternal condition in the parental generation to maximize fitness in the current one (Jablonka & Lamb 2005) and are likely to evolve if parents can process environmental signals more accurately than the current generation (Uller 2008). One example of a predictable environment is winter moths Operophtera brumata living on oak Quercus robur. Oak trees have consistent phenology from one year to the next (Crawley & Akhteruzzaman 1988) and so the moth can use her parental experience to help ensure her eggs hatch when food (young oak leaves) is most abundant (van Asch, Julkunen-Tiito & Visser 2010). Such effects confer fitness benefits: seeds of the monocarpic herb Campanulastrum americanum grown in the same light environment as experienced by their mother had higher survival probability than those grown in the opposite environment (Galloway & Etterson 2007). Similarly, relative fitness of the bryozoan Bugula neritina was increased if warmer temperatures were experienced by both mother and offspring (Burgess & Marshall 2011).
The interplay between within-and transgenerational plasticity matters because the stage when an individual processes an environmental cue affects its biological response (Ricklefs & Wikelski 2002). Environmental change propagates through phenotypes to affect mean population fitness in theoretical (Greenman & Benton 2005), experimental (Petchey 2000) and empirical (Garcia-Carrera & Reuman 2011) settings, as well as across generations (Galloway & Etterson 2007;Burgess & Marshall 2011). In ecological scenarios, environmental change is often positively autocorrelated (Halley 1996;Vasseur & Yodzis 2004), meaning environments in successive time intervals are more similar than would be expected by chance alone. In this scenario, a clear adaptive benefit of transgenerational plasticity is expected (Uller 2008), yet the maternal phenotype will likely have been in part determined by a phenotypically plastic response within the preceding generation. The two types of plasticity are interdependent.
Relatively little is known about how additive genetic variation, within-generation phenotypic plasticity and transgenerational maternal effects interact during adaptation to a changing environment (Chevin, Collins & Lefevre 2013). Wolf & Brodie (1998) showed how stabilizing selection on a maternally influenced trait favours a genetic correlation among direct genetic and indirect maternal effects that is opposite in sign to the maternal effect coefficient. Lande (2009) showed a two-phase adaptation, first by phenotypic plasticity and subsequently by genetic assimila-tion. Vedder, Bouwhuis & Sheldon (2013) parameterized Chevin, Lande & Mace's (2010) mechanistic model of adaptation via genetic change and phenotypic plasticity for the Wytham woods great tit population. Here, we use a quantitative genetic model (Hoyle & Ezard 2012) to identify the fitness implications of phenotypes constructed from different combinations of additive genetic variance, within-and transgenerational plasticity.

Materials and methods
We use a quantitative genetic model of adaptation via an additive genetic component, within-generation phenotypic plasticity and transgenerational maternal inheritance via a constant maternal effect coefficient m, which was derived by Hoyle & Ezard (2012) and merges ideas from Kirkpatrick & Lande (1989), Lande & Kirkpatrick (1990) and Lande (2009). m (Kirkpatrick & Lande 1989) represents the path from maternal phenotype after selection in the previous generation z Ã tÀ1 to offspring phenotype without direct genetic transmission. This is an indirect genetic effect (Moore, Brodie & Wolf 1997;McGlothlin & Brodie 2009;Hadfield 2012) because the maternal effect is determined by the mother's phenotype, which has a heritable genetic component. We assume therefore that the phenotypes in previous generations contribute to the current phenotype under selection and so the fitness is assigned to the current generation (Wolf & Wade 2001).
Within Wolf & Wade's (2009) definition of a 'true' maternal effect as one with a causal link between maternal genotype or phenotype and the offspring phenotype, various subcategories have emerged: transgenerational effects that enhance fitness have been termed 'adaptive maternal effects' (Marshall & Uller 2007) while R€ as€ anen & Kruuk (2007) define positive maternal effects as ones that (could) accelerate microevolution because a positive m means that, on average and all other inheritance forms held equal, larger-than-average females produce larger-than-average offspring (Kirkpatrick & Lande 1989). Under the same conditions, a negative m means larger-than-average females produce smaller-than-average offspring. Using the same model as here, Hoyle & Ezard (2012) showed how a negative maternal effect is 'adaptive' in a relatively stable environment, but a positive maternal effect is 'adaptive' during adaptation following a rapid shift in environment.
Our aim here is to quantify adaptive maternal effects (sensu Marshall & Uller 2007) under simple models of environmental change and then investigate how additive genetic variance, phenotypic plasticity and the maternal effects combine to underpin phenotypic evolution. We focus on a single, normally distributed phenotypic trait z, such as body size, asking how it influences itself in future generations (e.g. Falconer 1965). Generations are discrete and non-overlapping. The phenotype of an individual at time (≡generation) t is: where a t is the additive genetic component (elevation, breeding value) of the phenotype and b t the slope of the linear reaction norm of the plastic phenotypic response to the environment e t (Lande 2009). z Ã tÀ1 is the phenotype after selection of a selected parent contributing to the next generation. Note that we do not consider fertility selection here (see Hadfield (2012) for discussion). There is a lag, measured in fractions of a generation s, between development and selection. If s is small, then this lag is short and selection and development are closer together. Large s potentially increases the within-generation mismatch between expressed and target phenotype if the environment is not constant. Finally, e t is the residual component of phenotypic variation, which we assume to have mean zero without loss of generality.
The phenotypic variance, r 2 zt , of z t is Following Lande (2009), we assume a constant, equilibrium reference environment e = 0 where variance in the phenotypic plasticity reaction norms is minimized such that the covariance between a t and b t is necessarily 0. We also assume constant additive genetic variances of a t and b t (G aa and G bb , respectively). If z t is normally distributed with variance r 2 z , then the mean fitness (with respect to phenotype distribution) in the offspring generation (Lande 1979) is: where c ¼ 1=ðx 2 þ r 2 z Þ and the optimum phenotype, h t = A+Be t , is assumed to be a linear function of the environment at time t. x is the 'width' of the fitness function so that c represents the strength of stabilizing selection. A, B, W max and x are constants.
Assuming that a t and b t are bivariate normally distributed, the per generation change in their population means, a t and b t , is (Lande 1979 where the average phenotype after selection in the previous generation is where we have assumed that there is no fertility selection, in order to set In the next subsection, we will see how these components evolve in response to a changing environment.

E N V I R O N M E N T A L C H A N G E A N D E X P E C T E D D Y N A M I C S
To study the interplay of within-and transgenerational plasticity under micro-and macroenvironmental change, we allow the environment to follow a noisy long-term variation e t = U t +ξ t , where U t is the long-term variation and ξ t a Gaussian stationary autocorrelated random process with mean zero, variance r 2 n and autocorrelation q s over the interval s. We assume q s is negligible across generations.
Due to the good predictive ability of the expectation of the stochastic dynamics (Hoyle & Ezard 2012), we take the stochastic expectation of D a, D b and z Ã t averaged over the distribution of ξ (treating a t , b t and z Ã tÀ2 as fixed). We therefore find: Under the same conditions, the expected phenotypic variance at time tEðr 2 zt Þ satisfies: Assuming small noise relative to the width of the fitness function (r 2 n ( x 2 ), we can derive an expression for expected mean fitness based on eqn (3) by averaging over the phenotypic distribution in the expected environment: approximates c (Hoyle & Ezard 2012). In deriving eqn (5), we have again treated a t , b t and z Ã tÀ2 as fixed when averaging over the distribution of environments. Setting W max = 1, eqn (5) can be written as the product of three terms: We consider these three components as the variance load (eqn 6, Lande & Shannon 1996), adaptation (eqn 7) and a fluctuation load (eqn 8). Note that the fluctuation load is caused purely by microenvironmental fluctuations (cf. Lynch & Lande 1993) and it vanishes when r 2 n ¼ 0. The biological relevance of accurate adaptation is clear. Studying the variance load is important because its consequences for fitness depend on the type of environmental change: it increases extinction risk and reduces fitness in constant or unpredictable environments, but has the opposite effect in highly variable but predictable ones (Lande & Shannon 1996). There is no guarantee that an expressed, plastic phenotype will match its target ('perfect plasticity') because of the lag (here, s) between the processing time (environment of development) and the point when selection acts (Gavrilets & Scheiner 1993;Tufto 2000;Lande 2009). If the correlation between the environments of development and selection is weak, perhaps to due longer s, inaccurate processing of the cue and/or greater environmental variation, the fitness costs of mismatched plasticity are larger (Reed et al. 2010). Note that the fluctuation load requires microenvironmental variation, (r 2 n [ 0) and includes both within-and transgenerational components (eqn 8), emphasizing the benefits of studying both simultaneously.
We generated sequences of environmental change around the reference environment e = 0. We considered slow-and fast-changing environments, which either cycle deterministically as a sine wave or flip stochastically as a Poisson process. We superimposed microenvironmental fluctuations with r 2 n ¼ 2 on top of longerperiod shifts between environmental extremes, d. We selected d = 10 as the distance between environmental extremes to represent a distinct macroenvironmental signal beyond the normal, background range experienced by different generations through the microenvironmental noise from r 2 n and to match the assumption of Lande (2009) that G bb d 2 /(G aa +G bb d 2 ) is near 1. d>10 would push the optimal m towards lower values and vice versa. This distinction lets us consider the role of long-scale, seasonal changes on top of microenvironmental fluctuations and lets us investigate fitness costs of mismatched plasticity (eqn 8).
We generated deterministic, cyclic environments using sine waves with peak-to-peak amplitude d, that is, sinð2pmtÞ d 2 for the environment of selection and sin 2pmðt À sÞ ð Þ d 2 for the environment of development. We used m = 0Á1 and m = 0Á01, that is, 10 and 100 generations per complete sinusoidal cycle for fast-and slow-flipping environments, respectively. To obtain waiting times between successive flips from e = d/2 to e = Àd/2 in stochastic environments, we generated random numbers from the exponential distribution with mean l set at l = 50 and l = 5 to represent slow-and fast-flipping environments, respectively. The environment can therefore change before development and selection, between development and selection, or after development and selection. l represents the average time that the environment remains in a given state so that the environment is, on average, back in its original state after 2l generations. In these stochastic environments, the same environmental sequence was saved, and all combinations of other parameters under investigation were subjected to it. Results were time-averaged over 100 000 generations (after a burn-in of 50 000 generations) to find the Eð WÞ delivered by different phenotype constructions. We update the expected phenotypic variance (eqn 4) in each generation and calculate it using the environments experienced in both the present (offspring) and previous (maternal) generation. Since r 2 z ! AE1 as m?AE1 (r 2 z is undefined at m = AE1), we restrict ourselves here to À0Á7 <m <0Á7 and consider increments of 0Á05 across this range. Naturally, individuals do not knowingly act to optimize expected population mean fitness; our goal here is to examine the fitness costs that arise from phenotypes constructed using different combinations of additive genetic variance, phenotypic plasticity and maternal effects. The MATLAB scripts we used are available at http:// dx.doi.org/10.6084/m9.figshare.894438.

Results
Representative dynamics of phenotypic evolution are given in Fig. 1. The phenotypic response to the changing environment is delivered by both the additive genetic component (Eð a t Þ, Fig. 1c,d) and within-generation phenotypic plasticity (Eð b t Þ, Fig. 1e,f). In sinusoidal environments, the frequency of the evolution of Eð b t Þ is twice that of Eð a t Þ and also of the expected phenotype, Eð z t Þ. Eð a t Þ can either be in phase (Fig. 2a) or antiphase (Fig. 2b), which is determined by different combinations of the maternal effect coefficient m and the lag between the environments of development and selection s.
As expected, expected mean population fitness Eð WÞ is higher for more strongly positive m in slowly changing, more predictable environments (Fig. 3a). In slowly changing environments, s makes no real impact on Eð WÞ unlike in fast-changing ones (compare left with middle and right columns in Fig. 3). Longer s in fast-changing sinusoidal environments means highest Eð WÞ when m is nonnegative, but closer to zero (Fig. 3e). In the fast-changing stochastic environment, the same m maximizes Eð WÞ for all s considered, but larger s incurs greater fitness costs (Fig. 3i).
Eð WÞ consists of three components: variance load (eqn 6), adaptation (eqn 7) and fluctuation load (eqn 8). In all cases, the adaptation term (matching of the phenotype with its target) is the major influence on Eð WÞ (Fig. 3c,g,  k). The impact of s is through adaptation (Fig. 3g,k) and the fluctuation load (Fig. 3h,l), not the variance load (Fig.  3f,j). In fast-changing sinusoidal environments, the adaptation term is maximized by either positive or negative m if s is short or long, respectively (Fig. 3g). A negative m minimizes the variance load (Fig. 3b,f,j), but a positive m reduces the plasticity penalty because this term is then closer to 1 (Fig. 3d,h,l). Neither positive nor negative m is consistently optimal for each of the three components of Eð WÞ.
In sinusoidal environments, the observed phenotype is delivered via contributions from m and the time-averaged expected mean plasticity (Eð bÞ, Fig. 4a,c). Although the time average of the expected additive genetic component Eð aÞ is zero (Fig. 4), it changes across generations (Figs 1  and 2). In the stochastic environment, the observed phenotype is delivered via contributions from Eð aÞ, Eð bÞ and m (Fig. 4b,d), with Eð aÞ particularly important when m<0. The dependence of Eð aÞ on m in stochastic environments arises because the evolutionary dynamics are not ergodic under the stochastic forcing (Fig. 5). The optimal phenotype is not at highest observed levels of within-generational plasticity (Fig. 4d). In the fast stochastic environment, Eð bÞ is highest when s is shortest (Fig. 4). If the environment changes rapidly, then lower m and higher Eð bÞ compared to those for the slowly varying environments deliver higher Eð WÞ (Fig. 4).

Discussion
There can be many paths to the same phenotype (West-Eberhard 2003). Interrogating our quantitative genetic model of adaptation by within-and transgenerational plasticity indicates that neither a positive nor negative maternal effect coefficient is optimal to minimize fitness costs from the variance load, adaptation term and fluctu- ation load (eqn 5, Fig. 3), and the optimal phenotype is constructed away from the highest observed levels of within-generational phenotypic plasticity, particularly in slowly changing environments (Fig. 4). This is the first attempt, to the authors' knowledge, to simultaneously dissect fitness costs during adaptation to predictably changing environments via additive genetic variance, withingeneration phenotypic plasticity and transgenerational maternal effects.
The optimal strength of transgenerational plasticity depends on the within-generation processes of phenotypic plasticity b t and the lag from juvenile development to adult selection s. The interplay among the maternal effect coefficient m, expected mean plasticity Eð bÞ and s is particularly clear when comparing fast-changing stochastic and fastchanging sinusoidal environments (Fig. 3), where Eð bÞ is higher than in slowly varying environments (Fig. 4). The slope of the linear reaction norm Eð b t Þ evolves at twice the frequency of the phenotype and additive genetic  Fig. 1. Expected mean fitness is expressed as a proportion of W max with variance load, adaptation and fluctuation load as proportions whose product is expected mean fitness. component (Fig. 2) because there are two turning points in each sinusoidal environmental cycle. Next to the extreme environments, selection changes direction, so negative (or at least reduced) plasticity is favoured because the delay in selection response due to s means a change in the direction of selection between development and selection. van Dooren (2001) discussed this predictability near peaks and troughs in sinusoidal environments, albeit in a different genetic model of reaction norm evolution. In stochastic environments, Eð WÞ and s (Fig. 3i) are negatively correlated: longer s means more chance that the environment that triggers the plastic response is different from that when selection acts (Chevin, Collins & Lefevre 2013).
In predictable environments, the parental condition (here, the phenotype in the parental generation) is a reliable indicator of a desirable phenotype in the current generation, and positive maternal effects are expected to evolve (Uller 2008). Our model is consistent with this expectation: Eð bÞ is likely to be higher when s is shorter (Fig. 4d) and, although we assume m is fixed, transgenerational plasticity will likely be greater for species living in slow-rather than fast-changing environments (Figs 3 and  4). The latter case is the situation when light environments of an understorey herb (Galloway 2005;Galloway & Etterson 2007) or the timing of food availability for moths (van Asch et al. 2010) are predictable across generations. Were either environment to become less predictable from one generation to the next, the model suggests that m would be reduced, or even negative.
The time average of the additive genetic component Eð aÞ is not zero in stochastic environments (Fig. 4b,d) because the irregularity of the environmental sequence is amplified by 'evolutionary momentum' from the transgenerational time lags in the response to selection (Kirkpatrick & Lande 1989). The amount of delay, or inertia, depends on selection pressure, particularly when the trait is genetically correlated with others in the offspring (Kirkpatrick & Lande 1989;Lande & Kirkpatrick 1990;Townley & Ezard 2013). Kirkpatrick & Lande (1989) argued that the effect of maternal traits on phenotypic evolution decreases as a geometric series in m t when selection stops. Hoyle & Ezard (2012) suggested that the time-scale for adaptation to a new equilibrium is of the order of tens of generations, being slower for more negative m. There are two key differences between the assumptions of Kirkpatrick & Lande (1989) and our model here. First, selection does not stop in our simulations and acts at each time step of each trajectory. The chain of transgenerational inheritance runs all the way back to the first generation, since Eð a t Þ, Eð b t Þ and Eð z Ã tÀ1 Þ depend in part on the environment in the previous generation and therefore, in turn, all previous environmental states. Secondly, the environment is not constant, it switches indefinitely between two extremes with an average time between flips of either 5 (fast) or 50 (slow) generations. Given the time-scales for adaptation to a new equilibrium under this model (Hoyle & Ezard 2012), the biological system is in a permanent condition of adjustment to a new equilibrium. This is how the non-ergodicity is able to take effect. Further evidence in support of this argument is that Eð aÞ is very close to zero for m>0, when the maternal effects accelerate adaptation towards the shifting optimum (R€ as€ anen & Kruuk Maternal effect coefficient Fig. 4. The time-averaged contributions of the expected additive genetic component (Eð aÞ, black) and expected phenotypic plasticity [Eð bÞ grey] to the phenotype after selection differ substantially among slow (top), fast (bottom), sinusoidal (left) and stochastic (right) environments. Short, intermediate and long lags between juvenile development and adult selection (s = 0Á05,0Á20,0Á35) are denoted in dotted, dashed and long-dashed lines, respectively. The points denote the combination of components that maximizes expected mean fitness with circles, diamonds and crosses for s = 0Á05,0Á20,0Á35, respectively. The units of Eð aÞ and Eð bÞ are absolute numbers and should be considered relative to the target A = 0 and B = 2. Eð aÞ should be considered in absolute terms; it is negative because the random initial environmental state was Àd 2 . Parameters as Fig. 1. 2007) and the system is close to ergodic (Figs 1 and 4). Our model also includes within-generation phenotypic plasticity assuming linear reaction norms (Lande 2009), but supports the general conclusion of (Kirkpatrick & Lande 1989): maternal inheritance fundamentally alters phenotypic evolution. Our model does not incorporate any costs of phenotypic plasticity (van Tienderen 1991;. Taking this step would lower expected mean plasticity Eð bÞ and could slow adaptation, unless m increased to compensate. We chose not to include costs of phenotypic plasticity to facilitate comparison of within-and transgenerational plasticity (Uller 2008), the latter of which we also assumed cost-free (Kirkpatrick & Lande 1989). Selection is also assumed to be density-independent, and we ignore demographic processes, even though density-dependent compensation can buffer populations against any phenological mismatch (Reed et al. 2013). We also make the strong assumption of a constant m. Empirical and laboratory systems have documented substantial context dependence in the strength of m, mediated through, for example, elevated physiological levels (Sheriff, Krebs & Boonstra 2010), altered phenology (Galloway, Etterson & MnGlothlin 2009) or changing strengths of delayed density dependence (Inchausti & Ginzburg 1998) via interactions with the population's age structure (Plaistow & Benton 2009). Evolvable maternal effects would be far more biologically realistic than constant m, but is non-trivial in this framework because of, for example, the need to account for the covariance between a t and the product of the maternal coefficient and offspring phenotype, which in turn complicates the selection gradient calculations and necessitates additional approximations beyond those used here. One of those assumptions concerns the inclusion of background environmental noise, r 2 n , and not doing this would prevent the decomposition of Eð WÞ (eqns 6-8). Using the expected values rather than the actual changes allows us to look at the decomposition of fitness into its three components and so investigate the interplay between within-generation and transgenerational plasticity in more detail. A fixed m also makes comparison with statistical approaches easier (McGlothlin & Brodie 2009). Parameterizing empirical models of the genotypephenotype map can, in principle, be achieved with sensitivity analysis (Coulson et al. 2006).
Phenotypes are complex structures, built from various genetic and non-genetic components. Here, we studied how within-and transgenerational plasticity facilitate adaptation to a changing environment. While within-and transgenerational plasticity are flexible ways of delivering dynamic phenotypes, the fitness costs or benefits of each pivots on the type of environmental change experienced by individuals in the population. Our results emphasize the flexibility of these underlying phenotypic components, and so help explain the wide range of maternal effect coefficients reported empirically (R€ as€ anen & Kruuk 2007). The non-ergodicity of the expected additive genetic component Eð aÞ with the stochastic forcing used is clear since the cumulative mean Eð aÞ over 40 000 generations is not converging on zero (dotted), unlike the cumulative mean Eð aÞ at generation 50 001 from 40 000 randomly initiated sequences (at either e = d/2 or e = Àd/2, dashed lines). Note that this the cumulative mean of the randomly initiated sequences contains no information on the evolutionary process and is included as a visual aid.