## Abstract

Mutations are crucial to evolution, providing the ultimate source of variation on which natural selection acts. Due to their key role, the distribution of mutational effects on quantitative traits is a key component to any inference regarding historical selection on phenotypic traits. In this paper, we expand on a previously developed test for selection that could be conducted assuming a Gaussian mutation effect distribution by developing approaches to also incorporate any of a family of heavy-tailed Laplace distributions of mutational effects. We apply the test to detect directional natural selection on five traits along the divergence of Columbia and Landsberg lineages of *Arabidopsis thaliana*, constituting the first test for natural selection in any organism using quantitative trait locus and mutation accumulation data to quantify the intensity of directional selection on a phenotypic trait. We demonstrate that the results of the test for selection can depend on the mutation effect distribution specified. Using the distributions exhibiting the best fit to mutation accumulation data, we infer that natural directional selection caused divergence in the rosette diameter and trichome density traits of the Columbia and Landsberg lineages.

A great deal of progress has been made modeling selection on molecular characters, especially DNA sequences (Clarke *et al.* 1988; Nielsen and Yang 1998; Nielsen 2005; Yang 2006). While progress has been made in understanding the short-term effects of selection on phenotypic traits (Sorensen and Hill 1982; Skibinski and Shereif 1988; Kingsolver *et al.* 2001; Hereford *et al.* 2004; Siepielski *et al.* 2009), and methodologies for identifying the contribution of evolved genetic variation to specific traits are well developed (Liu 1998; Lynch and Walsh 1998; Doerge 2002), little progress has been made developing methods to assess the degree of selection on quantitative traits between populations whose traits have diverged over longer intervals. Studies of quantitative traits, which have been conducted to determine the genetic architecture of phenotypic differences between divergent lineages, are a possible source of data that should be informative with regard to this question. Indeed, one of the major questions arising—but often left unanswered—in these studies is the degree to which QTL evolved under neutrality or natural selection.

This dearth in answers arises because, in contrast to examination of molecular characters, there have been far fewer approaches developed to use QTL data to test the strength of selection for divergent phenotypes. Orr (1998a) pioneered usage of QTL locus effect sizes for detection of selection by developing the QTL sign test and the QTL effect size test of selection. These tests integrate QTL data and population genetic theory and use a null model of neutral phenotypic evolution (Orr 1998a). However, the QTL sign test has a high false-positive rate (Anderson and Slatkin 2003) and is highly sensitive to the variance in QTL effects (Rice and Townsend 2012b). In particular, it has negligible power for datasets with low QTL effect variance (Rice and Townsend 2012b). The QTL effect size test also provides very low power in detecting the presence of selection (Anderson and Slatkin 2003). Neither of these two tests explicitly models selection (Orr 1998b). Riedel *et al.* (2015) developed a population genetic test for selection using multiple-line crosses that has higher statistical power in detecting selection, compared to the two-line test. The multiple-line test is successful in statistically distinguishing the effect of lineage-specific selection from the bias resulting from testing traits with large phenotypic differences, while the two-line test is not. However, neither the Orr (1998b) tests nor the Riedel *et al.* (2015) test incorporate any assessment of the distribution of mutational effects (Orr 1998b). Because sums of mutational effects represent the baseline change expected for a quantitative trait in the absence of natural selection, they are critically important to the attribution of phenotypic change to neutral mutations or selected mutations. Consequently, Rice and Townsend (2012a) developed an approach that determined the magnitude and direction of selection on the QTL by comparing the mutational effect distribution of a trait from a mutation accumulation study to the realized distribution of QTL via a specific population genetic model of mutation and selection. Applying their test to a difference between lineages of *Drosophila melanogaster* selected in the laboratory for higher and lower counts of sternopleural and abdominal bristles, their test assumed a Gaussian distribution of mutational effects on sternopleural and abdominal bristles. Using the estimated distribution of mutational effects and the probability of fixation, the distribution of mutational effects that go to fixation was calculated. Given a set of QTL effect sizes, the test yields a maximum likelihood estimate of the strength and direction of selection.

To characterize the distribution of mutational effects, Rice and Townsend (2012a) assumed an analytically tractable Gaussian distribution of mutational effects, while Mackenzie *et al.* (2005) assumed a gamma distribution of mutational effects. Here we extend this theory by developing approaches to conducting the test with alternate assumptions regarding the form of the distribution of mutational effects. In particular, we incorporate the Laplace distribution and extended versions of the Laplace distribution in addition to the Gaussian distribution to characterize the distribution of mutational effects. We experiment using various forms of the Laplace distribution, because of the conjectured under-recognition of the potential of the Laplace distribution to better fit data with wide tails and a sharp peak than does the Gaussian distribution (Kotz *et al.* 2012). If these properties are characteristic of the distribution of the mutational effects of a quantitative trait, then any of a number of tractable Laplace-related distributions might represent more appropriate models for the effects of mutations on phenotypes than would the Gaussian, or displaced Gaussian, distribution.

Furthermore, to apply this test to phenotypes of lineages subject to natural selection rather than artificial selection, we have assessed the level of selection required to explain divergence in several quantitative traits for which QTL were determined that differentiate two lineages of *Arabidopsis thaliana*. We quantified selection on five quantitative traits involved in life history, plant architecture and morphology: bolting time, rosette diameter, number of rosette leaves at bolting, number of elongated axils, and trichome density (Supplemental Material, Table S1 in File S1). For each of these five traits, we examined the performance of the proposed mutational effect distributions, and tested for the presence of selection using each distribution. Contemporary selection has been detected in *A. thaliana* for all of these traits in experiments in the field (Mauricio 1998; Scheiner *et al.* 2000; Callahan and Pigliucci 2002; Donohue 2002). From our results, we considered the implications of selection to the evolution of *A. thaliana*, and the evolution of quantitative traits and phenotypes.

### Theory

#### Overview:

Selection strength was inferred by integrating mutational effect distributions and mutation accumulation data (Figure 1, A and B), population genetic theory, and the empirically obtained QTL effect size distributions. In our example, mutation accumulation is shown to result in an increase in variance of the phenotype with no directional bias (Figure 1A) or result in an increase in variance of the phenotype with downward bias (Figure 1B). The mutational effect distributions were obtained by fitting the models to these phenotypic observations along the time course of the mutation accumulation experiments for each quantitative trait (Figure 1, C–J). We use the distribution of mutation effects obtained from mutation accumulation experiment(s), in which mutations are allowed to accumulate in a very small experimental population in which the effects of selection have been minimized (Halligan and Keightley 2009). The effects of mutations on the quantitative trait can be inferred from assays of the trait at regular intervals or at the end of the experiment to obtain the change in the quantitative trait. With the change in trait means and the number of generations between measurements, we estimate a Gaussian distribution, the Laplace distribution, and seven extensions of the Laplace distribution with varying assumptions.

We then obtain a pointwise product of the specified distribution with the probability that a mutation fixes—which is, in turn, a function of the strength of selection on the phenotypic effect (Figure 1, K–R). For a neutral phenotype, fixation probabilities are constant across all phenotypes for a neutral phenotype (Figure 1, K, M, O, and Q). For a phenotype under selection, fixation probabilities depend on the phenotype (Figure 1, L, N, P, and R). The constant fixation probabilities, reflecting drift processes, were one multiplicand of the pointwise product with corresponding mutation effects distributions (Figure 1, C, E, G, and I), while the fixation probabilities that depend on the phenotype were one multiplicand in the point-wise product with the corresponding mutation effect distributions (Figure 1, D, F, H, and J). The resulting pointwise products provide the distributions of phenotypic differences arising from the compositional substitutions at the QTL studied (Figure 1, S–Z). Next, we sample substitutions from these distributions, total their effects, and compare them with the summed effects from the QTL (Figure 2, A–D). The distribution of mutational effects and the selective regime affect the agreement between the sampled and observed effects. This agreement is stronger when the level of selection is appropriate for generating fixed effects where the cumulative effects are similar to the effects of the observed QTL. Analysis of the match of sums of sampled substitutions to QTL differences observed, quantified within a likelihood framework, provides a means for testing the model fit, and provides a hypothesis test of neutrality (Figure 2, E and F). One can conclude that the studied trait has been subject to selection if neutrality is rejected by a hypothesis test.

#### Estimating mutational effect distributions:

Data from a mutation accumulation experiment were used to approximate potential mutation effect distributions for each phenotype. Because no single mutational effect distribution is known to be appropriate for quantitative traits (Keightley and Lynch 2003), we evaluated nine distributions (a Gaussian distribution and a family of eight Laplace distributions) for their fit to the data. As in Turelli (1984), Sawyer *et al.* (2003), Jones *et al.* (2007), and Rice and Townsend (2012a), we first specified a Gaussian distribution,(1)with mean *μ* and SD *σ*, to characterize the mutation effect distribution. The sum of *k* independent normally distributed random variables with the same mean and variance has a Gaussian distribution,(2)with a mean *k* × *μ* and SD *k* × *σ*. Consequently, assuming that the effects of mutations are independent, the distribution of the total effect of *k* mutations is the Gaussian distribution in Equation 2.

We next fit a family of eight Laplace distributions (*cf*. Kotz *et al.* 2012) to comparatively evaluate against each other and the Gaussian distribution. We began with a Laplace distribution of mutational effects characterized by symmetric submodal and supermodal exponential decay, equally proportionate submodal and supermodal probability masses, and a mode at zero (LSPZ). Then, extensions of the Laplace distribution are developed by breaking the assumptions that the submodal and supermodal exponential decay is symmetric (LA** *vs.* LS**; Figure 1, I and J), that the distribution has equal submodal and supermodal probability masses (L*D* *vs.* L*P*; Figure 1, E and F), and that the mode is zero (L**D *vs.* L**Z; Figure 1, G and H). These additional distributions (LSPZ, LSPD, LSDZ, LSDD, LAPZ, LAPD, LADZ, and LADD) are of interest because the distribution of mutational effects could have asymmetric submodal and supermodal exponential decay, disproportionate submodal and supermodal probability masses, and/or a nonzero mode. For each version of the Laplace distribution, a different set of assumptions hold, which results in a family of eight Laplace distributions. To select among these alternate distributions, assessments of fit to the mutation accumulation data are performed.

Assuming that the distribution has symmetric submodal and supermodal exponential decay, equally proportionate submodal and supermodal probability masses, and a mode at zero, we started with a symmetrical Laplace distribution of mutational effects(3)featuring a single exponential decay parameter *b*. By Taillie *et al.* (1980), the sum of *k* Laplace distributions with exponential decay parameter *b* is distributed(4)To provide for the possibility of higher-complexity distributions of mutational effects, we then relaxed the assumptions that enable the classic Laplace distribution to be fit with a single exponential decay parameter. We created a family of eight Laplace distributions by permutatively relaxing three parameters: asymmetry of the submodal and supermodal exponential decay (*b*), asymmetry of supermodal and submodal probability masses, and displacement of the mode from zero. Tallying all permutations of these parameters as fixed or free defines the Laplace distribution and seven extensions of the Laplace distribution (Table 1, Table 2, and Table 3).

Out of these extended distributions, only the sum of *k* LSPD distributions with exponential decay parameter *b* and displacement parameter *a* has a closed form solution and is characterized by(5)The other distributions in the family of Laplace distributions do not have a sum that can be characterized by a closed form solution. However, they can be approximated by sampling a vector of mutations from each distribution *k* times, adding the mutational effects, and fitting a kernel smoothing density function on the summed effects (Bowman and Azzalini 1997).

Equations 2, 4, and 5 provide the distribution of mutational effects for *k* mutations. However, we do not know *k* when we observe QTL data. Rather, we would like to calculate the distribution of the total effect of the mutations that might accumulate in *g* generations. Assuming independence of the mutations, we apportion the likely effects by the Poisson probability of each value of *k*(6)where *u* is the per-generation rate of mutations that affect the trait. Then, the sum over *k* of the distribution of the total effect of *k* mutations is weighted by the probability of *k* mutations occurring in *g* generations. Thus, the total effect of *k* mutations weighted by the probability of *k* mutations occurring in *g* generations using the Gaussian distribution is provided by(7)and the total effect of *k* mutations weighted by the probability of *k* mutations occurring in *g* generations using the Laplace distribution is provided by(8)and the total effect of *k* mutations weighted by the probability of *k* mutations occurring in *g* generations using the Laplace distribution with the displacement parameter *a* is provided by

#### Inference of historical selection from QTL and MA data:

With the expected distribution of the total effects of *k* mutations over *g* generations under no selection, we would like to compare and contrast to the expected distribution of the total effects of *k* mutations over *g* generations in the presence of selection. To do so, we must specify how the strength of selection is related to the probability of fixation of a mutation. To quantify the likelihood that a mutation will fix, we use the probability given by the diffusion equation result of Kimura (1962):(10)where *s* is defined to be the coefficient of selection, and *N* is defined as the size of the population.

We assume the selective value of a mutation is proportional to its effect on a quantitative trait, as in Lande (1983) and Chevin and Hospital (2008).

Furthermore, we assume that the fitness gradient has not changed over the course of evolution of the QTL, freeing us from making any assumptions regarding the optimal value of the trait. Perhaps counterintuitively, this assumption of a constant fitness gradient is likely conservative (for the goal of detecting the action of selection). This conservatism arises because under a typical Fisherian adaptive walk, a few fixations of mutations of large effect are followed by a larger number of mutations of smaller effect as the optimum is approached (Fisher 1930; Orr 1998b, 2005). These potentially numerous mutations are of smaller effect than our model incorporates, based on an unchanging fitness gradient. These mutations of small effect will tend toward consistency with neutral fitness effects when evaluated under our model, and will decrease the sensitivity of our model to the signature of selection. We also assume that sequential fixation of beneficial mutations (Atwood *et al.* 1951) occurs, instead of concurrent selection at different loci, driving adaptation. This assumption is also conservative, because if selection occurs at many loci, then the effective strength of selection on an individual locus decreases over time (Chevin and Hospital 2008). Lastly, we assume that mutation only affects the trait in question, and therefore that there are no additional effects on fitness. Because pleiotropic effects are considered most often to be antagonistic to directional selection on a quantitative trait, this assumption is also conservative (Griswold and Whitlock 2003). Under these assumptions, we parameterize the selection coefficient by(11)where *x* is the phenotypic effect of a mutation, and *c* is a constant whose sign determines the direction of selection and magnitude determines the strength of selection. If *c* is positive, it indicates that an increase in the trait was selected in the focal population. A negative value of *c* indicates that a decrease in the trait was selected in the focal population. A value of *c* = 0 is consistent with an absence of directional selection of the trait over the time course of divergence, but does not in our implementation rule out any number of other scenarios, such as stabiliz-ing selection. Substituting Equation 11 into Equation 10 yields(12)which provides the probability of fixation of a mutation with size *x*.

The product of the distribution of effect sizes of mutations and the corresponding fixation probability provides the probability distribution of the effects of substitutions over the course of selection:(13)where is one of the Gaussian or Laplace distributions. The normalization constant *R* in Equation 13 is included to ensure that *W*(*x*) is a probability density that integrates to one over all *x*.

## Methods

### Overview of the data required

We used the change in trait value and the number of generations between each measurement from a mutation accumulation experiment to calculate the parameters of our nine mutational effect distributions (Rutter *et al.* 2010). To calculate the probability of fixation of mutations, we used estimates of the historical effective population size of the population. The effect sizes and SEs of QTL, obtained from crosses of individuals from divergent lineages with different trait values, were used to construct a likelihood function for the parameters of the distribution and the mutation rate for the trait.

### Approximating mutational effects distributions

We were not able to derive sums of Laplace distributions LSDZ, LSDD, LAPZ, LAPD, LADZ, and LADD in closed form. Therefore, we approximated the distributions of the summed mutational effects by randomly sampling *S* = 2.5 × 10^{4} mutational effects from the summand distribution *k* times, and fitting a smoothing function on *S* sums of *k* effects. This kernel smoothing was performed using the ksdensity function in MATLAB, which returns a probability density estimate for the sampled data.

### Likelihood function for mutational effects distributions

Using the change in trait value in each line of the mutation accumulation experiment, the moments of the mutation effect distributions were inferred. A likelihood function for the parameters of each distribution and the mutation rates was built given the change in trait means, and the number of generations between assays in the mutation accumulation experiment. The likelihood function for the Gaussian distribution is(14)and the likelihood functions for the LSPZ and LSPD distributions are(15)(16)where the function *Y*(…) is provided by Equations 7–9, and *x _{i}* and

*g*are the change in the mean value of the trait, and number of generations between each observation. The likelihood functions for the remaining versions of the Laplace distributions are similar in form to Equations 14–16. However, instead of closed form equations, approximated functions are used in the likelihood functions.

_{i}These likelihood functions are infinite sums, so we were not able to calculate them with perfect precision. However, we were able to obtain accurate estimates for small *g* × *u* by using a small number of terms normalized by the total of all the terms. We included enough terms so that the total weight of the terms was at least 0.999 for each *g* × *u*.

### Estimating SEs of the QTL effect sizes

The SEs of the observed QTL were calculated based on QTL marker and phenotypic data using a recombinant inbred line mapping population. First, the subjects were split into two groups according to their genotype at each marker. The estimated difference in phenotypic values of the two groups is the difference of the sample means of the two groups,(17)We assumed that the observations were independent, and that, based on the central limit theorem, the sample populations are approximately normally distributed. Thus, the difference between sample means would be normally distributed, and the SE of the difference between sample means can be estimated using the SEs and number of observations within each group:

(18)### Estimating selection

We used a maximum-likelihood approach to infer the strength of historical selection on the quantitative trait. Substitutions sum to yield the phenotypic effects of the QTL, so we maximize the joint likelihood of *c* and the vector *J* = (*J*_{1}, *J*_{2}, …, *J _{n}*), where

*J*is the number of substitutions at the

_{i}*i*th locus, and

*n*is the number of total loci. The number of mutations at each locus and the direction of the mutational effects are unknown and are not constrained; however, there is a constraint on the total effect provided by the QTL. Because of this, our approach is able to handle cases where multiple QTL are detected as one QTL (

*cf*. Mackay 2001; Mackay and Lyman 2005; Studer and Doebley 2011) even if the QTL are in opposite directions. The likelihood of

*c*and

*J*given each locus individually is(19)Here,

*q*and

*e*represent the effects and SEs of the observed QTL. The joint likelihood of

*c*and

*J*for each a single locus(20)where

*G*(…) is one of the Gaussian, Laplace, or extended Laplace probability density functions, and

*F*(…) is the probability density function of the sum of

*J*

_{i}substitutions under a strength of selection

*c*. The probability density function of the sum is calculated by

*J*

_{i}convolutions of W(

*x*) from Equation 13 with itself.

The integral in Equation 20 does not have a closed form solution. However, using Markov chain Monte Carlo sampling from the distribution *F*(*x*), and the average value of *G*(*x*) evaluated for each sampled effect, we can approximate the likelihood, and are able to approximate the integral. This approximation is provided by(21)Angle brackets in Equation 21 denote averaging over the samples *x _{j}*. The right-hand side of Equation 21 better approximates the true value of the integral as the number of samples used increases.

The right-hand side of Equation 21 was evaluated for a range of vectors of *J* given a fixed value of *c*. Starting with the value 1 we sequentially increased *J _{i}* until the unimodal maximum likelihood value was found for each locus

*i*. We defined the product over

*i*of the likelihoods (Equation 19) to be the likelihood of

*c*. The maximum-likelihood

*ĉ*was calculated with a manually guided dense grid of values of

*c*. Then we performed a likelihood-ratio test with a null hypothesis of

*c*= 0, and an alternate hypothesis of

*c*=

*ĉ*to test neutral evolution.

### Sources of data

We applied our method to QTL and mutation accumulation data (Figure 4) for the bolting time, rosette diameter, number of rosette leaves at bolting, number of elongated axils, and trichome density of *A. thaliana* (Table S1 in File S1). The QTL effect size and SE data were extracted from Ungerer *et al.* (2002) and Symonds *et al.* (2005) for the five traits (Figure 3). QTL mapping was performed in a recombinant inbred line (RIL) mapping population derived from crossing *A. thaliana* accessions Columbia and Landsberg *erecta*. These accessions were part of a collection propagated by Rédei (1992a) and have been grown in laboratories or stock centers since 1955. The two accessions likely represent distinct populations from Poland and Germany (Nordborg *et al.* 2005). Further evidence that the two genotypes originate from different populations is their difference at an inversion polymorphism that differs between European populations (Zapata *et al.* 2016). Epistasis was not observed in the number of rosette leaves at bolting, bolting time, rosette diameter, number of elongated axils, and the trichome density for Col × Ler (Ungerer *et al.* 2002; Symonds *et al.* 2005). We therefore we do not consider epistatic effects in our model.

Mutation accumulation data for the five traits was provided from greenhouse experiments described in Rutter *et al.* (2010). The five traits considered here were not discussed in Rutter *et al.* (2010), so their measurements are briefly described here. Rutter *et al.* (2010) obtained the 24th generation *A. thaliana* mutation accumulation lines derived from a single founder from the Columbia accession from Shaw *et al.* (2000, 2002), and propagated them for an additional generation (resulting in a grand total of 25 generations of mutation accumulation). For the number of rosette leaves at bolting, bolting time, rosette diameter at bolting, number of elongated axils, 100 out of 120 mutation accumulation lines (randomly chosen) were used by Rutter *et al.* (2010), and six lines that were direct offspring of the progenitor plant were used to represent the premutation genotype. For the trichome density trait, 50 mutation accumulation lines were used, and five lines were used to represent the premutation genotype. Seed for all measured plants were generated at the same time, and multiple maternal plants were used for each line to reduce maternal effects. Leaves were sampled for trichome measurement at bolting. Trichome density was calculated as the number of trichomes per square centimeter on the adaxial surface of rosette leaves. Rosette leaves at bolting, rosette diameter at bolting, and bolting time were recorded during the experiment. The number of elongated axils was measured when plants senesced. All traits were measured for the premutation lines and for the 25th generation (Figure 4).

### Estimating the distribution of mutational effects

To estimate the parameters for each of our distributions and the corresponding mutation rates, we calculated their maximum likelihood values. We let *x _{i}* and

*g*be the change in the trait and the number of generations between measurements for all lines. A hill-climbing algorithm was used to calculate the joint maximum likelihood estimates for the parameters and the mutation rate for the trait in all the distributions.

_{i}The likelihood for a range of selection coefficients was calculated for the QTL. We used an effective population size *N* = 2.3 × 10^{6} with the maximum likelihood estimates of all the parameters of the mutation effect distributions. The effective population size was calculated with the rate of mutation per site per generation, and the amount of genetic variation at drift-mutation equilibrium, using(22)where *θ* quantifies the genetic variation, *N* is the effective population size, and *μ* is the rate of mutation per site per generation (Wang 2005). This method assumes a simple situation with an isolated population, and a constant effective population size. The amount of genetic variation was estimated by(23)where *H* is the heterozygosity of the population at equilibrium (Wang 2005). We used *H* = 0.06 (Stenoien *et al.* 2005) and *μ* = 7 × 10^{−9} (Ossowski *et al.* 2010).

We obtained 1.0 × 10^{5} Markov chain Monte Carlo samples from the distribution of mutation effects for every value of *c* and *H _{i}*. We then conducted a likelihood ratio test with the test statistic(24)where

*D*is approximately distributed chi-square with one degree of freedom (Hudson 1971). Both the hill climbing algorithm and the sampling of mutations were carried out in MATLAB. The source code for executing these analyses is released under the GPLv3 license at https://github.com/Townsend-Lab-Yale/QTL-MA-test-for-selection.

### Data availability

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

## Results

### The distribution of mutation effects

The maximum likelihood of the parameters of the Gaussian distribution and the family of Laplace distributions of mutation effects were estimated for each trait (Tables S2–S6 in File S1). The Akaike information criterion (AIC) was then calculated for each of these models for each trait to quantify their fit to the data (Figure 5 and Table 4). The Gaussian and LSPD (Laplace with symmetric submodal and supermodal decay) distributions consistently yielded the lowest AIC values. Mutational effect distributions for these traits—and perhaps other traits as well—are best fit with these distributions. In contrast, the most highly parameterized distribution, LADD (Laplace with asymmetric submodal and supermodal exponential decay, unequal submodal and supermodal probability masses, and potentially nonzero mode) produced consistently high AIC values relative to the other distributions. Estimating free parameters for both the asymmetry of exponential decay and the asymmetry of submodal and supermodal probability mass was not necessary to achieve good fit to mutational effect data.

### Likelihood ratio tests

Likelihood ratio tests were carried out for each trait using the QTL detected and the estimated mutational effect distributions. A significance threshold of *<* 0.05 was specified to determine whether a neutral hypothesis for the change in each trait should be rejected. The corresponding chi-square statistics, *P* values, and the maximum likelihood estimates of c were calculated for each trait (Tables S7–S11 in File S1). The largest sampled likelihood was used to identify the maximum likelihood estimate of the selection *c*. The confidence intervals of the selection values were also obtained for each trait (Figure 6).

Likelihood ratio tests for rosette diameter and trichome density rejected the neutral hypotheses for the increases in the traits. All nine likelihood ratio tests rejected a neutral hypothesis for an increase in trichome density, while seven out of nine likelihood ratio tests rejected a neutral hypothesis for an increase in rosette diameter. All likelihood ratio tests except those that utilized the LSPZ and LSDD distributions rejected the neutral hypothesis for an increase in the rosette diameter. A neutral hypothesis could not be rejected for the number of elongated axils, bolting time, or the number of rosette leaves, based on any analyzed model of the mutation effect distribution.

## Discussion

We have extended both the theory and the application of a test for historical selection developed by Rice and Townsend (2012a), applying it to five traits of *A. thaliana*. We first expanded the method by fitting the Laplace distribution and a family of related distributions—in addition to the Gaussian distribution—to observed mutational effects. We demonstrated that testing of multiple potential distributions of mutational effect can be important, because the results of the test can vary depending on which distribution was chosen for the analysis. For instance, the neutral hypothesis for an increase in rosette diameter could only be rejected when the best-fitting distributions were chosen. These results demonstrate the importance of the fit of the mutational distribution in our model. Selection was also detected on trichome density—for all distributions. These tests imply that natural selection led to the divergence in rosette diameter and trichome density of the Columbia and Landsberg populations of *A. thaliana*.

Trichome density is known to vary across natural populations, and has been associated with resistance to natural enemies (Mauricio 1998; Handley *et al.* 2005). Variation in the abundance of natural enemies alters the pattern of selection on trichome density (Mauricio and Rausher 1997). Molecular analyses of selection have suggested it is possible that positive selection on loci with major effects on trichome production has occurred during the diversification of *A. thaliana* populations in its recent history (Bloomer *et al.* 2012). Rosette diameter in Arabidopsis has been found to be under selection under some environmental conditions but not others. For example, rosette diameter was under positive selection for fall germinating plants, but not spring-germinating plants (Donohue 2002). When selection is observed on rosette diameter, there are frequently large selection gradients (Callahan and Pigliucci 2002). Although the genetic basis for the timing of bolting is not identical in Columbia and Landsberg *erecta* (Lee *et al.* 1994), they are both rapid-flowering lineages. Thus, there may not be extensive selective differentiation between them for this trait (Koornneef *et al.* 1991). There is evidence that the number of elongated axils is under consistent directional selection—potentially in both of these populations (Donohue 2002). Consequently, there may not be much selective differentiation for this trait. Rosette leaf number may be under similar consistent selection through indirect selection for larger plant size (Scheiner *et al.* 2000). Differences in selection history between Landsberg *erecta* and Columbia may also be due to differences in artificial selection during recent decades of cultivation. However, the small effective population size, selfing habit, and similar growth conditions (greenhouse and stock center) employed during propagation of the two accessions may have limited the effects of divergent artificial selection.

The original Landsberg lineage of *A. thaliana* underwent a program of mutagenesis shortly after its collection (Rédei 1992b), yielding the Landsberg *erecta* line. As a result, all phenotypic differences between Landsberg *erecta* and Columbia do not necessarily arise from the natural mutation process—a shortcoming that could compromise the validity of the test. However, as long as the distribution of mutation effects arising from mutagenesis is consistent with the distribution of effects from natural mutations, the mutagenesis of Landsberg should only weaken our power to detect directional selection on these traits and would not create a bias. Even if the distributions are different, if the functional form of their summed distribution matches one of the distributions specified in this paper, the prior mutagenesis should only weaken our power to detect directional selection on these traits and would not create a bias.

Application of this test to data from populations of *D. melanogaster* (Rice and Townsend 2012a) rejected neutrality for both traits analyzed. In contrast, application of the test in this paper yielded results that did not reject neutral evolution in three of the five traits studied. This discrepancy is perhaps consistent with expectation: the QTL estimated in the *Drosophila* study were known to be the product of strong, artificial directional selection. Although the power to detect the action of selection in the traits of the *Arabidopsis* populations was lowered by the comparatively small number of QTL identified (Ungerer *et al.* 2002; Symonds *et al.* 2005; Rice and Townsend 2012a), the lack of rejection of neutrality could indicate neutral evolution of these traits in these natural populations.

As in Rice and Townsend (2012a), our tests of selection assumed that a new mutation had a probability of fixation dictated by the diffusion-equation results of Moran (1960), Kimura (1962), and Gillespie (1974), which can accommodate a range of evolutionary dynamics. However, some evolutionary dynamics yield different probability of fixation equations, such as cyclically varying population size (Otto and Whitlock 1997), population subdivision with symmetric migration (Pollack 1966), and the presence of deleterious mutations at completely linked loci (Charlesworth 1994). A discussion of other equations for the probability of fixation can be found in Patwa and Wahl (2008).

Other methods to estimate the genetic variation at drift-mutation equilibrium could be explored in order to improve the estimate of the effective population size. Our method of calculating the genetic variation using heterozygosity of *A. thaliana* assumes a single stepwise mutational model for microsatellites (Ohta and Kimura 1973). However, a mutation can cause insertion or deletion of multiple repeat units (Di Rienzo *et al.* 1994). The estimator of genetic variation is upwardly biased when mutations allowing multiple repeat units occur (Xu and Fu 2004). Therefore, it can provide an overestimate of the effective population size. An overestimation of the population size can in turn lead to an underestimation of the probability of fixation for mutations drawn from the mutational effect distribution.

Other mutational effect distributions have been contemplated or analyzed in other contexts (*e.g.*, Gamma distributions; Keightley 1994; Shaw *et al.* 2002; Silander *et al.* 2007). Though they are less likely to be analytically tractable in this analysis, these distributions could also be tested against mutational effect data, and might improve the fit of the model to the data. Testing alternate distributions could yield more precise results, and help to relax the strength of distributional assumptions. For instance, inferred mutation rates and effect distributions can vary even with the genetic background of the founder used in the MA experiment. For example, Behringer and Hall (2016) inferred a higher mutation rate of *Schizosaccharomyces pombe* from MA data than Farlow *et al.* (2015), who evaluated MA data derived from a different founder. This difference in rates could be attributed to differences in the founder’s genetic background.

Lastly, with regard to the model framework presented here, we could conduct explicit accounting of time instead of implicitly incorporating time (Rice and Townsend 2012a) to improve our model. Since it takes time for mutations to generate and fix, the time a trait has had to evolve affects the amount of substitutions that cause each QTL effect in our model. Accounting for the amount of time a trait has had to evolve could affect the likelihood of the number of substitutions (given the level of selection) that are responsible for each QTL effect. An explicit incorporation of time would help determine the number of substitutions at each locus, which would provide additional information that could lead to a more precise measure of selection strength.

Our expansion of the repertoire of distributions to which mutational effects can be fit has facilitated more accurate detection of phenotypic selection from quantitative genetic data on divergence between lineages and observed effects of mutation accumulation on phenotype. Increases in the paired availability of QTL and MA data of phenotypic traits for diverse organisms of interest for the selective history on phenotype would enable additional application of this quantitative test for directional selection, and therefore continued progress in understanding the history of directional selection leading to the abundant phenotypic divergence observed between diverse lineages and populations.

## Acknowledgments

The authors thank two anonymous reviewers for their helpful comments on the manuscript.

## Footnotes

*Communicating editor: C. D. Jones*Supplemental material is available online at www.genetics.org/lookup/suppl/doi:10.1534/genetics.116.199190/-/DC1.

- Received January 31, 2017.
- Accepted May 22, 2017.

- Copyright © 2017 by the Genetics Society of America