Rapid adaptation of a polygenic trait after a sudden environmental shift

# Rapid adaptation of a polygenic trait after a sudden environmental shift

Kavita Jain and Wolfgang Stephan

Theoretical Sciences Unit,
Jawaharlal Nehru Centre for Advanced Scientific Research,
Jakkur P.O., Bangalore 560064, India
Natural History Museum Berlin, Invalidenstr. 43,
10115 Berlin, Germany

Keywords: polygenic selection, unequal effects, rapid adaptation

Correspondence: jain@jncasr.ac.in, stephan@bio.lmu.de

Abstract: Although a number of studies have shown that natural and laboratory populations initially well-adapted to their environment can evolve rapidly when conditions suddenly change, the dynamics of rapid adaptation are not well understood. Here a population genetic model of polygenic selection is analyzed to describe the short-term response of a quantitative trait after a sudden shift of the phenotypic optimum. We provide explicit analytical expressions for the time scales over which the trait mean approaches the new optimum. We find that when the effect sizes are small relative to a scaled mutation rate, the genomic signatures of polygenic selection are small to moderate allele frequency changes that occur in the short-term phase in a synergistic fashion. In contrast, selective sweeps, i.e., dramatic changes in the allele frequency may occur provided the size of the effect is sufficiently large. Applications of our theoretical results to the relationship between QTL and selective sweep mapping and to tests of fast polygenic adaptation are discussed.

In “What Darwin got wrong”, J. B. Losos (2014) persuasively argues that we can observe evolution in action and, in particular, that evolution can be so rapid that evolutionary and ecological time scales are confluent. The examples range from the peppered moth (Cook et al., 2012), insecticide resistance in Drosophila (Ffrench-Constant et al., 2002), color of field mice (Vignieri et al., 2010), beak size in Darwin’s finches (Grant and Grant, 2008), guppies in Trinidad (Reznick, 2009), and Anolis lizards (Losos, 2009) to name but a few. The rapid changes are responses to natural and human-induced shifts in the environment. The genetic architecture underlying these traits ranges from few genes of major effect to highly polygenic systems (van’t Hof et al., 2011; Linnen et al., 2013; Lamichhaney et al., 2012, 2015). In this paper, we study a model that encompasses a wide range of genetic architectures. Our aim is to understand the genomic signatures of positive selection in these systems that occur after environmental changes with an emphasis on rapid adaptation.

There is a growing body of literature on the detection of adaptive signatures in molecular population genetics. Following pioneering work of Maynard Smith and Haigh (1974), the impact of positive selection on neutral DNA variability (Òselective sweepsÓ) has attracted much interest. This theory has been applied to huge datasets that emerge from modern high-throughput sequencing. A large number of statistical tests have been developed to detect sweep signals and estimate the frequency and strength of selection (Kim and Stephan, 2002; Nielsen et al., 2005; Pavlidis et al., 2010). To find sweep signatures in the genome, strong positive selection is required (with , where is the effective population size and the selective advantage of a beneficial allele) (Kaplan et al., 1989; Stephan et al., 1992). Thus, sweeps are characteristic signals of fast adaptation. They are expected to be found at individual genes or at major loci if a trait is controlled by multiple genes.

This raises the question whether fast evolution is also possible in highly polygenic systems and, if so, which genomic signatures arise. A number of genome-wide association studies (GWAS) have shown that quantitative traits are typically highly polygenic (e.g., Turchin et al. (2012)) and there is growing evidence that the molecular scenario of sweeps only covers part of the adaptive process and needs to be revised to include polygenic selection (Pritchard and Di Rienzo, 2010). As GWAS yield information about the distribution of single nucleotide polymorphisms (SNPs) relevant to quantitative traits (Visscher et al., 2012), it is important to understand the models of polygenic selection in terms of the frequency changes of molecular variants, i.e., with reference to population genetics (Bürger, 2000).

Although the equilibrium structure of the allele frequencies and the stationary variance have been the subject of research in a large number of such studies (Bürger, 2000), the dynamics have received relatively little attention. Perhaps the most well studied model in this context is Lande’s model (Lande, 1983) in which the changes in the allele frequency at a major locus in the background of a large number of minor loci are considered (Chevin and Hospital, 2008; Gomulkiewicz et al., 2010). However in Lande’s model, the background is not explicitly modeled and it is assumed that the background variance does not evolve. Pavlidis et al. (2012) have studied a more detailed model but their analysis is restricted to a very small number of loci.

Following our preliminary study (Jain and Stephan, 2015), here we perform a population genetic analysis of the dynamics of a polygenic trait after a sudden environmental shift of the phenotypic optimum. We consider a quantitative trait that is determined additively by a large number of diallelic loci of unequal effects. The population is assumed to be infinitely large and to evolve under stabilizing selection and mutation. This model was first proposed by Wright (1935) and more recently revisited by N. H. Barton (Barton, 1986; de Vladar and Barton, 2014). We carry out a detailed study of this model to understand the dynamics of allele frequencies as well as of the trait mean and variance. In particular, to describe rapid adaptive evolution, we concentrate on the short-term period after the environmental change, which may be defined as the time until the phenotypic mean reaches a value close to the new optimum.

We find that the short-time dynamics are well described by a directional selection model which is analytically tractable. Working within the framework of this simpler model, we reproduce some results obtained using different models or assumptions (Lande, 1983; Chevin and Hospital, 2008; Jain and Stephan, 2015) when most effects are small relative to a scaled mutation rate (de Vladar and Barton, 2014). In addition, we obtain several new results, in particular when most effects are large.

MODEL

We consider an infinitely large population of diploids in which a trait is determined by diallelic loci. At the th locus, let the allele have an effect and frequency and, correspondingly, the allele have an effect and frequency (de Vladar and Barton, 2014; Jain and Stephan, 2015). The effect at each locus is assumed to be exponentially distributed with mean , as is often the case in quantitative genetic studies (Mackay, 2004; Goddard and Hayes, 2009). Then, on neglecting dominance and epistasis, the trait is determined additively and its mean can be written as

 c1=ℓ∑i=1(γipi−γiqi)=ℓ∑i=1γi(2pi−1) . (1)

For a phenotypic trait under stabilizing selection, assuming that the deviation of the fitness from its optimum is quadratic (Bürger, 2000), we can write the phenotypic fitness as , where measures the strength of stabilizing selection on the trait. Averaging over the phenotypic trait distribution (which is not necessarily Gaussian), one finds the average fitness to be

 ¯¯¯¯w = 1−(s/2)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(z2−c21)+(c21−2zz′+z′2) , (2) ≈ exp[−(s/2)(c2+(c1−z′)2)] , (3)

where is the variance of the trait and is given by (see Appendix A)

 c2=2ℓ∑i=1γ2ipiqi . (4)

When selection is weak and recombination rate is high (linkage equilibrium; see also Appendix B) and mutations are absent, the allele frequency evolves according to Wright’s equation as (Chapter 6, Bürger (2000))

 ˙pi=pi(1−pi)2¯¯¯¯w∂¯¯¯¯w∂pi , (5)

where dot denotes the derivative with respect to time. For the model described above, we then have (de Vladar and Barton, 2014; Jain and Stephan, 2015)

 ˙pi=−sγipiqi(c1−z′)−sγ2i2piqi(qi−pi)+μ(qi−pi) , i=1,...,ℓ . (6)

In the above equation, the first two terms on the RHS follow from (3) and (5). The first term tends to stabilize the phenotypic mean to the optimum trait value and the second term results in the fixation of one of the alleles thus depleting genetic variation. The last term on the RHS models the mutation process restoring genetic variance (Barton, 1986), and describes changes between the and allele at locus with locus-independent, symmetric rate . In the following, we refer to the model defined by (6) as the full model.

In this article, we are interested in the response of the phenotypic mean and variance when the phenotypic optimum is suddenly shifted to a new value in a population initially equilibrated to . A preliminary investigation of such a situation has been carried out numerically by de Vladar and Barton (2014), and our goal here is to provide an analytical understanding of the dynamics. We will thus focus on the dynamics of the allele frequency that evolves according to (6) when the phenotypic optimum , starting from the stationary state frequency of the population that is equilibrated to a phenotypic optimum at zero.

Unless stated otherwise, in this article, we assume that so that the initial mean deviation from the phenotypic optimum, , is negative. We also restrict the magnitude of the shift such that ; this is because the first term on the RHS of (6) acts to minimize the deviation of the trait mean from its optimum. But, from (1), the magnitude of the mean can not exceed and therefore a shift beyond the total effect of all loci is not within evolutionary reach.

STATIONARY STATE OF THE FULL MODEL

Before proceeding further, we briefly review the known results for the stationary state that are pertinent to the later discussion. In the stationary state where the allele frequencies are independent of time, setting the LHS of (6) to zero and performing some simple algebra, we find that the equilibrium allele frequency is a solution of the following equation,

 (7)

where and is the deviation of the mean from the optimum in equilibrium.

When the stationary mean deviation is zero, equation (7) for the equilibrium allele frequency has three solutions that are given by (de Vladar and Barton, 2014)

 p∗i=⎧⎪⎨⎪⎩12, γi<^γ12±12√1−^γ2γ2i, γi>^γ . (8)

While the solution is stable when the effect is smaller than the threshold effect , the two latter ones hold when the effect is larger than . The threshold between large- and small-effect alleles arises because of mutation. When selection is much weaker than mutation, the equilibrium frequency is one half; i.e., when the stationary trait mean is at the optimum, mutation and selection balance each other at an intermediate allele frequency. The value of one half arises because we assumed that mutation in both directions is symmetric (see above). Relaxing this symmetry assumption would lead to different equilibrium frequencies. However, since we focus in the following on the short-term behavior of the dynamics, detailed mutation models are of relatively little importance compared to selection and will not be considered here. For exponentially distributed effects, as the fraction of loci with large effects is given by , most effects are large when , while most effects are small in the opposite parameter regime.

Furthermore, on separating the contribution from loci with small and large effects and using the stationary state frequency (8), we obtain the stationary state variance to be (de Vladar and Barton, 2014; Jain and Stephan, 2015)

 c∗2=ℓ¯γ2[1−(1+^γ¯γ)e−^γ/¯γ] . (9)

The above result shows that when most effects are small, while in the opposite case, it is well approximated by .

The correction to the equilibrium allele frequency in (8) when the stationary mean deviation is nonzero but small can be found in the Appendix B of de Vladar and Barton (2014). For future reference, we note that these corrections are negligible except when the effect is close to the threshold (see Fig. 2 below).

DYNAMICS WHEN MOST EFFECTS ARE SMALL

We now study the dynamics of the allele frequency when most effects are small (). Since the population is initially equilibrated to the phenotypic optimum at zero, the initial frequency at the th locus is well approximated by (8) and close to one half when most effects are small. Then, as at short times, we can neglect the last two terms on the RHS of (6) compared to the first term to get (Jain and Stephan, 2015)

 ˙pi=Sipiqi , (10)

where and . When the final mean deviation is negligible, the final allele frequency is also close to one half and we expect that (10) can describe the bulk of the allele frequency dynamics as is indeed supported by Figs. 1 and 2.

Directional selection model: In the following, we will refer to the model defined by (10) as the directional selection model. Appendix C details the prescription for calculating relevant quantities analytically in this model. As shown in Appendix C, the trait mean can be written as

 c1(t)=ℓ∑i=1γi−ℓ∑i=12γi1+pi(0)qi(0) eβγi/¯γ , (11)

where the parameter defined in (C.6) is proportional to the logarithm of the ratio of allele frequencies and evolves in time as

 ˙β=−s¯γΔc1. (12)

Thus a closed equation for can be obtained by eliminating the mean using (11) and on plugging the result for back in (11), the mean can be found. Furthermore, it is shown in Appendix C that in the directional selection model, the th cumulant of the effect evolves according to (Bürger, 1991)

 ˙cn=−sΔc1cn+1 , n≥1 , (13)

which shows that the th cumulant can be found once is known. Thus to describe the short-time dynamics, we will focus on the key equations (11) and (12).

We now return to (11) for the mean and analyze it when the number of loci is large. When most effects are small (), as mentioned above, the initial allele frequency is one half. Furthermore, as explained in DISCUSSION, the sum over the contribution from individual loci in (11) amounts to averaging over the distribution of effects when the number of loci is large. For large , on approximating the sum in (11) by an integral, we thus get

 c1(t)=ℓ¯γ−ℓ¯γ∫^γ0dγ2γe−γ/¯γ1+eβγ/¯γ . (14)

Writing in the above integral, we get

 c1(t)=ℓ¯γ(1−2Is(β)) , (15)

where

 Is(β)=∫^γ/¯γ0dx xe−x1+eβx . (16)

Our main task is to find the time dependence of using (12) and (15) which is dealt with in Appendix D.

Qualitative patterns of the cumulant dynamics: Before turning to explicit expressions, we first make some general remarks on the behavior of the cumulants in the directional selection model. Since the variance is always nonnegative, it follows from (13) that the magnitude of mean deviation decreases with time. As we are considering the scenario where , the mean deviation increases monotonically with time towards zero as shown in Fig. 1. Due to (12), this also means that initially increases and then saturates to a constant.

Furthermore, due to (12) and (13), the variance can be written as

 c2=−2ℓ¯γ2I′s(β) , (17)

where prime denotes the derivative with respect to . To arrive at the above expression, we have used that due to (15). Equation (17) gives the rate of change of variance as . It is easy to check that the integrand of is always positive from which it follows that . Therefore the directional selection model predicts that when most effects are small, the variance decreases with time as indeed verified in Fig. 1. From this result and (13), it immediately follows that the skewness defined in (A.6) is always negative.

Quantitative dynamics of the mean and higher cumulants: As shown in Appendix D, the mean deviation vanishes exponentially fast,

 Δc1(t)=Δc1(0)e−sc2(0)t , (18)

where is the initial variance given by (9) when most effects are small. The above result has been previously obtained by assuming that the variance and higher cumulants do not evolve in time (Lande, 1983; Chevin and Hospital, 2008; Jain and Stephan, 2015). The basis of this assumption lies in phenotypic data (Lande, 1976) or is motivated by mathematical tractability (Chevin and Hospital, 2008; Gomulkiewicz et al., 2010; Jain and Stephan, 2015). Here, in contrast, we have obtained (18) without any additional assumptions.

Using the solution (18) in (13) for cumulant dynamics, we find that the variance stays at its initial value and the higher cumulants vanish. The corrections to these behavior at late times are given by (D.13) and (D.14), and we find that the variance remains a constant at short times () but decreases thereafter to a smaller value:

 c2(t)c2(0)≈1−3(zfℓ¯γ)2(1−e−sc2(0)t)2 . (19)

The above expression shows that the drop in the variance is larger for large shifts in the optimum and can be quite significant as shown in Fig. 1 for a moderate shift (relative to the maximum mean). A comparison between the full model, directional selection model and the analytical results is also shown in Fig. 1 for mean and variance. The reason for the difference observed in Fig. 1 between the cumulants obtained using (10) and the large- approximation (15) is explained in DISCUSSION.

The main panel of Fig. 1 shows that the mean is very close to the stationary state at time . However, as shown in the inset, the mean continues to change albeit very slowly until when the true stationary state is reached. A similar pattern is seen for variance (and allele frequency shown below in Fig. 2). These observations suggest that the dynamics of the mean can be divided into a short-term phase in which the mean approaches a value close to the new optimum and a long-term phase in which it reaches the exact stationary state. To explore fast evolution, we will concentrate in the following on the short-term phase.

Dynamics of the allele frequencies: From (10) for allele frequency dynamics, we first note that for , all the allele frequencies increase in the short-term phase (note that so far, we have argued that (10) holds for small-effect loci only but in the next section, we will find this to be true for large-effect loci also). This simple result has potentially useful application as explained in DISCUSSION.

Using (18) in (10), the time dependence of the allele frequency for small deviations in the phenotypic optimum can be obtained analytically (Chevin and Hospital, 2008; Jain and Stephan, 2015):

 pi(t)=pi(0)pi(0)+qi(0)exp(−zfγi(1−e−sℓ¯γ2t)ℓ¯γ2) . (20)

Thus in the directional selection model, the allele frequency reaches the stationary state over a time period . These shifts in allele frequency are small to moderate, as verbally predicted by several authors (e.g. Pritchard and Di Rienzo (2010)). However, this directional behavior may change in the long-term period, as we will discuss next.

The allele frequency dynamics in the full model and directional selection model are compared for some loci in Fig. 2, and we see that while the two match well at short times, in the full model, the allele frequency at large times can vary in a nonmonotonic fashion (see inset of Fig. 2) and approaches a stationary state value given by the solution of (7). The long-time dynamics of the allele frequency can be understood by solving (6) with zero mean deviation and the initial condition given by the stationary state solution of the directional selection model (Jain and Stephan, 2015). A succinct way of expressing such a prescription is to write

 Δc1(t)=Δc∗1+Δc1(0)e−sc2(0)t (21)

in (6). Numerical integration of the resulting equation matches well with the full model as shown in Fig. 2. An analytical understanding of the long-time dynamics can also be obtained as described in Appendix E.

DYNAMICS WHEN MOST EFFECTS ARE LARGE

We now turn to the situation when most effects are large (). Due to (8), the initial allele frequencies are close to either zero or one. When the final mean deviation is close to zero, either the final allele frequency stays close to the initial one or sweeps to fixation as shown in Fig. 3. Thus, unlike in the case where most effects are small, in the full model the factor in the last two terms on the RHS of (6) is not negligible. However, the directional selection model defined by (10) describes the short-time dynamics when most effects are large as explained below.

Since we are dealing with the case where , we can ignore the effect of mutations. Furthermore, as can not exceed one, (6) shows that as long as , the second term on the RHS can be neglected to yield (10) for the allele frequency at the th locus. Therefore we expect the directional selection model to apply at short times. But as time progresses, the mean deviation decreases in magnitude and the above condition for the validity of the directional selection model fails.

Using the initial frequency given by (8) for , equation (11) for the mean yields

 c1(t)=ℓ∑i=1γi−∑i∈P+2γi1+p+i(0)q+i(0)eβγi/¯γ−∑i∈P−2γi1+p−i(0)q−i(0)eβγi/¯γ , (22)

where refers to the set of loci with initial frequencies given in (8). From symmetry considerations, we expect the number of loci with initial frequency above and below one half to be equal. Furthermore, for , the initial allele frequency given by (8) can be approximated as . We thus arrive at

 c1(t)≈ℓ∑i=1γi−∑i∈P+2γi1+(2γi^γ)2eβγi/¯γ−∑i∈P−2γi1+(^γ2γi)2eβγi/¯γ . (23)

Then, as in the last section, approximating the sums in (23) by integrals for large , we get

 c1(t)=ℓ¯γ(1−I+l−I−l) , (24)

where

 I+l = 12¯γ∫∞^γdγ2γe−γ/¯γ1+(2γ^γ)2eβγ/¯γ (25) = ∫∞2/αdx xe−x1+(αx)2eβx , (26)

and

 I−l = 12¯γ∫∞^γdγ2γe−γ/¯γ1+(^γ2γ)2eβγ/¯γ (27) = ∫∞2/αdx xe−x1+(αx)−2eβx . (28)

In the above equations, the parameter .

Qualitative patterns: In the directional selection model, as (12) and (13) hold irrespective of whether most effects are small or large, the mean deviation and the allele-frequency-dependent variable increase with time as discussed in the last section, and the variance is given by . However, unlike when most effects are small, here the rate of change of variance, is not a monotonic function of time; this is because the sum is negative when is small and positive for larger . Thus the directional selection model predicts that initially the variance increases with time and then decreases towards the stationary state at the new phenotypic optimum.

Quantitative dynamics: Equation (24) along with (12) forms a closed system of equations and allows us to find the mean. The integrals are estimated in Appendix F; we find that the integral does not contribute substantially to the mean since it includes the loci with initial frequency close to one and we are dealing with the case when the initial mean deviation is negative (see also (32) below), and is well approximated by given by (F.8). When the initial deviation from the phenotypic optimum is small, using (F.14) and (F.21), we find that at large times, the mean deviation decreases exponentially fast,

 Δc1(t) ∼ e−sℓ¯γ2(lnα)−1ρ(lnρ)2t , (29) = e−szf¯γ(lnα)−1(lnρ)2t , (30)

where is the phenotypic shift relative to the maximum mean magnitude. The above result shows that time scale over which the mean approaches the stationary state depends weakly on the number of loci unlike the case when most effects are small (see (18)). However as the mean of the distribution is large here, the relaxation in the directional selection model occurs faster than in the small-effects case.

When most effects are large, the variance changes considerably even when the relative phenotypic shift is quite small as Fig. 4 shows. The peak variance can be estimated by the stationary state variance in the directional selection model. At large times, using (F.16) and (F.21), we find that

 c(max)2=ℓ¯γ2(lnα)−1ρ(lnρ)2 (31)

and therefore . This result can also be seen by noting that the allele frequencies at loci with initial frequency close to zero shift to intermediate values at short times. Then from (4), the maximum variance is obtained when the allele frequency is close to one half leading to .

The dynamics of the allele frequency at long times are discussed in Appendix E and at short times in the next section.

WHEN DO SELECTIVE SWEEPS OCCUR?

Below we obtain a criterion on the size of effects for which the allele frequency at a major locus can sweep. When , at short enough times, keeping only the lowest order terms in on the RHS of the full model (6) and neglecting the effect of mutations, we get

 ˙pi≈−sγipi2(2Δc1(0)+γi) . (32)

When the initial mean deviation , the above equation shows that only the loci with effect can potentially sweep since their frequency increases with time. On repeating the above analysis for loci with initial frequency close to one, we find that the frequency at such loci does not sweep to zero. Thus, when , only the loci with small initial frequency and effect are likely to undergo large changes in the frequency. These criteria are necessary but as Fig. 3 shows, they are not sufficient for a selective sweep to occur.

For the major loci that meet the above necessary conditions for selective sweeps, the second and the third term on the RHS of the full model (6) can be neglected in comparison to the first term thereby reducing (6) to the directional selection model (10). For the rest of the (major) loci, the allele frequency changes are not appreciable and the frequencies remain close to zero or one, a solution satisfied by (10). Let denote the time when the allele frequencies evolving according to the directional selection model equilibrate. Figure 3 shows that the mean deviation obtained using the full model is close to zero at time . Then for subsequent times , we can ignore mutations and set in the full model (6) to obtain (Jain and Stephan, 2015)

 ˙pi=sγ2i2piqi(2pi−1) , t>~t . (33)

The above differential equation is subject to the initial condition which is the allele frequency in the steady state of the directional selection model.

From the definition (C.6) of , we have

 pi(~t)1−pi(~t)=pi(0)qi(0)exp(γi~β¯γ) (34)

which simplifies to give the frequency as

 pi(~t)=[1+qi(0)pi(0)exp(−γi~β¯γ)]−1 , (35)

where . It is clear from (33) that if the frequency is greater than one half, the allele frequency at the th locus will increase monotonically towards one. Thus we predict that the allele frequency at a locus would sweep when . Combining this criterion with the necessary condition mentioned above, we find that a sweep can occur when the effects lie in the following range:

 ¯γ~βln(q−i(0)p−i(0))<γi<2zf. (36)

This expectation is indeed seen to be consistent with the numerical solution of the full model as shown in Fig. 3 except for one case.

When most effects are small: While considering the dynamics of mean and variance when most effects are small in an earlier section, we ignored the contribution of (few) loci that have large effect. However the results obtained earlier can be utilized to understand the dynamics of the allele frequency at major loci also. As (D.11) shows, the auxiliary parameter in the stationary state of the directional selection model is . Using this in (36), we find that a selective sweep can occur at a major locus if its effect

 2zf>γi > ℓ¯γ2zfln(q−i(0)p−i(0)) , (37) = c2(0)|Δc1(0)|ln(q−i(0)p−i(0)) , (38)

which matches the criterion (26b) of Chevin and Hospital (2008). Furthermore, using (8), we find that

 2c2(0)|Δc1(0)|ln(2γi^γ)<γi<2zf . (39)

For the parameters in Fig. 1 where the average effect size and the threshold effect , the above criterion requires that the effect size at a major locus be . For exponentially distributed effects, as is assumed in this article, the probability for the effect size to exceed this criterion is extremely low ( for the parameters in Fig. 1) and therefore selective sweeps at major loci are prevented when most effects are small (Chevin and Hospital, 2008; Pavlidis et al., 2012).

When most effects are large: The criterion for short-time sweeps is more involved in this case; using (8) in (36), we obtain

 2zf>γi > 2¯γ~βln(2γi^γ) , (40) ≈ 2¯γln(2¯γ^γ)ln(ℓ¯γzf)ln(2γi^γ) , (41)

where we have used (F.20) for in the last expression. For the parameters in Fig. 4, using in (40), we find that for a selective sweep to occur, is required whose probability is . In other words, we expect as many as sweeps to occur for these parameter values over the time scale the mean deviation approaches zero.

DISCUSSION

Dynamical regimes in the full model: Figure 3 shows an example of the dynamics of the mean and allele frequencies when an infinitely large population in linkage equilibrium and equilibrated to a phenotypic optimum is suddenly shifted to a new optimum. We observe a short-term phase where the mean deviation increases quickly to a value close to zero followed by a long-term phase where the mean deviation undergoes minor changes at a slow pace. The former phase where the adaptation process is rapid is the focus of this article since: (i) much of the action happens in this early phase (for example, in Fig. 3, the allele frequencies at out of loci were close to the stationary state at the end of the short-term phase), (ii) these time scales are experimentally observable, and very importantly, (iii) all alleles increase in frequency in a synergistic manner.

Directional selection model: Although the full model defined by (6) can be investigated numerically (de Vladar and Barton, 2014), as the allele frequency dynamics at a single locus are related to all the other loci via the phenotypic mean (first term on the RHS of (6)), it is difficult to obtain analytical results. Moreover, the full model results in equations for the cumulants that do not close, and therefore one resorts to an often unjustified and uncontrolled approximation of setting all the cumulants above a certain number to zero (see Barton and de Vladar (2009) for a discussion of this approach). This approximation, however, works quite well when most effects are small (Jain and Stephan, 2015) but fails completely when most effects are large (see supplemental file S1 for a detailed explanation).

Here we have devised an approximate model (directional selection model) defined by (10) in which the allele frequencies are coupled but nonetheless, one can make analytical progress without truncating the cumulants arbitrarily. The central equations in the directional selection model are (11) and (12) that relate the mean deviation and an auxiliary parameter which is a function of a single allele frequency. As explained in Appendices D and F, it is possible to obtain explicit expressions for the time dependence of and thereby that of the mean. Higher cumulants such as the variance can be then found using (13).

Distribution of effects: All the numerical examples shown in the figures use a single realization of the effect at a locus, i.e., one set of effects are chosen from an exponential distribution. For additive quantities such as mean and variance that are obtained by summing the contribution over a number of loci, we expect that a single realization of effects is a good representative when the number of loci are sufficiently large (in physics literature, this concept is known as self-averaging (Castellani and Cavagna, 2005)). A similar idea has been used recently in a model of stabilizing selection with the same effects at all loci in which the phenotypic mean for one set of initial allele frequencies is approximated by the corresponding average (Charlesworth, 2013). Figure S1 shows that the quantitative agreement between the results for the mean deviation obtained for a single realization of phenotypic effects for and the analytical result for large gets better as the number of loci increase.

Main results: The key conclusions of this article are discussed below and summarised in Table 1.

Response of the mean to a sudden shift of the phenotypic optimum: How does the time scale over which the mean reaches the stationary state depend on various factors such as the number of loci, initial phenotypic mean deviation and size of the effects? We find that the mean deviation becomes negligible on a time scale that is inversely proportional to the number of loci () when most effects are small but depends weakly on when most effects are large (cf. (18) and (30)). This difference in the behavior may be understood as follows: when most effects are small, there is sufficient genetic variation initially for selection to act on because a large number of loci are initially polymorphic. Therefore, with increasing , the equilibration time decreases. In contrast, when most effects are large, the initial genetic variation is negligible and the adaptation dynamics are determined by the size of the effects. As the availability of large-effect loci depends on the mean of the effect distribution, it follows that the larger the mean, the faster the adaptation dynamics.

Figure 5 shows that the relationship between the equilibration time and the initial mean deviation depends on the size of effects (see also Gomulkiewicz et al. (2010)). When most effects are large, selective sweeps occur at some loci that result in drastic changes in the allele frequencies thus accelerating the approach to the phenotypic optimum. But when most effects are small, the equilibration time remains roughly constant as selective sweeps do not occur but small changes in the allele frequencies at a large number of loci drive the adaptation process. Figure 5 also shows that when the number of loci contributing to a trait are the same, adaptation is faster when the effect size is larger, which makes intuitive sense.

However one may also compare situations in which a large number of small effects control a phenotypic trait with the one in which few loci with large effects contribute to a trait. An example shown in Fig. 6 suggests that the adaptation at short times is faster when effects are small but at longer times, larger effects lead to rapid adaptation. This conclusion is borne out by our analysis as well: on comparing (18) and (30), we find that for identical and , the time scales are inversely proportional to and therefore it takes shorter time to equilibrate in the large-effect case.

An important omission in the above discussion is random genetic drift as we have considered infinitely large populations. Some progress in this regard has been made recently (Matuszewski et al., 2015; Bod’ová et al., 2016); however, a detailed analysis is currently not available and would be interesting to consider in the future.

Dynamics of the variance and allele frequencies: Besides the dynamics of the mean, we also studied the time dependence of the variance. We find that at short times, the variance does not vary much and decreases monotonically with time when most effects are small. This behavior can be explained by the fact that the allele frequencies show small to moderate changes. In the opposite case, the variance changes considerably and its variation is nonmonotonic in time - it increases at short times and then decreases to the stationary state value. The selective sweeps at major loci are responsible for this effect. The rate of change of average fitness given by (C.9) receives contributions from the square of the mean deviation and the rate of change of variance. Thus when most effects are large, at short times, the average fitness increases more slowly in the short-term phase than at longer times (Gomulkiewicz et al., 2010).

The directional selection model describes the bulk of the mean dynamics in the short-term phase. Therefore at late times where the mean deviation is close to zero, the allele frequency at each locus evolves independently (the first term on the RHS of (6) is zero) and the dynamics are much slower as described in Appendix E. At short times, the allele frequency dynamics involve the same time scales as for the mean and variance.

Applications: The analysis of our model relates to various important questions that are currently discussed in evolutionary genetics.

Rapid adaptation: The approximations we have presented here describe the short-term evolution of a phenotypic trait after a sudden environmental change very well. We have shown that the mean of a phenotypic trait may respond very fast after an environmental shift. In the case that most effects are small, this is possible because the time to the new optimum is inversely proportional to the number of loci controlling the trait. In the opposite case of mostly large effects, rapid fixations (leading to selective sweeps) may produce fast phenotypic responses. In the examples of fast adaptation mentioned in INTRODUCTION both of these extremes and combinations thereof have occurred.

QTL and selective sweep mapping: Selective sweep mapping has been used to dissect QTL (Svetec et al., 2011; Rubin et al., 2012; Axelsson et al., 2013; Wilches et al., 2014). However, the success of this approach was varied. Nonetheless, there seems to be a tendency that it worked better in domesticated than in natural populations (Rubin et al., 2012; Axelsson et al., 2013), probably due to the action of artificial selection during domestication. In artificial selection, the shift to the new optimum may be larger than under natural conditions; our criterion (36) thus indicates an enhancement of sweeps in domesticated populations.

More generally, our analysis provides also some new insights into the question whether selective fixations (and thus sweeps) occur at QTL. While Chevin and Hospital (2008) predicted that the probability of selective sweeps is extremely low at QTL (based on a model with one major locus and infinitely many minor loci), others have found sweeps at appreciable frequencies using simulations of various multi-locus models (Pavlidis et al., 2012; Wollstein and Stephan, 2014). The prediction of Chevin and Hospital is consistent with our study for mostly small-effect loci and small shifts in the phenotypic optimum.

Testing for fast polygenic adaptation: Using standardized frequencies, one can construct tests to detect SNPs that deviate strongly from neutral population structure (e.g. Coop et al. (2010)). However, this approach is only working if there exist relatively large extended gradients of ecological variables, which may not be the case in rapid adaptation. Shortly after a population occupies a new habitat, we expect that the allele frequency shifts between the parental and derived populations are relatively small. This also means that available software, such as BayeScan-like methods (Foll and Gaggiotti, 2008; Riebler et al., 2008), is not able to detect significant frequency shifts for individual SNPs between populations. Therefore, for detecting small allele frequency shifts after environmental changes in fast adapting populations, it may be best to consider the frequency shifts of alleles simultaneously at all loci involved (instead of individual SNPs). Our results suggest that this may be a promising approach, as all + alleles shift their frequencies in the same direction in the short-term phase, which should increase the power of the test. In human population genetics, this approach has been used to analyze recent polygenic adaptation for instance of height (Turchin et al., 2012).

Acknowledgements: We thank the ICTS, Bangalore, for their hospitality during the Second Bangalore School on Population Genetics and Evolution (ICTS/Prog-popgen/2016/01) that facilitated our discussions. Furthermore, we are grateful to Graham Coop and two anonymous reviewers who made valuable suggestions to improve the paper. The research of WS was supported by the German Research Foundation DFG (grant Ste 325/17-1 from the Priority Program 1819).

APPENDIX

## Appendix A Cumulants of the effects and their dynamics

The effect at the locus is distributed according to a Bernoulli distribution

 p(ei)=piδei,γi2+qiδei,−γi2 . (A.1)

The generating function of the distribution of the effects at the th locus is given by

 Fi(ξ)=∑eip(ei)e−eiξ=pie−γi2ξ+qieγi2ξ . (A.2)

The logarithm of the generating function defines the th cumulant at locus as (Sornette, 2000)

 lnFi(ξ)=∞∑n=1c(i)n(−ξ)nn!=ξγi2+ln(qi+pie−γiξ) . (A.3)

The th cumulant is then obtained as

 c(i)n=(−1)ndnlnFi(ξ)dξn∣∣ξ=0 . (A.4)

In linkage equilibrium, since the cumulants for the -locus problem are additive, we get

 cn=2ℓ∑i=1c(i)n , (A.5)

where the factor of on the RHS takes care of diploidy. Using the last three equations above, we find that the mean and variance are given by (1) and (4), respectively. The skewness which is a measure of asymmetry of a distribution can also be found and given by

 c3=2ℓ∑i=1γ3i(qi−pi)piqi . (A.6)

## Appendix B Quasi-linkage equilibrium approximation

For the additive genotype-phenotype map, the trait where if the allele is present at the locus . Following Barton and Turelli (1991), we first find the association coefficients by writing where and expanding the relative fitness to lowest orders in ,

 w({X,X∗})¯¯¯¯w = 1−(s/2)(z−z′)2¯¯¯¯w ≈ 1+ℓ∑i=1ai,ϕ(δi−⟨δi⟩)+ℓ∑i=1aii,ϕ(δ2i−⟨δ2i⟩) + ℓ∑i=1a∗i,ϕ(δ∗i−⟨δ∗i⟩)+ℓ∑i=1a∗ii,ϕ(δ∗2i−⟨δ∗2i⟩) + ∑i

In the above equation, the average is taken with respect to the joint distribution of the frequencies at all the loci and the coefficients

 ai,ϕ=−s(c1−z′)γi, aii,ϕ=−(s/2)γ2i, aij,ϕ=ai,j=−sγiγj . (B.3)

The average allele frequency change after selection (and recombination) is given by (see (10a), Barton and Turelli (1991))

 ˙pi=ai,ϕ⟨δ2i⟩+aii,ϕ⟨δ3i⟩+∑j≠iaj,ϕ⟨δiδj⟩+∑j≠iajj,ϕ⟨δiδ2j⟩+∑j≠kakj,ϕ⟨δiδjδk⟩ . (B.4)

All but the first two terms on the RHS vanish in linkage equilibrium (LE); here we focus on the lowest order corrections to LE which are contained in the third and fourth terms. Using and , it is easy to show that . Inserting these expressions in (B.4), we verify that the first two terms on the RHS match the corresponding ones in (6). The two-locus linkage disequilibrium evolves according to (see (12), Barton and Turelli (1991))

 ˙⟨δiδj⟩=−rij⟨δiδj⟩+2aij,ϕ⟨δ2i⟩⟨δ2j⟩+higher order terms ,i≠j , (B.5)

where denotes the recombination rate between the loci and and we have dropped terms that are nonlinear in the parameters and on the RHS of the above equation. Then at quasi-linkage equilibrium (QLE) where the LHS is of or higher, we have

 ⟨δiδj⟩QLE ≈ 2aij,ϕrij⟨δ2i⟩% LE⟨δ2j⟩LE (B.6) = −2sγiγjrij(piqipjqj)|LE , (B.7)

where the allele frequencies in the last equation are in LE.

Using the above result in (B.4) and assuming that all the recombination rates are equal to , we find that within QLE, the allele frequency obeys the following equation:

 ˙pi≈−sγipiqi(c1−z′)(1−src2,LE)−sγ2i2piqi(qi−pi)(1−src3,LEγi) (B.8)

where and are the variance and skewness, respectively, in linkage equilibrium. Thus, as expected, the corrections to the model defined by (6) are of order which can be neglected when selection and linkage are weak.

## Appendix C Directional selection model

The dynamics of the allele frequency are described by (10) where the effective selection coefficient depends on all the allele frequencies through the mean . This property makes it difficult to solve for the allele frequencies; however, we note that for two arbitrary loci and , (10) gives

 ˙pi˙pm=dpidpm=SipiqiSmpmqm , (C.1)

which, on using the identity , simplifies to

 dlnϱidlnϱm=γiγm , i≠m , (C.2)

where . As the above equation does not involve the phenotypic mean, it can be easily solved to give

 (C.3)

We remark that the simple idea used in (C.1) to eliminate the global variable is quite similar to that employed in single-locus directional selection models to get rid of the average fitness (see, for example, p. 73, Charlesworth and Charlesworth (2010)).

Equation (C.3) is useful since it allows us to express the mean in terms of a single allele frequency, say, at locus . Using (C.3) in (1), we can write the mean as

 c1(t) = ℓ∑i=1γi(ϱi−1ϱi+1) (C.4) = ℓ∑i=1γi−ℓ∑i=12γi1+ϱi(0)(ϱm(t)ϱm(0))γiγm , (C.5)

which simplifies to give (11). In the above equation,

 β(t)=¯γγmln(ϱm(t)ϱm(0))=¯γγmln(pm(t)pm(0)qm(0)qm(t)) . (C.6)

Note that in the expression (C.5), although the dynamics of the mean deviation are determined by that of the allele frequency , the rest of the frequencies are not redundant as the mean depends on their initial state also.

Furthermore, from the evolution equation (10) for the allele frequency, we also have , where the mean is given by (C.5); we have thus obtained a closed equation for . Using the result for in (C.5), the dynamics of the phenotypic mean can be calculated. This knowledge then allows us to find higher cumulants using (13) which is derived below.

To obtain the evolution equation (13) for the th cumulant, we take the derivative with respect to time in (A.3) and use (10) for the allele frequency dynamics to obtain (Bürger, 1991)

 ∞∑n=1˙c(i)n(−ξ)nn!=sγipiqiΔc1(1−e−γiξ)qi+pie−γiξ . (C.7)

But, on taking a derivative of (A.3) with respect to , after some simple algebra, we also have

 ∞∑n=1c(i)n+1(−ξ)nn!=−γipiqi(1−e−γiξ)qi+pie−γiξ . (C.8)

On comparing the last two equations and using (A.5), we arrive at the promised result (13). The rate of change of the logarithm of the average fitness can be found using the above results and given by

 dln¯wdt=˙¯w¯w=−s2˙c2+c2(sΔc1)2 . (C.9)

## Appendix D Cumulant dynamics when most effects are small

Here we study the differential equation (12) obeyed by ,

 ˙β=−sℓ¯γ2[1−ρ−2Is(β)] , (D.1)

where and the integral

 Is = ∫^γ/¯γ0dx xe−x1+eβx (D.2) = ∫∞0dx xe−x1+eβx−∫∞^γ/¯γdx xe−x1+eβx . (D.3)

The first integral on the RHS is exactly solvable in terms of special functions and we have

 I(1)s = ∫∞0dxxe−x1+eβx , (D.4) = Ω2[ψ(1)(Ω+1/2)−ψ(1)(Ω+1)]