# Rapid Adaptation of a Polygenic Trait After a Sudden Environmental Shift

^{*}Theoretical Sciences Unit, Jawaharlal Nehru Centre for Advanced Scientific Research, Bangalore 560064, India^{†}Natural History Museum Berlin, 10115, Germany^{‡}Biocenter, University of Munich, 82152 Planegg, Germany

- 1Corresponding author Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore, 560064, India. E-mail: jain{at}jncasr.ac.in

## 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 timescales over which the trait mean approaches the new optimum. We find that when the effect sizes are small relative to a scaled mutation rate, small to moderate allele frequency changes 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.

- polygenic selection
- unequal effects
- rapid adaptation

IN “What Darwin got wrong,” Losos (2014) persuasively argues that we can observe evolution in action and, in particular, that evolution can be so rapid that evolutionary and ecological timescales 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; Lamichhaney *et al.* 2012, 2015; Linnen *et al.* 2013). In this article, 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 the 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 data sets 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 *s* is 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 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 *z* is determined by diallelic loci. At the *i*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 *z* is determined additively and its mean can be written as(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 *s* 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(2)(3)where is the variance of the trait and is given by (see Appendix A)(4)When selection is weak and recombination rate is high [linkage equilibrium (LE); see also Appendix B], and mutations are absent, the allele frequency evolves according to Wright’s equation as (Bürger 2000, Chap. 6)(5)where the 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)(6)In the above equation, the first two terms on the right-hand side (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 *i* 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 cannot exceed and therefore a shift beyond the total effect of all loci is not within evolutionary reach.

## Results

### 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 left-hand side (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):(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)(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 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 Figure 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 *i*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)(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 Figure 1 and Figure 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(11)where the parameter *β* defined in (C6) is proportional to the logarithm of the ratio of allele frequencies and evolves in time as(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 *n*th cumulant of the effect evolves according to (Bürger 1991)(13)which shows that the *n*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 the *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(14)Writing in the above integral, we get(15)where(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 toward zero, as shown in Figure 1A. 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(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 Figure 1B. From this result and (13), it immediately follows that the skewness defined in (A6) is always negative.

#### Quantitative dynamics of the mean and higher cumulants:

As shown in Appendix D, the mean deviation vanishes exponentially fast,(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 behaviors at late times are given by (D13) and (D14), and we find that the variance remains a constant at short times but decreases thereafter to a smaller value:(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 Figure 1B, 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 Figure 1 for mean and variance. The reason for the difference observed in Figure 1 between the cumulants obtained using (10) and the large approximation (15) is explained in the *Discussion*.

The main panel of Figure 1A 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 Figure 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 also find this to be true for large-effect loci]. This simple result has potentially useful application as explained in the *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):(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 Figure 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 Figure 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(21)in (6). Numerical integration of the resulting equation matches well with the full model as shown in Figure 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 Figure 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 cannot 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 *i*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(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(23)Then, as in the last section, approximating the sums in (23) by integrals for large we get(24)where(25)(26)and(27)(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 toward 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 (F8). When the initial deviation from the phenotypic optimum is small, using (F14) and (F21), we find that at large times the mean deviation decreases exponentially fast,(29)(30)where is the phenotypic shift relative to the maximum mean magnitude. The above result shows that the timescale 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 Figure 4B shows. The peak variance can be estimated by the stationary state variance in the directional selection model. At large times, using (F16) and (F21), we find that(31)and therefore This result can also be seen by noting that the allele frequencies at loci with an 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(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 frequency. These criteria are necessary, but as Figure 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)(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 (C6) of *β*, we have(34)which simplifies to give the frequency as(35)where It is clear from (33) that if the frequency is greater than one half, the allele frequency at the *i*th locus will increase monotonically toward 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:(36)This expectation is indeed seen to be consistent with the numerical solution of the full model as shown in Figure 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 used to understand the dynamics of the allele frequency at major loci too. As (D11) 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(37)(38)which matches the criterion (26b) of Chevin and Hospital (2008). Furthermore, using (8), we find that(39)For the parameters in Figure 1 where the average effect size and the threshold effect the above criterion requires that the effect size at a major locus be >0.81. 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 Figure 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(40)(41)where we have used (F20) for in the last expression. For the parameters in Figure 4, using in (40), we find that for a selective sweep to occur, is required whose probability is In other words, for these parameter values, we expect as many as sweeps to occur over the timescale the mean deviation approaches zero.

### Data availability

The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article.

## Discussion

### Dynamic 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 LE 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 Figure 3, the allele frequencies at 13 out of 20 loci were close to the stationary state at the end of the short-term phase), (ii) these timescales 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 Material, 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 Appendix D and Appendix 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 in File 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 increases.

### Main results

The key conclusions of this article are discussed below and summarized in Table 1.

#### Response of the mean to a sudden shift of the phenotypic optimum:

How does the timescale 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 timescale 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 Figure 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 timescales are inversely proportional to and therefore it takes a 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 (C9) 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 timescales 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 quickly 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 the 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 also provides 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 a neutral population structure (*e.g.*, Coop *et al.* 2010). However, this approach only works if there are 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 of height, for instance (Turchin *et al.* 2012).

## Acknowledgments

We thank the International Centre for Theoretical Sciences, 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 article. The research of W.S. was supported by the Deutsche Forschungsgemeinschaft (grant STE 325/17-1 from the Priority Program 1819).

## Appendix A: Cumulants of the Effects and Their Dynamics

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

The generating function of the distribution of the effects at the *i*th locus is given by(A2)The logarithm of the generating function defines the *n*th cumulant at locus *i* as (Sornette 2000)(A3)The *n*th cumulant is then obtained as(A4)In LE, since the cumulants for the *n*-locus problem are additive, we get(A5)where the factor of two 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

## Appendix B: Quasi-LE Approximation

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

In the above equation, the average is taken with respect to the joint distribution of the frequencies at all the loci and the coefficients(B3)The average allele frequency change after selection (and recombination) is given by [see (10a) in Barton and Turelli 1991](B4)All but the first two terms on the RHS vanish in 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 (B4), 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) in Barton and Turelli 1991](B5)where denotes the recombination rate between the loci *i* and *j* and we have dropped terms that are nonlinear in the parameters *s* and *r* on the RHS of the above equation. Then at quasi-LE (QLE) where the LHS is of or higher, we have(B6)(B7)where the allele frequencies in the last equation are in LE.

Using the above result in (B4) and assuming that all the recombination rates are equal to *r*, we find that within QLE, the allele frequency obeys the following equation:(B8)where and are the variance and skewness, respectively, in LE. 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 *i* and *m*, (10) gives

which, on using the identity simplifies to(C2)where As the above equation does not involve the phenotypic mean, it can be easily solved to give(C3)We note that the simple idea used in (C1) 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, *e.g.*, Charlesworth and Charlesworth 2010, p. 73).

Equation C3 is useful since it allows us to express the mean in terms of a single allele frequency, say, at locus *m*. Using (C3) in (1), we can write the mean as(C4)(C5)which simplifies to give (11). In the above equation,(C6)Note that in the expression (C5), 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 also depends on their initial state.

Furthermore, from the evolution Equation 10 for the allele frequency, we also have where the mean is given by (C5); we have thus obtained a *closed* equation for Using the result for in (C5), 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 *n*th cumulant, we take the derivative with respect to time in (A3) and use (10) for the allele frequency dynamics to obtain (Bürger 1991)(C7)But, on taking a derivative of (A3) with respect to *ξ*, after some simple algebra, we also have(C8)On comparing the last two equations and using (A5), 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

## Appendix D: Cumulant Dynamics When Most Effects Are Small

Here we study the differential Equation 12 obeyed by *β*,

where and the integral(D2)(D3)The first integral on the RHS is exactly solvable in terms of special functions and we have(D4)(D5)(D6)where and are the second derivatives of the logarithm of the gamma function [(6.4.1) in Abramowitz and Stegun 1964]. Since when most effects are small, the second integral on the RHS of (D3) can be estimated by carrying out an integration by parts and we have(D7)where The above integral is thus exponentially small in and can be neglected in comparison to Using the series expansion for the polygamma function for small arguments [(6.4.10) and (6.4.11) in Abramowitz and Stegun 1964], we obtainAs both the initial and stationary state frequencies are close to one half, the function given by (C6) stays close to zero. Using the result (D8) in the Equation D1 for we get(D10)where As ignoring the term on the RHS of the above equation and integrating, we obtain the zeroth order solution given by(D11)where we have used that An approximate solution of (D10) can be found iteratively by writing in (D10) and retaining the leading order term in to obtain(D12)which yields Using this result in (D8), we get Equations 15 and 13 then yield the first two cumulants to be(D13)and

(D14)## Appendix E: Long-Time Dynamics of the Allele Frequency

When most effects are small, as shown in Figure 1A, the mean deviation is close to zero when But the exact stationary state for the mean and the allele frequency (shown in Figure 2) is reached much later at A similar qualitative behavior of the allele frequencies when most effects are large is seen in Figure 3. To understand the slow dynamics at long times when the effect size can be small or large, here we study (6) for the allele frequency dynamics in the full model approximating the mean deviation by its stationary state value We thus have

An exact solution of the above equation requires solving a cubic equation which is possible to do but not particularly enlightening.

To find the behavior of the allele frequency close to the stationary state, we write in the above equation and obtain(E2)where is the exact solution of (7) with nonzero and the constants and are functions of As we are interested in large times where the allele frequency is close to the stationary state, the first term on the RHS of the above equation can be neglected and we arrive at a first-order differential equation for which can be easily solved. For large times, we then obtain where(E3)Thus after the directional selection phase, the allele frequency at the locus with effect reaches a stationary state over a timescale For small mean deviation from the optimum, which results in negligible we can use (8) to obtainThe above equation shows that when most effects are small, the dynamics at long times are driven by mutations (recall that while in the opposite case, the effect size continues to play a role even at long times.

## Appendix F: Cumulant Dynamics When Most Effects Are Large

Consider first the integral defined by (26) that includes the contribution from loci with initial frequency close to one and given by

where is large. The last integral can be written as an exponential integral and we have(F3)(F4)on using that the argument of is small (Abramowitz and Stegun 1964). We next consider the integral defined by (28) that contains the contribution from loci with small initial frequencies,(F5)(F6)(F7)where we have used that to obtain the last expression and defined(F8)As shown in Table S1 in File S1, the integral decreases as or slower. Thus, for large *α*, Equation (24) for the mean gives(F9)In the above equation, the integral does not contribute to the mean since, as discussed after (32), the allele frequency at loci with high initial frequencies does not change significantly while the mean is evolving.

The time dependence of *β* is obtained from Equation 12,(F10)The function *β* increases with time and saturates to its steady state value where Correspondingly, as Table S1 in File S1 indicates, the integral is close to zero when and approaches *ρ* as Thus at short times, we can neglect the term on the RHS of (F10), and at large times, we can expand about *ρ* to write (F10) as Denoting the time below and above which these two solutions are valid by and setting we obtainwhere and Using this result in (F9), we find the mean deviation to beUsing the above expressions in (13), we find that the variance is given byTo express the mean and the variance in terms of the model parameters, we first need to determine which is a solution of the equation The integral is not exactly solvable but we can estimate it by a saddle-point method [see (F17) below] (Arfken 1985). Denoting the integrand of in (F8) by and the maximum of by we obtain(F17)where is a solution of the equation (F18)As our numerics indicate that is large and constant for a given *α*, using the above equation, we obtain for large *α*. Using this solution in (F17), after some simplifications, we find that(F19)Setting in the above expression yields(F20)Taking the derivative with respect to *β* in (F19), we arrive at(F21)Our numerics are consistent with the above functional form and a numerical fit shows

## Footnotes

Supplemental material is available online at www.genetics.org/lookup/suppl/doi:10.1534/genetics.116.196972/-/DC1.

*Communicating editor: G. Coop*

- Received October 18, 2016.
- Accepted March 7, 2017.

- Copyright © 2017 by the Genetics Society of America