Estimating Selection on Nonsynonymous Mutations

Abstract

The distribution of mutational effects on fitness is of fundamental importance for many aspects of evolution. We develop two methods for characterizing the fitness effects of deleterious, nonsynonymous mutations, using polymorphism data from two related species. These methods also provide estimates of the proportion of amino acid substitutions that are selectively favorable, when combined with data on between-species sequence divergence. The methods are applicable to species with different effective population sizes, but that share the same distribution of mutational effects. The first, simpler, method assumes that diversity for all nonneutral mutations is given by the value under mutation-selection balance, while the second method allows for stronger effects of genetic drift and yields estimates of the parameters of the probability distribution of mutational effects. We apply these methods to data on populations of Drosophila miranda and D. pseudoobscura and find evidence for the presence of deleterious nonsynonymous mutations, mostly with small heterozygous selection coefficients (a mean of the order of 10⁻⁵ for segregating variants). A leptokurtic gamma distribution of mutational effects with a shape parameter between 0.1 and 1 can explain observed diversities, in the absence of a separate class of completely neutral nonsynonymous mutations. We also describe a simple approximate method for estimating the harmonic mean selection coefficient from diversity data on a single species.

SURVEYS of DNA sequence polymorphisms in many species have revealed substantial variation in the aminoacid sequences of proteins, although the nonsynonymous nucleotide site diversity is usually much less than that for silent variants (Li 1997). This is consistent with the action of purifying selection on protein sequences, removing deleterious amino acid mutations from the population while neutral or nearly neutral silent variants can persist (Kimura 1983; Li 1997). It is clearly important to characterize the distribution of fitness effects of new nonsynonymous mutations. This distribution is relevant to a broad range of problems in evolutionary genetics, and a variety of methods have been used to characterize it, including direct estimates from mutation-accumulation experiments and indirect estimates from comparisons of DNA sequences among related species (Keightley and Eyre-Walker 1999). The extent, nature, and magnitude of selection on amino acid variants are also relevant to understanding the relation between human disease and genetic variation (Sunyaev  et al. 2001; Wright  et al. 2003).

Several methods have been used to detect purifying selection from data on variability in natural populations and to estimate the parameters describing such selection. An important method was introduced by Sawyer and Hartl (1992), based on the McDonald-Kreitman test (McDonald and Kreitman 1991). It compares the ratio of the number of within-species nonsynonymous polymorphisms to the number of synonymous polymorphisms and the corresponding ratio for fixed differences between a pair of related species. This ratio of ratios is called the “neutrality index” (Rand and Kann 1996). It is expected to be greater than one if there is predominantly purifying selection against nonsynonymous variants, since selection is less effective in reducing the level of polymorphism for deleterious variants than in preventing their fixation (Kimura 1983). This approach can be incorporated into maximum-likelihood or Bayesian methods for estimating N_es, where N_e is the effective population size and s is the selection coefficient against nonsynonymous mutations. The selective effects of mutations are usually assumed to be codominant (in the case of nuclear genes), and constant across sites within a gene (e.g., Nachman 1998; Bustamante  et al. 2002). Recent extensions allow for a distribution of selection coefficients at different sites (Piganeau and Eyre-Walker 2003) or varying degrees of dominance (Williamson  et al. 2004).

Overall, this method has been more successful in detecting and estimating purifying selection on nonsynonymous variants in mitochondrial genomes than in the nuclear genome (Nachman 1998; Rand and Kann 1998; Weinreich and Rand 2000), with a few exceptions such as Arabidopsis thaliana (Bustamante  et al. 2002), Drosophila miranda (Bartolomé  et al. 2005), and humans (Bustamante  et al. 2005). The fixed selection coefficient model was fitted by Nachman (1998) to 17 animal mitochondrial DNA data sets with neutrality index values >1 and gave N_es estimates predominantly between 1 and 3; Bustamante  et al. (2002) estimated N_es on a gene-by-gene basis and obtained a mean value of ∼1 for 12 nuclear genes in A. thaliana (in both cases, N_e is the haploid effective population size, as is appropriate for the ploidy and breeding systems in these cases). A modification of this approach, allowing a normal distribution of N_es values across different genes, was applied to data on D. melanogaster, indicating a mean N_es of ∼3.5 (Sawyer  et al. 2003). The fits of a gamma distribution of selection coefficients across individual nonsynonymous sites to animal mitochondrial data sets (Piganeau and Eyre-Walker 2003) gave much larger but noisier estimates of mean N_es.

An alternative approach is to fit the observed frequency spectra of nonsynonymous and silent/synonymous variant frequencies to the distributions expected under mutation-selection-drift equilibrium (Sawyer 1994; Sawyer  et al. 1987). Variants of this approach have been developed by Hartl  et al. (1994), Akashi (1999), Fay  et al. (2001), Bustamante  et al. (2003), and Williamson  et al. (2004, 2005). Either a fixed given N_es value or a distribution of N_es values across loci is assumed.

These studies have yielded large differences among estimates of the scaled selection parameter N_es (or its arithmetic mean) for nonsynonymous sites under purifying selection, ranging from values of the order of 1 to several hundred, depending on the methods and species used. Varying conclusions about the proportion of sites subject to positive, as opposed to purifying selection, have also been reached; for example, compare Sawyer  et al. (2003) with Bierne and Eyre-Walker (2004), who estimated that 94 and 25% of nonsynonymous mutations distinguishing D. simulans and D. melanogaster were fixed by positive selection, respectively.

Methods that fit details of frequency spectra to equilibrium models are clearly highly vulnerable to departures from equilibrium, and there is increasing evidence that many of the model systems used for the study of natural variation, as well as human populations, are subject to such effects (Andolfatto and Przeworski 2000; Williamson  et al. 2005). Incorporating even the simplest model of demographic change into selection models is computationally extremely demanding (Williamson  et al. 2005). Methods that ignore the details of the distribution of variant frequencies may thus be preferable to potentially more powerful methods that exploit all the features of the data. Another problem is that many of the methods outlined above assume that amino acid mutations in a given gene are unidirectionally from wild type to selectively deleterious or vice versa. However, unless the magnitude of N_es for all mutations is much greater than one, there will be a flux of amino acid substitutions over evolutionary time, such that some fraction of sites will be fixed for selectively deleterious alleles and can therefore back mutate to create fitter variants, as in models of codon usage bias (Li 1987; Bulmer 1991; McVean and Charlesworth 1999). Only the model of Piganeau and Eyre-Walker (2003) has explicitly incorporated this possibility.

In this article, we develop a model of nonsynonymous site variation and evolution that includes reversible mutation, as in the standard models of codon usage evolution. We use this to estimate the strength of selection on amino acid variants, by exploiting the difference in the responses of nonsynonymous and synonymous/silent variants to differences in effective population sizes between related species. The basic idea is that variants subject to sufficiently strong purifying selection will not increase much in abundance as effective population size increases, whereas neutral or nearly neutral diversity is expected to increase in proportion to N_e. The extent to which nonsynonymous diversity differs between species with different synonymous diversity values should thus shed light on the prevalence and strength of purifying selection. We also show how to provide useful bounds on selection parameters when data on only one species are available.

MATERIALS AND METHODS

Source of data:

We used published DNA sequence information on X-linked and autosomal genes for D. miranda (Yi  et al. 2003; Bartolomé  et al. 2005), removing genes for which there was either evidence for departure from neutrality (Annx, swallow) from the HKA test or no sequence data from D. affinis, the outgroup species used to estimate interspecific divergence from D. miranda and D. pseudoobscura (Bartolomé  et al. 2005). Polymorphism data on D. pseudoobscura were compiled from published population surveys; the gene exu2 was not used, since it showed evidence for selection on the basis of haplotype structure (Yi  et al. 2003) and the HKA test (this study; data not shown). These data were obtained from GenBank accessions, and sequences were aligned using the program SeAl (http://evolve.zoo.ox.ac.uk/). Details of the genes used and the relevant references are provided in supplemental material at http://www.genetics.org/supplemental/.

In total, 17 D. miranda genes (13,309 nonsynonymous sites and 9077 silent sites, 51% in introns) and 14 D. pseudoobscura genes (10,828 nonsynonymous sites and 10,245 silent sites, 65% in introns) were used. Sample sizes were 11 or 12 alleles per gene for D. miranda and 7–139 per gene for D. pseudoobscura. No adjustments for different effective population sizes for X-linked vs. autosomal genes were made, as mean diversities are similar for these two categories in both species, consistent with the action of sexual selection on males (see supplementary material at http://www.genetics.org/supplemental/ and Yi  et al. 2003; Bartolomé  et al. 2005). Polymorphism and divergence estimates were calculated using DnaSP (Rozas  et al. 2003). The estimates of divergence from D. affinis for the D. miranda loci were obtained from Bartolomé  et al.'s (2005) Table 3. Unfortunately, only three loci are in common between the two data sets, so that we have to treat the two sets of genes as representing more or less independent random samples from the genomes of the two species.

Computational methods:

The case of “arbitrary purifying selection” described in theoretical  framework (Equations 5) requires integration of the formulas for the nonsynonymous nucleotide site diversities π_A and rates of substitution K_A, over an assumed distribution of selection coefficients, ϕ(s). The relevant expressions involve integrals representing the sojourn times of codominant autosomal mutations, given the heterozygous selection coefficient s, the breeding adult population size N, the effective population size N_e, and the mutation rate per nucleotide site per generation u. These were obtained from the known solutions to the relevant diffusion equations (Kimura and Ohta 1969b; McVean and Charlesworth 1999). To predict π_S, we used similar equations as for π_A, but setting s = 0. All computations were implemented in the statistical programming language “R” (version 1.9) (Ihaka and Gentleman 1996; R-Project 2005).

Most mutations generated by a given ϕ(s) distribution have effects that can be handled by the diffusion methods employed here. However, for very strongly or weakly selected mutations, the formulas require more numerical accuracy than the 15 digits available in double-precision floating variables. The standard approximations for fixation probabilities and sojourn times of mutations, for the respective cases of either neutrality or strong selection, were then used (Haldane 1924, 1927; Kimura 1962; Kimura and Ohta 1969a,b).

To integrate over the distribution of selection coefficients, we partitioned the range of mutational effects of interest, from effectively neutral (s = 10⁻¹⁰) to lethal (s ≥ 1), into I groups small enough to assume constant mutational effects within each group. We found that π_A and K_A values computed from I = 30, 100, and 300 equidistant steps on a log scale were accurate to ∼2, 0.2, and 0.02% relative error, respectively. For each bin, we then independently computed (i) the probability P_i that a mutation will have an effect that belongs to bin i, obtained by integrating ϕ(s) from the lower to the upper limit of the bin, and (ii) all interesting quantities of the model, given the average mutational effect characterizing that bin. To get the overall result for a parameter of interest, we summed the corresponding results for the parameter over all bins, weighting the value for each bin by P_i.

Experimenting with different ϕ(s) functions, two special cases became obvious. Some mutations have effects smaller than the lower integration limit (s = 10⁻¹⁰). We added the probability mass of these mutations to the first bin. Since these are effectively neutral mutations, this amounts to full integration of the ϕ(s) down to neutrality. The distribution of s was truncated at s ≥ 1, keeping a record of the fraction of mutations that fall into this category; this represents dominant lethal mutations caused by amino acid substitutions, which are probably extremely rare.

For each quantity involved in the model, a plot over the whole range of values of s generated from a given ϕ(s) was produced, and the smoothness of transitions to neutrality and to strong selection equilibrium was used to verify the accuracy of the computations. At statistical equilibrium under drift, mutation, and selection at each site, we expect equal rates of substitution between preferred and unpreferred amino acids at sites with a given s, which was confirmed in our plots.

The d function was chosen to make small differences look large, a property needed for the simplex optimization routine (Amoeba) that we used (Nelder and Mead 1965), implemented in R. All ratios were rounded to at least six digits to ensure that bad fits were detectable. An optimization attempt was considered as successful if the data could be predicted with six-digit accuracy, using I = 100 integration steps; this excluded cases where optimization stopped without getting close to the data, wrongly suggesting that the data had been fitted. In all cases, the estimates resulting from the first Amoeba run were used as starting values for a second run, to make sure that the first result was not simply a local optimum. The parameter values that satisfied the optimization criterion after the second run were used to compute the results reported here, with an increased accuracy (I = 300).

To simplify the calculations using Equations 5 below, we assumed that the effective population size N_e is equal to the size of the breeding population, N. To obtain the relevant numerical values, we used Equation 1a below to estimate N_e from the observed mean silent diversity, assuming a mutation rate u. We mostly used u = 1.5 × 10⁻⁹ in the calculations, a value widely used for Drosophila (Powell 1997). Since this is not firmly established, analyses with other mutation rates were also done, to check sensitivity to u.

Once the parameters determining variability from Equations 5a were estimated, they were used in Equation 5b to predict the rate of amino acid substitutions arising from sites under purifying selection, assuming an arbitrary value for the unknown ancestral N_e. This in turn may be used to estimate the fraction of selectively advantageous amino acid substitutions, by comparing the prediction from Equation 5b with the observed divergence between species, similarly to Equation 4b.

Statistical analyses:

Preliminary analyses of the data showed that D. pseudoobscura had much greater silent diversity than D. miranda. Its genes also consistently show an excess of rare variants over neutral expectation (Machado  et al. 2002; Schaeffer 2002), in contrast to what is seen in D. miranda (Yi  et al. 2003; Bartolomé  et al. 2005). This suggests that D. pseudoobscura has undergone a recent period of population expansion and is therefore not likely to have reached its final level of neutral or nearly neutral diversity. The measure of neutral variability provided by Watterson's θ_w estimator, based on the number of segregating sites for silent diversities (Watterson 1975), is probably closer to the equilibrium value than the pairwise nucleotide site diversity estimator π, since new variants arising after an increase in N are predominantly rare (Tajima 1989). We therefore expect θ_w to provide a better estimate of the equilibrium neutral/nearly neutral diversity for D. pseudoobscura than π and have accordingly used θ_w for silent sites for both species. For sites under selection, there is no cogent reason to use θ_w, and so we used π as the diversity measure for nonsynonymous variants. We obtained very similar results if θ_w is used for both types of sites. For simplicity, we use the symbol π to denote diversity estimates for both cases in what follows.

Mean values across genes of the diversities and divergence statistics were used to estimate the parameters of the models. Weighted mean diversity estimates were obtained using the inverse of the sum of the estimated sampling and stochastic evolutionary variances of diversity as the weight for a given gene, obtained from the standard formulas under neutrality and free recombination (different formulas apply to the Watterson and pairwise diversity estimators as described in Chapter 10 of Nei 1987). This procedure is heuristic, given that there is some linkage disequilibrium among sites within genes and that nonsynonymous mutations are known to be subject to selection, but it provides a simple approximate way of accounting for different levels of noise across genes. Divergence values, K, were weighted by the number of sites involved. Their means across genes were then used to estimate K_A/K_S. For comparison, unweighted means of diversities and divergence were also computed and used for parameter estimation.

To assess the variability of our estimates for both methods of weighting, we generated 1000 “observations” by bootstrapping the diversity and divergence data across loci, as described by Bartolomé  et al. (2005). The upper and lower fifth percentiles of the distribution of each parameter were used as approximate 90% confidence intervals for the parameters in question; this provides a convenient basis for assessing 5% P-values for the null hypothesis that this parameter has a value of zero, in a one-tailed test.

THEORETICAL FRAMEWORK

General framework:

We assume that the sequences of population samples of many different genes are known, giving reliable estimates of diversities for nonsynonymous and silent or synonymous sites in each of two species and corresponding divergences from a third (Figure 1). Differences in N_e exist between the two species for which diversity data are available; the effective population size of species i (i = 1 or 2) is denoted by

⁠, where the smaller i denotes the species with the smaller N_e. We assume an infinites-sites model with autosomal inheritance (Kimura 1971). Silent mutations are assumed to be selectively neutral or sufficiently close to neutrality that their evolutionary fates are well described by the neutral model. Each site is assumed to evolve independently; i.e., recombination is sufficiently frequent that Hill-Robertson interference effects can be neglected. Mutation rates and selection coefficients at each site are assumed to be independent, and random mating is assumed for both species.

To keep our model simple, we assume that there are at most two types of amino acid at a given nucleotide site: “preferred” and “unpreferred,” with at most one of the four possible nucleotides corresponding to the preferred state (our methods do not require us to identify preferred and unpreferred states from the data). At some sites, all possible variants may be effectively neutral, meaning in practice that variants are subject to a similar level of selection to silent mutations (Figure 2D). The fraction of such “neutral” nonsynonymous mutations is denoted by c_n. We treat other sites as in Figure 2A, where nonsynonymous nucleotide site mutations can result in a change to a selectively deleterious amino acid (if the nucleotide in question codes for the preferred state), to a favorable amino acid (if the site is fixed for an unpreferred amino acid, and the mutation leads back to the preferred state), or else to a selectively neutral variant (if the site is fixed for an unpreferred amino acid, and the mutation causes a change to a different but selectively equivalent deleterious variant). Because some mutations at such sites are deleterious, we refer to these sites as undergoing purifying selection.

The selection coefficient s for selectively deleterious variants refers to the reduction in fitness of heterozygotes; this may vary according to the site and amino acid in question. For convenience of calculation, dominance is assumed to be intermediate. As in standard models of codon usage bias (Li 1987; Bulmer 1991; McVean and Charlesworth 1999), if selection at these sites is sufficiently weak, unpreferred mutations can become fixed and then mutate back to the preferred state, contributing to a flux of amino acid substitutions as well as to polymorphism. This model is not completely general (see Figure 2, B and C), but should suffice as a guide to the basic processes involved, especially when selection is fairly strong in relation to drift.

Occasionally, a new adaptive mutation (distinct from a back mutation at a site fixed for a deleterious mutation) may arise at a nonsynonymous site, perhaps in response to a change in the environment. This is assumed to spread rapidly to fixation if it survives initial stochastic loss. The substitution rate per generation per site for mutations that fall into this category is denoted by c_au, where u is the expected mutation rate per nucleotide site, and c_a measures the substitution rate as a fraction of all mutations. c_a is the product of the frequency of adaptive mutations and their fixation probability, integrated over all advantageous effects. We assume that c_a ≪ 1, since favorable mutations are likely to be rare.

Strong purifying selection:

Arbitrary purifying selection:

The validity of the assumption that N_es > 1 for all nonneutral nonsynonymous mutations is, however, questionable. If this assumption is relaxed, the formulas for the equilibrium diversity at nonsynonymous sites become more complex, and the possibility of a contribution to K_A from sites subject to purifying selection must also be considered, using a probability distribution ϕ(s) of selection coefficients for mutations subject to selection. We now consider this problem in detail.

As far as diversity is concerned, we note first that a nonsynonymous site that has become fixed for a mutant nucleotide coding for a deleterious amino acid can either mutate back to the original nucleotide or mutate to another nucleotide coding for a deleterious amino acid. If mutation rates among all four nucleotides were equal, the probability of back mutation to the original state would be

⁠; in general, however, inequalities in mutation rates are likely to make this probability different from

\(\frac{1}{3}\)

⁠, and so we represent it as 1/κ.

RESULTS

We now present the results of applying these methods to the diversity data on D. miranda and D. pseudoobscura described in materials  and  methods (D. miranda is designated as species 1 and D. pseudoobscura as species 2 in what follows). In view of such problems as the lack of overlap between the genes used in the two species, and the disparity in sample sizes between studies, the results based on these data should be regarded as merely provisional and illustrative of the methods.

Strong purifying selection:

We first present the results of applying the expectations for the case of strong purifying selection. Divergence values from D. affinis (species 3) were estimated for the 17 genes surveyed in D. miranda; the divergence values for these genes between D. affinis and D. pseudoobscura were almost identical (Bartolomé  et al. 2005), so that only the D. miranda results were used here. Both weighted and unweighted estimates of mean diversities were obtained, as described in materials  and  methods; statistical uncertainty was assessed by bootstrapping across genes.

The results are displayed in Table 1. The main conclusion is that an estimate of

substantially greater than one is supported by the bootstrapping results, even using the unweighted estimators, which yield lower values than the weighted estimators. The distribution of

\(N_{\mathrm{e}_{1}}s_{\mathrm{h}}\)

is, however, very wide, and infinite values were sometimes generated for the weighted data. The estimate of N_es_h for D. pseudoobscura was, of course, proportionately larger, corresponding to the large N_e-value estimated from silent site diversities (these suggested a 5.8-fold higher N_e for D. pseudoobscura).

The estimated proportion of neutral sites was smaller with the unweighted estimates than with the weighted ones, but both suggested a value of a few percent. There was a very wide distribution of values of the proportion of fixations due to adaptive mutations in both cases, and a zero value could not be ruled out by the data, although the value estimated from the data exceeded 50% for the unweighted estimate.

Arbitrary purifying selection:

We now relax the assumption that π_A is close to the deterministic prediction by applying the model behind Equations 5 to the same data. The assumptions of a gamma distribution of mutational effects, and that 0, 2.5, or 5% of all nonsynonymous mutations are neutral, yielded the parameters reported in Tables 2 and 3 for the weighted and unweighted data, respectively. With values of c_n ≥ 7.5% we rarely, if ever, found parameters for a gamma distribution that fitted the data, and the fraction of bad fits for c_n = 5% was 48.5% for 1000 bootstraps from the weighted data, whereas the mean of the unweighted data could not be fitted assuming c_n ≥ 5%. With the gamma distribution, a significant fraction of mutations had such small selection coefficients that N_es < 1 or even 0.5. The results of McVean and Charlesworth (1999) suggest that, with N_es < 0.5, both diversities and substitution rates are nearly equivalent to those for neutral sites; in addition, the intensity of selection on synonymous mutations that change codon usage from preferred to unpreferred codons in D. miranda seems to be close to N_es = 0.5 (Bartolomé  et al. 2005). It thus seems reasonable to use this value as the boundary for designating mutations as effectively neutral; the sum of c_n and the fraction of effectively neutral mutations generated by the fitted gamma distribution is denoted by c_ne in Tables 2 and 3. For completeness, we also show results using N_es = 1 as the boundary.

Despite the time-consuming nature of the multiple integrations involved in implementing this model, we attempted 1000 bootstrap replications to assess the reliability of our parameter estimates. Out of these bootstrap computations for the weighted or unweighted data in Table 2 or 3, only 80.3 and 90.1%, respectively, could be fitted assuming c_n = 0%; 75.6 and 65.6% could be fitted assuming c_n = 2.5%. Technically, our estimates of the distributions of mutational effects involve only the shape and the location parameter for the gamma distribution of s, together with values for

⁠. However, in Tables 2 and 3 we also report more intuitively meaningful measures of the underlying distribution of all selection coefficients at sites capable of mutating to nonlethal mutations that fall above the threshold of effective neutrality, scaled by N_e. These include the arithmetic and harmonic means, the coefficient of variation, and the lower and upper fifth percentiles of the distribution of s-values for these mutations.

As reviewed in more detail in the discussion, the bulk of the nonneutral mutations segregating in the population come from the more weakly selected tail of the distribution (Sunyaev  et al. 2001). The arithmetic mean of the distribution of selection coefficients among these polymorphic variants is much closer to the harmonic than to the arithmetic mean of the distribution for new mutations, and so the harmonic means in Tables 2 and 3 are more relevant than the arithmetic means to the properties of mutations found in populations. It is notable from Tables 2 and 3 that the means for D. pseudoobscura were smaller than might be expected from the ratio of effective population sizes and the corresponding means for D. miranda; this reflects the lower fraction of mutations that fall into the effectively neutral class when the effective population size is larger, reducing the average selection coefficients for mutations in the other class.

Given an ancestral N_e-value, we can also predict the substitution rate for mutations under purifying selection. This can be compared with the observed rate to estimate the proportion of adaptive substitutions, using Equation 4. We report three alternative values:

⁠,

\(P_{\mathrm{a}_{2}}\)

⁠, and

\(P_{\mathrm{a}_{3}}\)

⁠, which assume ancestral N_e-values equal to the estimated current N_e for D. miranda, the current N_e of D. pseudoobscura, and the mean of these, respectively. For each value, we give the approximate lower and upper fifth percentiles obtained by bootstrapping.

The main result of Tables 2 and 3 is solid support for the conclusion that ∼90% or more of all amino acid mutations have significantly deleterious effects. It is also remarkable that the estimates of the biologically important harmonic mean selection coefficient are close to those using the strong selection assumption, taking into account the statistical noise and the uncertainty regarding c_n for the arbitrary purifying selection model. Similar conclusions can be drawn for comparisons between c_ne in the arbitrary purifying selection model and c_n in the strong selection model. Comparing the weighted and unweighted estimates shows that the weighting procedure influences the estimates, but the noise in the data is larger than the noise from uncertainty about the best weighting procedure. While our definitions of effective neutrality (N_es < 0.5 or < 1) have little influence on the arithmetic mean and the upper fifth percentile of the distribution of mutational effects, their influence on the harmonic mean and the lower fifth percentile is greater. In no case, however, did the ∼1–3% of mutations that fall in the range 0.5 < N_es < 1 change our conclusions regarding the prevalence of deleterious mutations.

To visualize the gamma distributions of mutational effects estimated from the weighted data, we plotted our best-fitting estimates with assumed values of c_n = 0, 2.5, and 5% (Figure 3). This involved estimating effective population size from Equation 1a, using the “standard” mutation rate of 1.5 × 10⁻⁹. The distributions for the two smaller values of c_n peaked at selection coefficients ∼10⁻⁴, fairly close to the arithmetic means obtained from Table 2 for mutations that are not effectively neutral, whereas the harmonic mean was about one-tenth of this.

DISCUSSION

Robustness of the strong purifying selection model:

We found quite good agreement between the results of the strong purifying selection model and the model with a gamma distribution of selection coefficients for the estimates of the harmonic mean of the selection coefficient. This suggests that estimates of the magnitude of this important parameter are robust to the assumptions used, providing that a sizeable fraction of nonsynonymous polymorphisms is nonneutral. Some approximations that relax the assumptions of the strong purifying selection model, but that do not depend on the details of the distribution, are explored in the appendix. These provide very general methods for estimating N_es_h and again suggest that the magnitudes of the estimates of the main parameters of interest are fairly robust to details of the assumptions.

This result is remarkably robust, since we do not need to know N_e, u, s_h, or dominance coefficients. However, caution is necessary if Equation 8b gives values near 1, since this suggests that drift is probably too strong to be neglected. In this case, the true strength of selection for the deleterious mutations will be larger than predicted by this approach. The values estimated from this method are N_es_h ≈ 2.9 and 4.9 for D. miranda and D. pseudoobscura, as estimated for the respective weighted data. These values are little more than a third of the respective values estimated from the most precise methods and lie outside the confidence intervals for the latter. This very simple formula is therefore too crude for precise estimates, but seems to work reasonably well as a rough first estimate for the lower bound of N_es_h. It can be applied only when there is evidence that a substantial proportion of segregating nonsynonymous variants experience purifying selection, as in the present case (Bartolomé  et al. 2005).

One caveat should be noted concerning the estimates of N_e for D. pseudoobscura that we have been using. This assumes that silent sites are effectively neutral, which we have taken to mean that N_es is of the order of ≤0.5 McVean and Charlesworth (1999). While D. miranda synonymous sites seem to satisfy this condition (Bartolomé  et al. 2005), this condition clearly cannot apply to D. pseudoobscura, given that mean silent diversity is four to five times higher than that in D. miranda, unless selection on synonymous sites is much weaker in D. pseudoobscura. Akashi and Schaeffer (1997) estimated N_es against unpreferred codons to be 4.6 (95% confidence interval 2.4–12.1) for Adh plus Adhr in D. pseudoobscura, (although selection for preferred codons was negligible); this result may be in part confounded by the effects of population expansion. If N_es for synonymous sites in D. pseudoobscura is indeed higher than that in D. miranda, then mean silent site diversity (which includes a large contribution from synonymous sites) will yield an underestimate of 4N_eu for this species. Accordingly, N_es for nonsynonymous sites will be underestimated and the proportion of effectively neutral nonsynonymous mutations overestimated, since the ratio of nonsynonymous diversity relative to effectively neutral diversity will be overestimated for this species.

Sensitivity to mutation rates, mutational bias, and recombination rates:

Estimates derived from the arbitrary purifying selection model are surprisingly insensitive to plausible changes in mutational bias and mutation rate. Drake  et al. (1998, p. 1673) estimated from laboratory experiments that 8.5 × 10⁻⁹ mutations per base per generation happen in D. melanogaster. Powell (1997, p.369–371) reported rates between 0.67 × 10⁻⁹ and 3.3 × 10⁻⁹, estimated from divergence between the D. melanogaster and D. obscura groups, assuming that divergence happened 30 MYR ago and that each year represents ∼10 generations. McVean and Vieira (2001) estimated a rate of 1.5 × 10⁻⁹ (95% C.I. = 1.0 × 10⁻⁹–2.5 × 10⁻⁹), assuming 2–4 MYR divergence between D. melanogaster and D. simulans. Others (Andolfatto and Przeworski 2000; Przeworski  et al. 2001) found rates of between 0.6 × 10⁻⁹ and 4.75 × 10⁻⁹. Thus 0.5 × 10⁻⁹ and 8 × 10⁻⁹ appear to be reasonable choices for the most extreme lower and upper credible limits for mutation rates in Drosophila. As can be seen in Table 4, when c_n = 2.5%, a mutational bias of κ = 1 leads to the largest differences from our other estimates, while large changes in mutation rate seem to have only minor effects. Other assumed values of c_n yield the same conclusion, although the resulting estimates differ because of the strong influence of c_n (see Tables 2 and 3).

Three genes in our D. pseudoobscura data set (Amy1, eve, and exu1) are located on Muller's C, a genomic region that is segregating for paracentric inversions (Dobzhansky and Powell 1975) and therefore has a highly reduced recombination rate. These three genes violate the assumption of independence of sites much more than the other genes. We ran 500 bootstraps for a variance weighted data set without these genes under the assumption of c_n = 0%. Results indicate that the inclusion of these genes does not strongly affect our parameter estimates, as the confidence intervals of our estimates for the reduced data set mostly overlap our estimates for the full data set (data not shown).

The reliability of estimates of the proportion of adaptive mutations:

As can be seen from Figure 4, there is a large influence of the value of the (unknown) ancestral N_e on the estimate of the proportion of adaptive mutations. As expected from the fact that more slightly deleterious mutations can be fixed in smaller populations, fewer adaptive mutations are inferred with smaller ancestral N_e-values (Eyre-Walker 2002). Again, the mutation rate and mutational bias do not greatly affect the estimates. The combination of Figure 4 with the wide confidence intervals for P_a from Tables 1–3 raises the question of whether this approach can determine the presumably small fraction of adaptive mutations with any precision. Larger data sets and the use of the same sets of genes in the two species being compared may help to narrow the error bounds on the estimates.

However, potential difficulties still remain. One is the assumption of free recombination among variants. Linkage increases the rate of fixation of deleterious mutations while decreasing that for advantageous mutations (Birky and Walsh 1988); close linkage can have important effects even with weak selection (McVean and Charlesworth 2000; Kim 2004). But this assumption seems reasonable for D. pseudoobscura and its relatives, with their high rates of recombination and lack of linkage disequilibrium within genes (Dobzhansky and Powell 1975; Schaeffer and Miller 1993; Yi  et al. 2003). In addition, departures from equilibrium due to demographic effects may introduce errors into the estimates. As discussed in materials  and  methods, our choice of θ_W instead of π as a synonymous diversity estimator is intended to minimize the effects of the population expansion that seems to have occurred in D. pseudoobscura (Machado  et al. 2002; Schaeffer 2002). However, we cannot exclude population bottlenecks in the distant past that could have led to a higher contribution of deleterious mutations to K_A/K_S relative to π_A/π_S and thus have elevated the estimate of the proportion of adaptive substitutions (Eyre-Walker 2002), a problem common to all methods used to date.

The reliability of estimates of parameters of the distribution of mutational effects:

Methods of estimating selection parameters for deleterious mutations that assume a distribution of selection coefficients make inferences about the distribution for new mutations prior to the action of selection, on the basis of the properties of mutations that are segregating in populations; these have been exposed to a long history of selection. The arithmetic mean s for segregating mutations is obtained by summing the products of the selection coefficient at each site i by the corresponding frequency of heterozygotes

and normalizing by the summed frequencies of heterozygotes (⁠

\({\bar{s}}_{\mathrm{seg}}{=}{\sum}s_{i}\mathrm{{\pi}}_{\mathrm{A}_{i}}/{\sum}\mathrm{{\pi}}_{\mathrm{A}_{i}}\)

⁠). For strong selection, Equation 2a implies that this is equal to the harmonic mean of the distribution of s-values for new mutations (Orr and Kim 1998; Sunyaev  et al. 2001).

Tables 2 and 3 show that the mean of the prior gamma distribution can be much larger than the harmonic mean for mutations above the effective neutrality threshold, indicating that much of the probability mass of the gamma distribution is far from the s-values representative of segregating mutations. This raises a serious issue concerning the meaning of inferences concerning the parameters of the gamma distribution; these are based on the properties of mutations that have little relation to the bulk of the mutations in the prior distribution.

With this caveat in mind, one of the strongest results from our arbitrary purifying selection model is not the exact set of parameter values themselves, but rather the exclusion of the large number of parameter combinations that are not compatible with the data. Our difficulties fitting distributions of mutational effects for c_n ≥ 5%, for example, probably suggest that <5% of all nonsynonymous mutations stem from a small set of mutational effects distinct from the continuous distribution, which behave as neutral. Similarly, Tables 2 and 3 allow us to restrict the credible range for the shape parameter of a gamma distribution to shapes >0.1 and usually <1, where 1 is equivalent to an exponential distribution; not many credible estimates have less leptokurtic shapes. This agrees with results for mitochondrial genes, based on a different approach (Piganeau and Eyre-Walker 2003).

If our estimates of the arithmetic and harmonic means of the distribution of s are even approximately correct (of the order of 10⁻⁴ and 10⁻⁵, respectively; see results), they imply that most deleterious nucleotide substitutions affecting protein sequences in our two species are subject to very weak selection. This seems inconsistent with the classical estimates of harmonic mean heterozygous selection coefficients of the order of 1% in D. melanogaster, obtained by comparing inbreeding loads for viability with the mutational decline in mean viability, as well as with estimates of mean homozygous selection coefficients of the order of ≥10% obtained from mutation-accumulation experiments (Crow 1993; Charlesworth and Hughes 2000; Charlesworth  et al. 2004). The former are, however, biased upward by being weighted by the selection coefficients themselves, and the latter are known to be biased upward when there is a wide distribution of homozygous selection coefficients (Crow 1993). In addition, it is likely that insertional mutations that effectively knock out gene function, like those caused by transposable elements, contribute substantially to these estimates (Keightley and Eyre-Walker 1999; Charlesworth  et al. 2004). In contrast, the estimates of heterozygous selection coefficients of the order of 10⁻³, obtained from a population screen for null alleles at enzyme loci in D. melanogaster by Langley  et al. (1981), are consistent with our estimates, assuming that they represent the effects of loss of function at nonvital loci. As pointed out to us by Allen Orr (A. Orr, personal communication), it is difficult to reconcile the results of the null allele screen with the classical estimates of average selection coefficients for deleterious mutations.

Perspectives for the future:

An obvious way to narrow the confidence intervals is the compilation of data sets with more genes. This does not, however, solve the conceptual problem that more degrees of freedom are needed to estimate more parameters. One way of obtaining additional degrees of freedom is to use shared polymorphisms to obtain an estimate of ancestral N_e (Wakeley and Hey 1997). Since ∼3% of the polymorphisms in D. miranda and D. pseudoobscura seem to be shared by both species (Charlesworth  et al. 2005), this is in principle possible. In the present study, we have ignored this approach, due to limited statistical power (only three suitable loci have been surveyed in both species). Another possibility for estimating additional parameters is to use sets of three species where significant differences can be observed in K_A/K_S as well as in π_A/π_S. This might eventually allow simultaneous estimation of P_a and c_n as well as of the shape and location parameters of the distribution of mutational effects. We also deliberately did not classify amino acid changes according to conservative, radical, etc., to keep the model simple. There is no reason why this could not be done with larger data sets, by dividing observations of π_A into several classes according to an a priori classification of their effects on protein function (Sunyaev  et al. 2001; Williamson  et al. 2005).

APPENDIX: EXTENDING THE STRONG PURIFYING SELECTION MODEL

A similar relation can be written for species 2 (the species with larger N_e), retaining the same values of s′ and ψ(s). The only change is that

is substituted for

\(N_{\mathrm{e}_{1}}\)

when specifying γ inside the integral, and c_ne in the first term on the right-hand side of the equation is replaced by c_neα, where α < 1. The α parameter reflects the fact that a larger fraction of mutations with s < s′ do not behave as effectively neutral in species 2, so that s′ does not constitute the border of effective neutrality in this species.

This approach provides a conservative method for improving on the assumption of strong purifying selection, without having to make specific assumptions about the distribution of mutational effects. Application of this method to the weighted data on D. miranda and D. pseudoobscura, assuming γ₁′ = 2 and a mutational bias of 2, gave estimates of c_ne of 4.5% (−1.6%/11.6%),

of 3.00 (1.76/16.7), and an estimate of P_a of 51% (−60%/117%) for D. miranda (the terms in parentheses give the approximate lower and upper bootstrap fifth percentiles). The assumptions used to derive (A2) mean that this approach cannot be used for the species with the larger N_e.

Use of Equation 3a of the text provides an upper-bound estimate of c_ne, since it assumes that all deleterious mutations outside the effectively neutral range have s-values sufficiently large that the deterministic expression in Equation 2a is valid, ignoring mutations with very small s-values. 1/s_h in Equation 2a must therefore be smaller than the values allowing for a wider distribution of s-values, and so c_n must be larger. The true value of c_ne is thus likely to lie between the estimates from Equations 3a and A2.

Given the uncertainties involved, the degree of concordance between the estimates of

from this approximate method and those in Tables 1–3 is encouraging, supporting the conclusion that there must be sufficiently strong selection against deleterious amino acid substitutions for the harmonic mean of N_es to be substantially >1, even in D. miranda with its relatively low effective population size. Unfortunately, the percentile intervals for P_a are so wide that no confidence can be placed in the relevant estimates.

An even more conservative lower bound on

for nonsynonymous mutations above the threshold of effective neutrality is given by setting c_ne to zero in expression (A3); this has the advantage of using polymorphism data on only one species, which makes it widely applicable. With a neutrality threshold of γ₁ = 2 and with κ = 2, I₁ = 0.787, so we get

\(N_{\mathrm{e}_{1}}s_{\mathrm{h}}\)

≥ 2.15 (⁠

\(N_{\mathrm{e}_{2}}s_{\mathrm{h}}\)

≥ 12.5) in the present case.

This estimate can be applied to any suitable data set on coding sequence polymorphisms, as long as a credible assumption about c_ne can be made. For example, Sunyaev  et al. (2000) reported an estimate of 0.33 for the ratio of diversity at nondegenerate coding sites to fourfold degenerate sites, in a large-scale survey of EST-based human SNPs. With κ = 2 and I = 0.787, we get N_es_h ≥ 1.18 for c_ne = 0. With an effective population size for humans of ∼10,000, this suggests a harmonic mean selection coefficient for nonneutral amino acid variants of at least 1.18 × 10⁻⁴. Other studies have suggested that ∼20% of amino acid variants in humans are effectively neutral (Fay  et al. 2001; Sunyaev  et al. 2001); if this value of c_ne is used in expression (A3), the estimate of N_es_h for nonneutral variants becomes 2.42.

Footnotes

Footnotes

Acknowledgement

We thank Deborah Charlesworth, Stephen Schaeffer, and two anonymous reviewers for comments on the manuscript. This work was supported by grants from the Biotechnology and Biological Sciences Research Council, the Leverhulme Trust, and the Royal Society.

References