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.
A depiction of our approach of inferring the probability distributions of mutational effects from a set of mutation accumulation data resulting from unbiased mutations, and another set resulting from downward-biased mutations. (A) A constructed example of changes in phenotype resulting from unbiased mutation in seven mutation accumulation lines. (B) A constructed example of changes in phenotype resulting from mutations with downward bias in seven mutation accumulation lines. (C, D) The Gaussian distribution characterizing mutational effects inferred from data depicted in (A). (E, F) The LSDZ distribution characterizing mutational effects inferred from data depicted in (A). (G, H) The LAPZ distribution characterizing mutational effects inferred from data depicted in (A). (I, J) The LAPZ distribution characterizing mutational effects inferred from (B). (K, M, O, Q) A depiction of a probability of fixation of mutations that does not change with increasing phenotypic values. (L, N, P, R) A depiction of a probability of fixation of mutations that increases with increasing phenotypic values. (S) The product of the distribution depicted in (C), and the function depicted in (K): in this case, a symmetrical distribution of fixed mutations with a positive mode. (T) The product of (D) and (L): in this case, an asymmetrical distribution of fixed mutations with a positive mode. (U) The product of (E) and (M): in this case, a symmetrical distribution of fixed mutations with a near zero mode. (V) The product of (F) and (N): in this case, an asymmetrical distribution of fixed mutations with a positive mode. (W) The product of (G) and (O): in this case, a symmetrical distribution of fixed mutations with a negative mode. (X) The product of (H) and (P): in this case, an asymmetrical distribution of fixed mutations with a positive mode. (Y) The product of (I) and (Q): in this case, an asymmetrical distribution of fixed mutations with a zero mode. (Z) The product of (J) and (R): in this case, an asymmetrical distribution of fixed mutations with a positive mode.
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.
A depiction of our approach of inferring the likelihood values of a range of strength of selection values from a distribution in Figure 1 and quantitative trait locus data for a phenotype. (A–D) Examples of a single sample of mutations, which are shown by the outlines, drawn from the distribution in (S), (T), (W), and (X) of Figure 1, respectively, and the experimentally observed positive and negative QTL, shown by the colored bars. (E) A plot of the likelihood values across a range of strength of selection values. Samples drawn from the distribution in (S) seem to fit better than samples drawn from distribution (T). The region of the plot above the dashed red line correspond to the 95% confidence interval. The point at the y-axis (black) corresponds to neutral evolution and lies inside the 95% confidence interval, while the point (red) in the first quadrant corresponds to evolution driven by natural selection, and lies outside the 95% confidence interval. Therefore, we cannot reject the neutral hypothesis. (F) A plot of the likelihood values across a range of strength of selection values. Samples drawn from the distribution in (X) seem to fit better than samples drawn from distribution (W). The region of the plot above the dashed red line correspond to the 95% confidence interval. The point at the y-axis (red) corresponds to neutral evolution and lies outside the 95% confidence interval, while the point in the first quadrant (black) corresponds to evolution driven by natural selection and lies inside the 95% confidence interval. Therefore, we can reject the neutral hypothesis in favor of the hypothesis of evolution driven by natural selection and estimate that the strength of selection corresponds to the strength of selection indicated by this point.
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
(9)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 × 104 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 xi and gi 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.
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 = (J1, J2, …, Jn), where Ji is the number of substitutions at the ith 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 Ji substitutions under a strength of selection c. The probability density function of the sum is calculated by Ji 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 xj. 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 Ji 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.
The effect sizes and the SEs of rosette diameter, number of elongated axils, number of rosette leaves at bolting, bolting time, and trichome density QTL. The sign of the effect size corresponds to the direction of effect of the Columbia allele on the phenotype.
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).
Frequency distribution of changes in quantitative traits during mutation accumulation. (A) Distribution of the changes in the number of elongated axils. (B) Distribution of the changes in trichome density. (C) Distribution of the changes in number of rosette leaves at bolting. (D) Distribution of the changes in rosette diameter. (E) Distribution of the changes in bolting time.
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 xi and gi 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.
The likelihood for a range of selection coefficients was calculated for the QTL. We used an effective population size N = 2.3 × 106 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 × 105 Markov chain Monte Carlo samples from the distribution of mutation effects for every value of c and Hi. 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.
The normalized AIC of the Gaussian and Laplace distributions for the number of rosette leaves at bolting, bolting time, rosette diameter, number of elongated axils, and trichome density of A. thaliana. The normalized AIC for the LSPZ distribution was excluded, because its measure was very high—beyond the scale of the other distributions.
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).
Maximum likelihood estimates and the corresponding 95% confidence intervals of the selection values in units of 10−5 for bolting time, number of rosette leaves at bolting, number of elongated axils, rosette diameter, and trichome density using all nine distributions. There were no upper bounds for the 95% confidence intervals for trichome density, because the log likelihood values described a positive (but diminishing) slope as the strength of selection became larger. This phenomenon is a consequence of all of the observed QTL for trichome density exhibiting a single direction. Although not all estimates of the strength of selection statistically reject a value of zero, all traits are estimated to be under positive selection toward high values of the five quantitative traits in Columbia, consistent with the observation of higher values of these traits in Columbia, and the apparent downward pressure of mutational effects on the traits.
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
