Polygenic adaptation to an environmental shift: temporal dynamics of variation under Gaussian stabilizing selection and additive effects on a single trait

Predictions about the effect of natural selection on patterns of linked neutral variation are largely based on models involving the rapid fixation of unconditionally-beneficial mutations. However, polygenic traits are expected to evolve by small allele frequency changes at many loci. Here, I use explicit forward simulations of a single trait with additive-effect mutations adapting to an optimum shift. In this model, sweeps of large-effect mutations occur from a nearly equal number of new mutations and standing variants, except when the trait is highly polygenic and fixations from standing variation dominate. Large-effect mutations contribute to adaptation even when the trait is not mutation-limited. However, such mutations do not sweep to fixation on the time scale of adaptation for polygenic traits. The signature of hitch-hiking, measured using standard summaries of linked variation, are approximately identical for the two classes of sweeps, which is due to the equilibrium frequencies of standing variants at mutation-selection equilibrium at the time of the environmental shift. I also show that the mean times to adaptation and patterns of linked selection are robust to certain details of the distribution of effect sizes, and that the proportion of new mutations of large effect is a critical parameter affecting patterns of linked selection.

The effect of natural selection on linked neutral variation has been extensively studied for the case of directional selection on mutations with direct effects on fitness (e.g., Kaplan et al. 1989;Stephan et al. 1992;Wiehe and Stephan 1993). This framework leads to a natural simulation scheme using the structured coalescent (Kaplan et al. 1988), which has been widely used to study the power of various approaches to detect recent sweeps from new mutations (Fay and Wu 2000;Kim and Nielsen 2004), from standing variation (Innan and Kim 2004;Hermisson and Pennings 2005;Przeworski et al. 2005), from new mutations occurring at a fixed rate in the genome (Braverman et al. 1995;Przeworski 2002), or to test methods to distinguish between various models of adaptation (Garud et al. 2015;Schrider and Kern 2016).
The model of Gaussian stabilizing selection around an optimal trait value differs from the standard model in that mutations affect fitness indirectly via their effects on trait values. For the additive model of gene action considered here, and considering a single segregating mutation affecting the trait, the mode of selection is under-or overdominant in a frequency-dependent manner (Robertson 1956;Kimura 1981). This model has been extended to multiple mutations in linkage equilibrium by several authors (Barton 1986;de Vladar and Barton 2014;Stephan 2015, 2017b).
The equilibrium conditions of models of Gaussian stabilizing selection on traits have been studied extensively (Bürger 2000, chapters 4 and 5). In general, the dynamics are quite complicated, with many possible equilibria existing for the case of many biallelic loci with equal effect sizes and no linkage disequilibrium (Barton 1986). It is also common to assume that the forward and backward mutation rates per locus are equal (Barton 1986;de Vladar and Barton 2014;Stephan 2015, 2017b). Under these assumptions, and assuming distributions of mutational effects symmetric 0, large-effect variants (e.g., those with effect sizes . 2 ffiffiffi 2 p ffiffiffiffiffi V S p , where V S is the variance of the Gaussian fitness function) will be near the boundaries while small-effect variants will be at frequencies near one-half (de Vladar and Barton 2014;Jain and Stephan 2017b). Here, large and small effect is with respect to the effect of a variant on the genetic load of a population (de Vladar and Barton 2014).
While the fitness effects of individual mutations on trait values affect their fixation probabilities, change in the mean phenotype of a population depends on the additive genetic variance (Robertson 1960). When most mutational effects are small and additive, fixations require on the order of the population size in generations because phenotypic change proceeds via the fixation of small-effect mutations, primarily by genetic drift (Robertson 1960). Recent theoretical work has attempted to clarify when sweeps should happen and when adaptation should proceed primarily via subtle allele frequency shifts. Chevin and Hospital (2008) considered the case of a single mutation with a large effect on fitness in a highly polygenic background evolving according to an infinitesimal model. The authors found that sweeps stall at intermediate frequencies because frequency shifts in the polygenic background contribute to adaptation. Under models of linkage equilibrium, additive mutational effects, and equal rates of forward and back mutation at a biallelic locus (Barton, 1986;de Vladar and Barton 2014), polygenic traits adapt quickly to a sudden shift in the optimum via directional selection (Jain and Stephan 2017b). In an infinitely large population, mutations that are rare at the time of the optimum shift may fix if their effect sizes are not overly large relative to the magnitude of the shift. The number of large-effect sweeps during adaptation depends on the magnitude of the shift and the average effect size of segregating variants (Jain and Stephan 2017b). After the directional phase, selection becomes disruptive, and mutations affecting fitness are fixed or lost to reduce the genetic load of the population.
Under a model of a trait with a small number of phenotypic classes, Höllinger et al. (2019) describe the dynamics of mutations following an optimum shift for traits with low mutation rates and for highly polygenic traits. The key parameter in their model is Q ¼ 4Nm, where m is the mutation rate relevant to the trait. When Q ≲ 1, adaptation primarily occurs via complete sweeps. At intermediate values ðQ 10Þ, partial and complete sweeps occur by the time the population has adapted. When Q 100, adaptation (defined as when mean fitness has recovered following the optimum shift) proceeds via frequency shifts at many loci.
While the work described above identifies the conditions where sweeps are expected, we do not have a picture of the dynamics of linked selection during adaptation to an optimum shift. In large part, the difficulty of analyzing models of continuous phenotypes with partial linkage among sites has been an impediment to a theoretical description of the process. In general, the standard model of a single trait with additiveeffect mutations and Gaussian stabilizing selection assumes linkage equilibrium (or quasi-linkage equilibrium) (Turelli 1984;Barton 1986;de Vladar and Barton 2014;Stephan 2015, 2017b). Höllinger et al. (2019) were able to accommodate partial linkage by simplifying how mutations affect phenotype and focusing on the dynamics up until a particular mean trait value was first reached. In their simplest model, an individual is either mutant or nonmutant, and therefore there are only two phenotypes possible.
Here, I use explicit forward-time simulations to describe the average dynamics of linked selection during the adaptation of a single trait under "real" stabilizing selection (Johnson and Barton 2005) as it adapts to a single, sudden shift in the optimum trait value. These simulations accommodate genetic drift and partial linkage, and are also able to track the dynamics of neutral variants over time. By restricting mutations affecting the trait to specific "loci" (within which linkage is still relatively loose) and allowing neutral mutations to occur over much larger genomic intervals containing the loci, I describe the physical distances over which hitchhiking during polygenic adaptation leaves detectable signatures. The simulations conducted here are therefore analogous to those used to study the spatial dynamics of linked selection via the structured coalescent (Kaplan et al. 1988;Braverman et al. 1995;Kim and Stephan 2002;Przeworski 2002). The key conceptual difference is that the model of adaptation is changed from constant directional selection to the sudden optimum shift models involving a continuous trait considered in de Vladar and Barton (2014) and Stephan (2015, 2017b). I also investigate the effect of the recombination rate on the time to adaptation and the fixation time of beneficial mutations with respect to the mean time required to adapt to the new optimum.

Materials and Methods
Modeling stabilizing selection I modeled a single trait under real stabilizing selection (Johnson and Barton 2005). Mutations affecting trait values arise at rate m per haploid genome per generation according to an infinitely many sites scheme (Kimura 1969). For the majority of results, the effect sizes of new mutations on trait values, g, are drawn from a Gaussian distribution with mean zero and SD s g . Mutations have additive effects on trait value and therefore an individual's genetic value, z, is the sum of all effect sizes in that individual.
Here, I use the term "locus" to refer to a continuous genomic region within which mutation and recombination events occur uniformly. Within a locus, mutations occur at positions according to a uniform continuous distribution according to an infinitely many sites scheme. Thus, each mutation results in a biallelic variant and, in the case of trait-affecting mutations, the derived allele affects trait values. What I refer to here as mutations are typically referred to as loci in much of the theoretical literature (Robertson 1956(Robertson , 1960Turelli 1984;Barton 1986;de Vladar and Barton 2014;Stephan 2015, 2017b).
Traits are under Gaussian stabilizing selection, such that 2V S , where z o is the optimal trait value and V S is the sum of the variance in fitness plus the environmental variance in phenotype (Bürger 2000, p. 160). Figure 1 shows a schematic of the model. For all simulations performed here, I use V S ¼ 1.
I modeled an environmental challenge as a sudden optimum shift, where the optimum trait value changed from It is important to note that I considered all of the heritable variation for the trait to be modeled in the genomic regions that are explicitly simulated. Thus, the approach is similar in spirit to that of de Vladar and Barton (2014), but with partial linkage. An alternative would be to allow for a genetic background that also evolves, for which we are not tracking mutation fates. Chevin and Hospital (2008) used the latter approach to mathematically model the dynamics of largeeffect mutations in an infinitesimal background and Stetter et al. (2018) used a simple version of this method to simulate the dynamics of quantitative traits evolving under truncation selection.
Forward simulation schemes I ran all simulations using two different Python packages (see Software availability below) based on the C++ library fwdpp (Thornton 2014). For a given diploid population size, N, I simulated for 10N generations with z o ¼ 0, at which point the optimum shifted and evolution continued for another 10N generations.
Simulating large genomic regions with only selected variants: To study the dynamics of mutations affecting trait values over time, I evolved populations of size N ¼ 5; 000 diploids, where mutations affecting trait values occur uniformly (at rate m) in a continuous genomic interval in which recombination breakpoints arise according to a uniform Poisson process with a mean of 0.5 recombination breakpoints per diploid. The mutation rates used were 2:5 3 10 24 , 10 23 , and 5 3 10 23 , which is the total mutation rate per haploid genome. The total mutation rate per diploid, U, was 2m. These mutation rates corresponded to Q ¼ 4Nm values of 5, 20, and 100, respectively, meaning sweeps were expected to be high frequency, mixes of partial and complete sweeps, and adaptation primarily by allele frequency changes, respectively, as the population approached the new optimum (Höllinger et al. 2019). The three postshift optima used were z o ¼ 0:1, 0.5, and 1. For all combinations of m and z o , V S ¼ 1 and s g ¼ 0:25. At mutation-selection equilibrium, these parameters result in an equilibrium genetic variance given by the "House of Cards" approximation, which is 4m for the definition of mutation rate and the V S used here, and ignoring the contribution of genetic drift (Turelli 1984). With drift, the expected V G differs from the deterministic approximation by a factor of 1=½1 þ V S =ðNs 2 g Þ (Bürger (2000), p. 270, Equation 2.8), which is 1 for the parameters used here. For the low m and low V S used here, the expected genetic variance is therefore small and new mutations are more likely to have large effects relative to standing variation.
For the mutation rates and s g defined above, the mutational variances of the trait are 2ms 2 g , or 3:25 3 10 25 , 1:25 3 10 24 , or 6:25 3 10 24 , respectively, for each mutation rate. In practice, mutational variances are often estimated with respect to the environmental variants, which poses a small issue in relating the parameters to available estimates. Here, I simulated all traits with V S ¼ 1 and did not explicitly model random effects on trait values. If we were to simulate a trait with environmental variance equal to the expected genetic variance and hold V S ¼ 1 instead, the heritability of the trait would be one-half and the evolutionary dynamics would be unaffected because the contribution of the environmental variance to V S would be small (because the genetic variances simulated here are small with respect to the total V S ). Assuming a hypothetical simulation of a trait with heritability equal to one-half, these parameters result in a ration of the mutational variance to the environmental variance of Oð10 22 Þ, which is the upper limit of the ranges reported based on experimental results [Lynch (1988) and Falconer and Mackay (1996), p. 349]. Below, I describe simulations varying the distributions of effect sizes, thus changing the mutational variance.
For all combinations of m and z o , various summaries of the genetic variation (V G ; z, etc.) in the population were recorded every generation. In total, I ran 1024 replicates of each parameter combination. For the first 256 replicates, the frequency trajectories of all mutations were recorded.
Simulating a 10-locus system with neutral and selected variants: For multilocus simulations, a locus has scaled neutral mutation rate u ¼ 4Nm n ¼ 1000 and scaled recombination rate r ¼ 4Nr ¼ 1000, where m n is the neutral mutation rate per gamete at a locus and r is the mean number of recombination events per diploid at a locus. Mutation and recombination events occur uniformly along a locus, and each locus is separated by 50 cM. For these simulations, I performed 256 simulation experiments per parameter combination. Figure 2 shows how a locus is broken up into windows for analysis. Mutations affecting the trait occurred in the sixth out of 11 equal-sized windows in a locus and I analyzed each window separately. Thus, each window had u ¼ r 90 and mutations affecting trait values were clustered in the middle of each locus (and were intermixed with neutral mutations). In these simulations, the total mutation rate affecting the trait, m, was the sum over loci and the rate per locus was equal ðm=10Þ.
At each locus, mutations affecting the trait occurred only in the middle window ( Figure 2); therefore, the mean number of recombination events per diploid was 0:0045 in the middle window where trait-affecting variants arose. Similarly, the mean number of new mutations per diploid at a given locus affecting the trait was m=5. For the largest mutation rate used here ðm ¼ 0:005Þ, the ratio of recombination events to new mutations affecting the trait in this window was nine to one. The entire genome consisted of 10 such loci, for a total mutation rate of m and a total u ¼ 10 4 .
For a model of a single trait under Gaussian stabilizing selection with a constant optimum, the selection coefficient Simons et al. (2018), see also Kimura and Crow (1978)]. Here, V S ¼ 1, and therefore the relevant scaled strength of selection acting on a segregating variant was Ng 2 . For many of the results presented here, it is helpful to treat the dynamics of strongly selected mutations separately. To do so, I define a large-effect variant as having Ng 2 $ 100, meaning that the deterministic force of selection is much stronger than that of drift. To vary the probability that a new mutation is of large effect, I performed a second set of simulations, also involving 10 unlinked loci, varying the distribution of effect sizes (DES) such that the probability was that Ng 2 $ 100 would take on values of 0.1, 0.5, or 0.75. For Gaussian DES, the mean g is zero, as above, and s g is found by numerical optimization using scipy (Jones et al. 2001) to give the desired PrðNg 2 $ 100Þ. I also used g distributions with shape parameters equal to either one or one-half, and then found a value for the mean of the distribution using scipy. These shape parameters gave probability density functions that were "exponential-like" in shape. For simulations with g DES, I used an equal mixture of g distributions with mean g and 2g such that the DES was symmetric around a value of zero. I performed 100 simulation replicates for each parameter combination. Using the argument from above, assuming hypothetical simulations of a trait with a heritability of one-half, the Gaussian distribution and the g with a shape of one gave a ratio of the mutational variance to the environmental variance of 2 3 10 23 to 3 3 10 23 when the proportion of new mutations with Ng 2 $ 100 was 0.1. These values were close to the mean of 10 23 reported for a variety of traits [Lynch (1988) and Falconer and Mackay (1996), p. 349].
In a third set of simulations, I varied r ¼ 4Nr, the recombination rate within each locus. I ran 256 replicates of these simulations using the tree sequence recording algorithm (Kelleher et al. 2018) implemented in fwdpy11 version 0.3.2. For these simulations, I recorded the entire population as nodes in the tree sequences for each of 200 generations after the optimum shift. Recording nodes from these time points allows them to be analyzed after the simulation has completed. Each replicate was simulated twice. The first run simply output metadata about mutations that reached fixation. The second run was performed with the same random number seed as the first and used the metadata from the first run to track linkage disequilibrium around fixations over time, outputting those data along with the tree sequence for the simulation.

Genome scan statistics from multilocus simulations
The 10-locus simulations described above were used to look at the temporal dynamics of several population-genetic summaries of a sample. Each of the 10 loci consisted of 11 nonoverlapping windows ( Figure 2) and all summary statistics were calculated on a per-window basis. I used pylibseq version 0.2.1 (https://github.com/molpopgen/pylibseq), which is a Python interface to libsequence (Thornton 2003), to calculate all genome-scan statistics. All statistics were obtained from 50 randomly chosen diploids.
z-scores for the nS L statistic: Individual values of the nS L statistic (Ferrer-Admetlla et al. 2014) from the first and last window of each locus were binned into intervals of size 0.1 based on derived frequency. These windows were used because they were the furthest from mutations affecting trait values, and thus the least affected by linked selection. The data from all loci were combined, and the means and SDs of each bin were used to obtain z-scores for markers from the remaining windows. Software availability: I used fwdpy version 0.0.4 (http:// molpopgen.github.io/fwdpy) compiled against fwdpp version 0.5.4 (http://molpopgen.github.io/fwdpp) for singleregion simulations. I used fwdpy11 versions 0.1.4, 0.2.1, 0.3.2, and 0.5.1 (http://molpopgen.github.io/fwdpy11) for all multiregion simulations. fwdpy11 is also based on fwdpp, and includes that library's source code for ease of installation. Both packages were developed for the current work, but only the latter will be maintained.
All of these packages are available under the terms of the GNU Public License from http://www.github.com/molpopgen. The specific software versions used here are available for Linux via Bioconda (Grüning et al. 2017), with the exception of fwdpy11 0.2.1, which must be installed from source. I have made all Python and R (R Core Team 2016) scripts for this work available at http://github.com/molpopgen/qtrait_paper.  (Hunter 2007;VanderPlas 2016), and seaborn (http:// seaborn.pydata.org). The sqlite3 library (www.sqlite.org) facilitated data exchange between Python and R via the pandas and dplyr libraries, respectively.

Data availability
The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article. Supplemental material available at figshare: https://figshare.com/articles/simaterial_pdf/10046279.

Single-region results
In this section, I describe simulations of a large contiguous region with mutations affecting the trait occurring uniformly throughout the region. The technical details of the simulation parameters are given in the Materials and Methods. Briefly, I evolved populations for 10N generations to mutation-selection equilibrium around an optimum trait value of z o ¼ 0, at which point z o was changed to 0.1, 0.5, or 1.0 and evolution continued for another 10N generations. These simulations may be viewed as similar to the numerical calculations in de Vladar and Barton (2014) and Jain and Stephan (2017b), but with loose linkage between selected variants, whereas the previous studies assumed linkage equilibrium and I allowed for new mutation after the optimum shift. They differ from the approach of Höllinger et al. (2019) in that I simulated continuous traits and did not stop evolution once a specific mean fitness was first reached.
The mean trait value, z, rapidly approached the new optimum, typically reaching the new optimum within 100 generations [ Figure 3A, see also de Vladar and Barton (2014), Jain and Stephan (2017b), and Höllinger et al. (2019)]. Prior to the optimum shift, the average genetic variance was given by 4mV S [Turelli (1984) and Figure 3B]. Following the optimum shift, the genetic variance spiked as the population adapted [see also de Vladar and Barton (2014) and Jain and Stephan (2017b)], and then recovered to a value near 4mV S within 200 generations when the mutation rate was small and took longer to return to equilibrium when the mutation rate was higher. Figure 4 shows examples of the dynamics of z, V G , and of mutation frequencies following the optimum shift. Each example is a single simulation replicate. The top row of plots Figure 1 Schematic of the model. A Wright-Fisher population evolves to equilibrium around an optimum trait under Gaussian stabilizing selection with mean zero, where the parameter V S represents the intensity of selection against extreme trait values ðw ¼ e 2z 2 =2VS Þ. At equilibrium, the mean trait value is z 0 and the genetic variance V G equals the phenotypic variance V z . Mutations arise at a constant rate with effect sizes, g, drawn from a Gaussian distribution with mean zero and variance s 2 g . The optimum then shifts to z o . 0, such that w ¼ e 2ðz2zoÞ 2 =2VS . During adaptation, z approaches z o due to allele frequency change and new mutations. At any point during adaptation, mutations with effect sizes g . ðz o 2 zÞ=2 will overshoot the optimum if they reach high frequency or fix. shows that z quickly reached z o for the individual replicates. The approach of z to z o corresponded with a substantial increase in the genetic variance, similar to what is shown for the average genetic variance over time in Figure 3B. The middle row of panels in Figure 4 shows the frequency dynamics of mutations that eventually fixed. Importantly, z typically reached z o before the first fixation had occurred (see Supplemental Material, Figure S1 for details over a shorter timescale). The legends the panels in Figure 4 contain the effect sizes of variants where Ng 2 $ 100. The legends also contain the origin times, o, of these large-effect mutations, measured as generations since the optimum shift.
For these examples, mutations with large effects on trait value fix first, as predicted by Robertson (1956). In Figure 4, fixations of large effect typically have origin times close to zero, meaning that the mutations arose close to the time of the optimum shift. This observation is expected as such mutations contribute significantly to genetic load, and thus their equilibrium frequencies prior to the optimum shift should be near the boundaries (de Vladar and Barton 2014;Stephan 2015, 2017b). Here, because of the one-way mutation model, such large-effect variants are at frequencies near zero.
The final row of plots in Figure 4 shows the dynamics of mutations that reached a frequency of $ 1% but were eventually lost from the population. Large-effect mutations only exist for a relatively brief period of time after the optimum shift, after which most segregating variation reaching appreciable derived allele frequencies are of relatively small effect. An important observation in the final row of Figure 4 is that, for a short time following the optimum shift, several intermediate-frequency mutations with large effects on trait values may be segregating. Many of these variants are adaptive ðg . 0Þ but will only make short-term contributions to adaptation prior to their loss. The dynamics of these mutations recapitulate results from de Vladar and Barton (2014): due to epistatic effects on fitness, some mutations that are initially beneficial later become deleterious and are removed. Figure  S1 shows the data from Figure 4 over a shorter timescale, allowing a more detailed look at the dynamics of mutations during adaptation. Figure 4 suggests that fixation times are rather long, in the order of N generations even for mutations with large Ng 2 . These long fixation times are in fact typical, and large-effect mutations typically fix in N=2 to N generations ( Figure S2), which is long relative to the deterministic expectation for strongly selected sweeps from new mutations (Stephan et al. 1992). Large-effect mutations that reach fixation arise close to the time of the optimum shift ( Figure S3) and typically show shorter fixation times ( Figure S4). In general, the numbers of sweeps from new mutations and from standing variants are similar, although fixations of smaller-effect standing variants are more common in simulations with higher m ( Figure S5). In Figure S5, a sweep from a new mutation is defined as a mutation arising within 100 generations of the optimum shift and then reaching fixation. While somewhat arbitrary, this definition is justified by the rapid mean time to adaptation (Figure 3). In this model, largeeffect standing variants that fixed after the optimum shift were rare at the time of the shift ( Figure S6). Small-effect mutations were also typically rare at mutation-selection balance, in particular when the mutation rate was small ( Figure  S6).
For the parameters simulated here, and for the genetic map simulated here (Figure 2), Figures S3, S4, and S6 suggest that large-effect fixations occur from both new mutations and from standing variation, with more large-effect fixations occurring when m is smaller and/or the optimum shift is larger. Thus, we may predict that large-effect fixations from new mutations may show signs of "hard sweeps," such as an excess of high-frequency-derived neutral variants (Fay and Wu 2000;Zeng et al. 2006). Given that large-effect fixations from standing variation are typically rare at the onset of directional selection ( Figure S6), we may also expect them to affect linked neutral variation (Przeworski et al. 2005;Berg and Coop 2015). For the parameters simulated here, fixations from variants that are common at the time of the optimum shift have small effects on trait values ( Figure S6). The fixation of such mutations are unlikely to generate the patterns of haplotype diversity associated with "soft sweeps" because such patterns require strong selection on mutations at intermediate frequencies (Garud et al. 2015).

Fitness effects of mutations during adaptation
In this section, I explore in more detail the strength of selection on individual mutations during the directional phase of selection. These dynamics are relevant to the long fixation times noted in the previous section and also to the extent to which hitchhiking will affect patterns of linked variation, which is the topic of the next section. As the focus of the remaining sections will be on patterns of variation during adaptation, we switch from simulating a single large region to simulating 10 unlinked regions. The only difference between these simulations and those described above is the genetic map, and the position of mutations affecting trait values (see the Materials and Methods for technical details). Figure 5 plots the dynamics of mutations in a 10-locus system for one replicate of each of the three mutation rates used here. In each column, the gray vertical line is the time the population first reaches a mean trait value of 0:9z o , which corresponds to a mean fitness of $ 0.9 for each replicate. For simplicity, we will call this the time of adaptation. The top row of Figure 5A shows the frequency trajectories of mutations that eventually fixed. These replicates were chosen because each had one fixation of a strongly selected mutation Figure 4 Trajectories of selected mutations. The three columns show results from a single simulation replicate for low-, moderate-, and high-mutation rate simulations. Parameter values are at the top of each column. The first row of plots shows the trait value and the genetic variance (multiplied by a constant for plotting purposes) over time, for up to N generations post optimum shift. As in Figure 3, the populations adapt quickly to the new optimum of z o ¼ 1. The middle row shows the frequency trajectories of fixations. Solid, darker (purple/blue) colors reflect larger effects on trait values, and more transparent colors in the green/yellow range reflect smaller effect sizes. Fixations with effect sizes Ng 2 $ 100 are indicated in the legend. The bottom row shows the frequency trajectories of mutations that are eventually lost. The coloration is as for the fixations, and any mutations that did not reach a frequency of 1% are excluded. A maximum of five mutations, corresponding to the five largest jgj, are included in the legend in the final row. Figure S1 shows the same data on a smaller timescale, showing the details on allele frequency change during the rapid adaptation to the new optimum.
with a similar effect size. As the mutation rate increases, the genetic background of these fixing variants becomes more polygenic. As a result, the initial rate of frequency change of the fixation lessens because other mutations are involved in the response to the optimum shift, some of which may contribute to adaptation but not fix in the long-term. For all replicates, the fixations are at different loci (separated by $ 50 cM) with one exception. For the high-mutation rate case, the locus with the large-effect fixation also fixed one mutation with small g. Figure 5B shows the frequency dynamics of mutations arising prior to adaptation that were eventually lost. As the mutation rate increases, there are more large-effect mutations increasing in frequency during adaptation. For the lowest mutation rate simulated here, two such mutations are decreasing in frequency prior to adaptation. At The second row in Figure 5 shows the mean deviation of a genotype with a given mutation standardized by the SD in trait values (the z-score). The mutations that fix ( Figure 5A) are all initially found in heterozygous genotypes with trait values multiple SDs greater than the mean. Such mutations are not necessarily the largest-effect variants present at the time of the optimum shift, which is seen for the two higher mutation rates in Figure 5. The mutations that did eventually fix were initially at higher frequencies and/or associated with higher-fitness genotypes than large-effect mutations that were eventually lost.
As the population adapts, the deviation in trait value (from the population mean) for a mutation with a given effect size decreases. These z-scores decrease because the genetic variance transiently increases following the optimum shift ( Figure 3B) (de Vladar and Barton (2014); Jain and Stephan (2017b) because mutations are increasing in frequency and the variance is a function of allele frequency times the squared effect size. Mutations causing larger deviations are expected to become lost, as seen most clearly in the first column of Figure 5B: the blue and green mutations over-and undershoot the optimum, respectively. At For mutations with Ng 2 $ 10, the legend shows Ng 2 , and the mutation's origin and fixation times in parentheses, scaled so that zero is the time of the optimum shift. Defining an a 2 genotype to be any genotype containing at least one copy of these "focal" mutations, the second row shows the mean deviation from the mean trait value for the focal genotypes, standardized by the phenotypic SD. The final row shows the mean relative deviation in fitness for a 2 genotypes. The horizontal line in the last row is placed at the reciprocal of the population size ð1=NÞ. (B) The dynamics of losses. Plotting is identical to (A), but the data are filtered to only include mutations arising prior to the population first crossing 0:9z o and eventually reaching a frequency of $ 0:05. low mutation rates, there is a tendency to slightly overshoot on average (Figure 3) because such mutations will have larger initial increases in allele frequency than smaller-effect variants.
Finally, we can turn to the long fixation times. These are, in part, due to the decreasing strength of selection on individual mutations during the time period where directional selection occurs. The final row of Figure 5 shows the relative deviation due to genotypes carrying each mutation over time. As expected, genotypes with fitness above the mean increase in frequency, and these genotypes are associated with trait values multiple SDs closer to the new optimum. As the mean trait value approaches the new optimum, the relative excess fitness of these genotypes declines, approaching the reciprocal of the population size. Once the population has adapted, these mutations have small effects on phenotypic variation and their long-term dynamics are governed by underdominant selection against phenotypic variance (Robertson 1956;Kimura 1981). The underdominant selection means that mutations with frequencies greater than one-half will be weakly favored and are expected to fix, and those with frequency less than one-half will most likely be removed from the population. The small fitness differences among genotypes at the time of adaptation predict that fixation times will be slow due to relatively weak selection ( Figure 5). Note that all of the sweeping alleles in Figure 5 are from standing variation (origin times , 0) and are rare at the onset of directional selection (also see Figure S6).
Finally, traits with higher mutation rates have larger numbers of small-effect mutations segregating prior to adaptation ( Figure 5). Once the population is adapted, the deviations from mean fitness tend to be small for most genotypes and the large-effect mutants are not yet fixed, implying that interference (Hill and Robertson 1966) may also increase fixation times when the mutation rate is higher. We will return to the role of interference below. The observation in Figure 4 and Figure  Dynamics of linked selection in a multilocus system I now describe the temporal dynamics of genetic variation over time in a 10-locus system. The technical details of the simulations are identical to the previous section, and are described in detail in the Materials and Methods. Figure 6 summarizes patterns of variation in the central window ( Figure 2) of each locus where large-effect mutations segregate during adaptation to the new optimum. The figure is based on the data from Figure 5. The first two rows plot the frequency trajectories of eventual fixations and losses, and the next three rows summarize patterns of variation calculated from a random sample of individuals. These summaries of variation only show deviations from equilibrium values consistent with positive selection at loci where largeeffect fixations occur. Further, the deviations are more pronounced when the mutation rate is smaller. The partial sweeps occurring at intermediate mutation rates (middle column of Figure 6) are not associated with strong signals of hitchhiking, at least when the sample size is relatively small, as is the case here. The time when a given statistic shows its maximum departure from equilibrium values differs for each statistic and, for the replicate with m ¼ 0:001, the maximum departure may occur 100 generations after the time to adaptation. However, visually one could argue that haplotype diversity tends to minimize closer to the time to adaptation than the summaries of the site frequency spectrum. Figure 7 shows patterns of variation along each of the 10 loci from an additional simulated replicate for each of the parameters shown in Figure 5 and Figure 6. Each line corresponds to a different time point in the approach to the new optimum value of z o ¼ 1, showing data for the first time the population mean trait value crosses the thresholds of z $ 0:1, $ 0:5, and $ 0:9. While the values are noisy along a genome, it is apparent that directional selection is affecting patterns of variation at linked sites in the replicates with smaller mutation rates. In the leftmost column, where m ¼ 2:5 3 10 24 , an excess of high-frequency-derived variants is seen at locus 4, along with a reduction in haplotype diversity. A standing variant of large effect swept to high frequency at this locus during adaptation. In the middle column ðm ¼ 10 23 Þ, one sees a less-dramatic reduction in haplotype Figure 6 Signals of directional selection in single replicates of a 10-locus system. The data shown are based on the same simulations as in Figure 5. The first two rows show frequency trajectories for fixations and losses, with the colors indicating the locus where the mutation is found. The vertical gray line is the generation when the mean trait value first crosses 90% of the optimal trait value. The remaining rows show Tajima's D Tajima (1989), H9 (Zeng et al. 2006), and haplotype diversity in a random sample of 50 diploids, calculated using genotypes taken from the central "window" of a locus where causal mutations are occurring (Figure 2). diversity at locus 10, where a strongly selected standing variant reached high frequency. For these two replicates, there is some evidence of reduced haplotype diversity at loci 8 and 5, respectively, that is not associated with any fixations. In the final column, where m ¼ 5 3 10 23 , there are no obvious temporal nor spatial patterns to variation in diversity levels, and the largest deviations from the background are not associated with the fixation of beneficial mutations.
Overall, Figure 6 and Figure 7 suggest that patterns of strong hitchhiking are more likely at loci where large-effect mutations fix. Moreover, such mutations must arise on average before the mean time to adaptation. Below, when looking at average patterns of variation over time and along genomes, we will distinguish patterns of variation where fixations meeting these conditions occur from the mean pattern expected from a randomly chosen locus.
The site-frequency spectrum over time The expected histogram of mutation frequencies in a sample (the site-frequency spectrum) is a geometrically decreasing function of increasing mutation frequency under the standard neutral model (Wakeley 2008). Departures from this expectation are often summarized as single numbers whose expectations are 0 under this null model. In this section, I describe the average dynamics of two widely used statistics (Tajima 1989;Zeng et al. 2006) as a function of both time since the optimum shift and of distance from trait-affecting mutations. Figure 8 shows the average behavior of Tajima's D (Tajima 1989) over time. Figure 8A shows the mean D per window, averaging across loci and across replicates. Prior to the optimum shift, the mean D is negative in the central window containing selected variants. For highly polygenic traits, the equilibrium D is 2 0:1 in this window due to a large number of rare deleterious alleles segregating. After the optimum shift, D becomes more negative when the optimum shift is large and the mutation rate is smaller. In linked windows, the magnitude of the change in averaged D decays rapidly with increasing genetic distance.
Averaging over loci experiencing large-effect fixations, Figure 8B shows a stronger hitchhiking pattern centered on the window containing selected variants. Although the deviation in D from equilibrium decays relatively quickly both along a chromosome and over time, large-effect substitutions generate sufficiently negative D values that such loci will be enriched in the tails of empirical distributions of the statistic. Qualitatively similar patterns hold for the overall reduction in diversity ( Figure S7) and the H9 statistic ( Figure S8). The latter statistic returns to equilibrium rather rapidly, consistent with previous results (Przeworski 2002).
Here, large-effect fixations from new mutations and from standing variants have similar average effects on statistics like D and H9 (Figure 8 and Figure S8). Figure 9 shows the number of haplotypes at a locus associated with sweeps from standing variation as a function of the effect size of the variant. Here, a haplotype is defined as a unique genotype at a locus, including all neutral and nonneutral variants. Large-effect sweeps from standing variation are either extremely rare (at high m) or are rare at the time of the optimum shift when m is small, and are usually associated with few (and often only one) haplotypes at the onset of directional selection (Figure 9).
Power to reject the null model using the site-frequency spectrum Figure S9A shows the power to detect a value of D more negative than expected under the standard neutral model, Figure 7 Patterns of genetic variation along genomes in a 10-locus system during adaptation to an optimum value of z o ¼ 1 and s g ¼ 0:25. The mutation rate shown is the sum over loci and individual loci mutate at equal rates ðm=10Þ. Each column corresponds to a single simulated replicate with the mutation rate given at the top. The three rows correspond to nucleotide diversity (p), H9 (Zeng et al. 2006), and haplotype diversity. The three point colors refer to statistics calculated from 50 randomly chosen diploids in the first generation that the population mean trait value first crossed values of at least 0.1, 0.5, or 0.9. The gray shades refer to the locations within each locus where mutations affecting trait values occur ( Figure 2). The triangles along the top of each panel show where fixations occurred. Triangles pointing up are fixations from standing variation. Magenta refers to fixations with scaled effect sizes Ng 2 $ 100 and yellow refers to 1 # Ng 2 # 10.
after applying a multiple testing correction such that the per-window rejection rate under the null model is 0.05. The overall power of the test is low due to the number of tests performed (one per window) and is consistent with previous work (Braverman et al. 1995;Przeworski 2002). However, the set of loci representing "significant" deviations from the null model are enriched for large-effect substitutions ( Figure S9B), of which there are relatively few per replicate ( Figure S10). When mutation rates are smaller, significant D values are most common at loci where large-effect mutations fix. As the trait becomes more polygenic and/or the optimum shift is less drastic, the enrichment shifts toward sweeps from standing variation.
The behavior of H9 is similar to that of D, but power decreases more rapidly with time since the optimum shift [Figure S11A; also see Przeworski (2002)]. The behavior of a related test, the composite likelihood ratio test of Nielsen et al. (2005), evaluated using SweeD (Pavlidis et al. 2013), is qualitatively similar to that of H9 ( Figure S12).

Haplotype homozygosity
Rapid increases in allele frequency due to selection will result in long stretches of homozygosity flanking the selected mutation (Kim and Nielsen 2004). Summaries of haplotype homozygosity are widely used to detect recent selection (Voight et al. 2006;Ferrer-Admetlla et al. 2014) and are indirect summaries of the underlying linkage disequilibrium in the data (Sabatti and Risch 2002).
The nS L statistic (Ferrer-Admetlla et al. 2014) measures the ratio of homozygosity on the ancestral allele to that on the derived allele for each variant in the data. A negative value of the statistic implies longer runs of homozygosity around the derived allele. Figure S13 shows the average behavior of z-scores obtained from binning nS L scores by derived allele frequencies (see the Materials and METHODS). The signal of strong positive selection, indicated by a negative z-score, is shortlived, and only observed when the mutation rate is smaller and the optimum shift is large. The signal is also restricted to regions closest to where selected mutations arise.
Shortly after the optimum shift, the mean z-score becomes positive ( Figure S13). This temporal dynamic is qualitatively similar to what is seen under standard models of selective sweeps, as the time since the sweep moves further into the past ( Figure S14). Thus, the positive z-scores in Figure S13 may be interpreted as either older sweeps from new mutations or strong sweeps from common variants. However, the latter class of sweeps does not occur in these simulations ( Figure  9). This difficulty in interpretation is a general issue arising from the fact that patterns of variation due to strong sweeps from standing variation overlap considerably with those of older sweeps from new mutations (Schrider et al. 2015).
A related class of statistics designed to detect strong sweeps from standing variation are based on the overall haplotype diversity in a window (Garud et al. 2015). The temporal patterns associated with these statistics are again short-lived and are all in the direction of reduced overall haplotype heterozygosity, which is a signal of strong sweeps from new mutations (Figures S15, S16, and S17).

Robustness to variation in the recombination rate
In this section, I explore the effect of varying the scaled recombination rate within a locus, r. At higher mutation rates, longer fixation times are more likely as r decreases ( Figure 10). In individual replicates, there is a tendency toward negative disequilibrium among beneficial mutations (g . 0, Figure S18), suggesting a role for interference among selected sites affecting times to fixation (Hill and Robertson 1966;Felsenstein 1974). In the previous sections, the ratio of r to u within loci was one, which is roughly "human"-like (Dumont and Payseur 2008;Ségurel et al. 2014). For species like Drosophila melanogaster, where r u (Haddrill et al. 2005), fixation times will be much shorter on average ( Figure  10). Note that the effect of recombination rate on fixation Figure 9 The number of haplotypes associated with fixations from standing variation of different effect sizes. Each panel shows the effect size of a fixation from standing variation (x-axis) and the number of unique haplotypes in the entire population containing that mutation. The number of haplotypes for each mutation is taken immediately prior to the optimum shift and excludes any mutations that arose that generation. Thus, all mutations found on a single haplotype are more than one generation old.
time is most dramatic for m ¼ 0:005, which is also the part of the parameter space explored here where fixations of largereffect ðNg 2 $ 1000Þ are rare.
The within-locus recombination rate has no discernible average effect on z nor on V G ( Figure S19). The differences in the height of the "spike" in V G when m ¼ 2:5 3 10 24 show no clear pattern with r and are thus attributable to Monte Carlo error in estimating a second-order statistic from 256 replicates.
Unlike the mean trait value and variance, the mean temporal dynamics of summaries of variation data are strongly affected by r ( Figure S20) as expected (Kaplan et al. 1989;Braverman et al. 1995). Figure S21 shows how the withinlocus recombination rate affects patterns of haplotype diversity in a 10-locus system with s g ¼ 0:25. When r is small, the impact of linked selection is much more apparent. These effects of the local recombination rate on patterns of hitchhiking are expected from standard theory of directional selection, because both the magnitude and extent along the genome of linked selection depend on the ratio of the recombination rate to the selection coefficient (Kaplan et al. 1989;Durrett and Schweinsberg 2004;Nielsen et al. 2005).

Varying the DES
The results described in the previous sections are based on a Gaussian DES whose SD is held constant. In this section, I vary the DES such that the fraction of mutations with Ng 2 $ 100 varies, and compare the average dynamics of adaptation and patterns of hitchhiking. I also compare a Gaussian distribution to a g distribution with different shape parameters. To simplify the presentation, I only show results for the case of a large optimum shift ðz o ¼ 1Þ, which is the case resulting in the most extreme hitchhiking signals. I compare the results of Gaussian distributions of effect sizes to two g distributions with shape parameters of one and one-half.
Varying the fraction of large-effect mutations has a weak effect on the mean time to reach the new optimum, with traits with low mutation rates adapting more slowly on average when the majority of variants are of small effect ( Figure S22A). This observation should be unsurprising as the population must wait longer for a strongly selected mutation in this case.
Patterns of variation expected due to hitchhiking are more extreme when PrðNg 2 $ 100Þ is small, as the population has to wait longer for strongly selected variants ( Figure S23). The overall pattern is that the average differences between DES are subtle, with g distributions showing less-extreme hitchhiking patterns (negative values) on average than the Gaussian DES. However, this difference between DES is only observed when both the mutation rate and the proportion of new mutations of large effect are both small.

Discussion
I have used simulations to describe the average behavior of selected and neutral mutations during the adaptation of a Figure 10 The sojourn times of fixations in a 10-locus system with varying recombination rates within loci. The x-axis is the effect size of a fixation and the y-axis is its time to fixation scaled by the population size. The expected fixation time due to drift alone is 4 and the distribution of fixation times under neutrality has a long tail including large values. The points are collected from 256 replicates for each parameter combination and colored by order-of-magnitude ranges of their scaled strength of selection ðNg 2 Þ. quantitative trait to a single, sudden shift in the optimal trait value. The genotype-to-phenotype model considered here is the classic model of evolutionary quantitative genetics, assuming strictly additive mutational effects on trait values with fitness determined by Gaussian stabilizing selection (Turelli, 1984;Barton 1986;Bürger 2000). The primary goal here was to merge this model of a phenotype with the simulation methods commonly used in population genetics to study the effect of natural selection on the dynamics of linked neutral variation (Kaplan et al. 1988(Kaplan et al. , 1989Braverman et al. 1995;Przeworski 2002;Innan and Kim 2004).
The simulations performed here have several important differences from recent theoretical treatments of adaptation to sudden optimum shifts (see below). However, the conditions for a selective sweep are consistent with predictions made using theoretical results from Jain and Stephan (2017b) and Höllinger et al. (2019). Direct comparison with the quantitative predictions from Jain and Stephan (2017b) is difficult because their expressions depend on the assumption of equal forward and backward mutation rates at each position. However, several qualitative comparisons can be made. First, the simulations presented here are comparable to the "most effects are large" case from Jain and Stephan (2017b) because the trait variance increases during adaptation [also see de Vladar and Barton (2014)] due to largeeffect mutations moving from low to intermediate frequency.
Mutations with large effects on trait values at the time of the optimum shift are most likely to rise in frequency (Figure 4 and Figure 5), although mutations that eventually fix are not necessarily those with the largest effect size. When several large-effect mutations cosegregate, those with the highest initial frequencies tend to reach fixation. If initial frequencies are similar, the variant with the highest initial fitness typically fixes. For a given DES, faster sweeps are more likely at lower mutation rates.
Regimes where the genetic variance decreases during adaptation are not possible for any of the simulations presented here. The decrease in variance is seen in the "most effects are small" domain where the equilibrium frequency of variants prior to the optimum shift is one-half, which maximizes the variance (de Vladar and Barton 2014;Jain and Stephan 2017b). Adaptation to the new optimum displaces allele frequencies, reducing the variance from its maximum value [see, for example, figure 9 of de Vladar and Barton (2014)]. However, the equilibrium frequency of one-half for small-effect mutations requires equal rates of forward and back mutation (de Vladar and Barton 2014; Jain and Stephan 2017b, and is therefore incompatible with the infinitely many sites assumption made here. When considering the pattern of hitchhiking at a locus, the presence or absence of a large-effect fixation at a locus is a reliable predictor of the magnitude of hitchhiking patterns. As expected, such fixations are more common when the mutation rate is smaller (Höllinger et al. 2019) and thus strong departures from equilibrium patterns of variation are not expected for more polygenic traits (Figure 8). For the optimum shift model considered here, the strength of selection is not constant over time [ Figure 5; see also Kimura (1981)]. Thus, genotypes containing variants that were initially strongly favored by selection are subject to much weaker selection by the time the population has reached the new optimum. This weakening of selection increases fixation times to the order of the population size ( Figure S2), which is much longer than the times N generations expected for directional selection in large populations (Stephan et al. 1992).
The exploration of hitchhiking signals here involved the simulation of 10 unlinked loci within which mutations affecting the trait were concentrated in a central window (Figure 2). While the ratio of recombination to mutation events is at least nine to one for the majority of the results shown here (see Materials and Methods), it is possible that signals of selection are made more pronounced by the localization of selected mutations and should be explored further.
Here, the number of selected mutations segregating over time ranged from dozens to several hundred, as a function of the underlying mutation rate ( Figure S24). At high mutation rates, the number of segregating loci are roughly the same as some of the results presented in de Vladar and Barton (2014). However, the partial linkage among sites in this work leads to some negative linkage disequilibrium ( Figure S18), which is a signal of interference (Hill and Robertson 1966;Felsenstein 1974). This interference has little effect on the mean time to adaptation, but fixation times are increased. The lack of effect on time to adaptation is driven by initial large fitness differences among genotypes [ Figure 5, also see Höllinger et al. (2019)]. Once the population is close to the new optimum, selection on individual genotypes is much weaker (Figure 5), setting up the conditions for interference to affect fixation times (Hill and Robertson 1966).
The DES has different effects on properties of the trait than on patterns of hitchhiking. The mutation rates used here span the parameter space from partial and complete sweeps being most common to the optimum being reached via allele frequency shifts of many mutations (Figure 4; Höllinger et al. 2019). In general, the mean time to adapt is not strongly affected by the DES if the fraction of new mutations of large effect is constant ( Figure S22A). For a given mutation rate, lowering the mutational variance lowers the probability of a strongly selected mutation, increasing the waiting time until such mutations arise, and thus resulting in stronger signals of hitchhiking Figure S23). When the trait is more polygenic, the average patterns of variation are not strongly dependent on the DES nor on the proportion of new variants with large effect ( Figure S23).
The genetic model assumed here does not lead to sweeps of large-effect mutations from common variants (frequencies greater than, say, 5%). Rather, the stabilizing selection around the initial optimum keeps large-effect mutations rare, such that sweeps from such standing variants start at low frequencies. Importantly, it is not possible to tune the model parameters to obtain sweeps from large-effect, but common, variants with high probability. Changing the strength of stabilizing selection ðV S Þ preserves the rank orders of fitness for all genotypes, merely changing how fit they are in an absolute sense. One could randomly reassign effect sizes at the time of the optimum shift in an attempt to approximate a gene-byenvironment interaction. However, such a procedure would be arbitrary, and thus not represent a principled model for generating detectable soft sweep patterns. Rather, it is tempting to invoke a need for pleiotropic effects to have large-effect mutations segregating at intermediate frequencies at the time of the optimum shift, with the shift itself accompanied by a change in the covariance between trait values and fitness.
It is important to note a key methodological difference between this work and that of other authors. Höllinger et al.  2019) show is the parameter range where adaptation occurs primarily by changes in allele frequency. These results are consistent with the theoretical predictions from Höllinger et al. (2019), as the fixations in the simulations described here take place on timescales longer than the mean time to reach the new optimum. In the rightmost column of Figure 4, the population has adapted quickly, with the fixations occurring over a much longer timescale ( Figure S1). Likewise, the leftmost column of Figure 4 corresponds to Q ¼ 5, where we observe a mixture of partial and complete selective sweeps by the time the new optimum is reached, which is expected from the theory presented in Höllinger et al. (2019).
This work [and that of Höllinger et al. (2019)] differs from the analytical and numerical work of de Vladar and Barton (2014) and Stephan (2015, 2017b) in several key aspects. First, we consider irreversible mutation here [the infinitely many sites model of Kimura (1969)], while de Vladar and Barton (2014) assumed equal rates of forward and reverse mutation [see also Barton (1986) and Stephan (2015, 2017b)]. The infinitely many sites model used here was chosen because it is the most commonly used mutational model for investigating the effects of linked selection during adaptation (e.g., Braverman et al. 1995;Przeworski 2002;Przeworski et al. 2005). I also allowed for partial linkage among sites, which is a key difference from the work based on the Barton (1986) framework, which assumes free recombination. As noted above, partial linkage affects the long-term dynamics of selected mutations ( Figure 10).
I have focused on standard summaries of variation data that have been widely applied to detect selection from sequence data. The behaviors of the majority of such summary statistics have only been tested using coalescent simulations of strong selection on a single sweeping variant, which is the dominant generative model used to make predictions about linked selection. Thus, it is unsurprising that these statistics show the strongest departures from equilibrium neutrality for traits with low mutation rates. However, an important observation here is that the mean behaviors of these statistics are similar for sweeps from new mutations and sweeps from standing genetic variation, which is a consequence of the standing variants being rare at the onset of selection [ Figure  S6; also see Orr and Betancourt (2001), Hermisson and Pennings (2005), Przeworski et al. (2005), and Berg and Coop (2015)]. The only test statistic based on patterns of SNP variation for detecting polygenic adaptation that I am aware of is the singleton density score (Field et al. 2016). I have not explored this statistic here, as it would be more fruitful to do so using simulations of much larger genomic regions applying tree sequence recording (Kelleher et al. 2018), and explicit modeling of trait architectures at or near the infinitesimal limit Robertson (1970Robertson ( , 1977. It also appears that the magnitude of selective effects on phenotypes attributable to changes in the singleton density by Field et al. (2016) were substantially overestimated due to uncontrolled-for population structure in the genome-wide association study data, and there was little evidence for selection on height when the analysis was redone using effect sizes from the UK Biobank data (Berg et al. 2019;Sohail et al. 2019).
I have only considered the equilibrium Wright-Fisher model here. However, it is well understood that departures from this demographic model affect patterns of neutral variation and thus the detection of regions affected by linked selection (Thornton and Andolfatto 2006;Jensen et al. 2007Jensen et al. , 2008; ). Demographic departures from constant population size indeed affect the prevalence of sweeps and the rate of phenotypic adaptation in optimum shift models (Stetter et al. 2018). Here, we are primarily interested in how the parameters affecting the trait's "architecture," mainly the parameters affecting the mutational variance of the trait, impact patterns of linked selection.
It is crucial to restate the assumptions of the genetic model assumed here, which involves strictly additive effects on a single trait under real stabilizing selection (Johnson and Barton 2005). This model is the standard model of evolutionary quantitative genetics (Turelli 1984;Barton 1986;Bürger 2000), which is why it is the focus of this work. However, a more thorough understanding of the dynamics of linked selection during polygenic adaptation will require investigation of models with pleiotropic effects (e.g., Zhang and Hill 2002;Simons et al. 2018). Because the adaptation to the new optimum is rapid when the mutation rate is large, the allele frequency changes involved are also small when mutational effects are pleiotropic (Simons et al. 2018). The question in a pleiotropic model is the role that large-effect mutations may play, which is an unresolved question.
The simulations here also model the entirety of heritable variation for the trait. An alternative approach would be to allow for an unlinked additive genetic background with its own mutational variance. Such an approach would be straightforward assuming an infinitesimal model for the background, as has been done previously (Chevin and Hospital 2008;Stetter et al. 2018). Stetter et al. (2018) simulated "domestication" traits evolving to a new optimum via truncation selection and a heritable background affecting the focal trait. They concluded that the contribution of genetic background to several outcomes of interest (speed of adaptation, fixations of beneficial mutations, etc.) was of overall less importance to the dynamics than the variance in mutational effect sizes, s g . Clearly, however, the details will depend on the specifics of the model, with Chevin and Hospital (2008) at one extreme and the current work at perhaps the other. Here, the simulations with high mutation rates imply that any single segregating variant finds itself in a mutation-rich genetic background of up to several hundred segregating variants, the majority of which have small fitness effects ( Figure  S24). Another appealing alternative would be to simulate entire genomes using an adaptation of Robertson's (1977) method to incorporate tree sequence recording (Kelleher et al. 2018) and large-effect mutations occurring at some rate. Such a scheme would generate large-effect genomic regions through two different mechanisms: the occasional large-effect mutation as well as via large-effect haplotypes arising from stochastic recombination events (Sachdeva and Barton 2018).
It may also be of interest to explore nonadditive genetic models in future work. In particular, models of noncomplementing recessive effects within genes are a specific class of model with epistasis that deserve consideration due to their connection with observations of allelic heterogeneity underlying variation in complex traits (Clark 1998;Gruber and Long 2009;McClellan and King 2010;Thornton et al. 2013;King et al. 2014;Long et al. 2014;Sanjak et al. 2017;Chakraborty et al. 2018). Acknowledging the focus on the standard additive model, the current work is best viewed as an investigation of a central concern in molecular population genetics (the effect of natural selection on linked neutral variation) having replaced the standard model of that subdiscipline with the standard model of evolutionary quantitative genetics. As laid out by several authors (Messer and Petrov 2013;Jain and Stephan 2017a,b), there are considerable theoretical and empirical challenges remaining in the understanding of the genetics of rapid adaptation. For models of phenotypic adaptation, our standard "tests of selection" are likely to fail, and are highly underpowered even when the assumptions of the phenotype model are closer to that of the standard model.