## 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^{−5} 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*_{e}*s*, 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*_{e}*s* estimates predominantly between 1 and 3; Bustamante *et al.* (2002) estimated *N*_{e}*s* 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*_{e}*s* values across different genes, was applied to data on *D. melanogaster*, indicating a mean *N*_{e}*s* 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*_{e}*s*.

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*_{e}*s* value or a distribution of *N*_{e}*s* values across loci is assumed.

These studies have yielded large differences among estimates of the scaled selection parameter *N*_{e}*s* (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*_{e}*s* 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^{−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^{−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 fit of the nucleotide site diversity from the two species to a given set of assumed parameters was assessed, using the numerical criterionwhere is the predicted nonsynonymous diversity value for species 1, is the corresponding observation, and the other corresponding values are for silent diversity and for species 2, using Equation 5a below.

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^{−9} 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*_{a}*u*, 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.

These assumptions lead to the following general equations for the expectations of the silent and nonsynonymous nucleotide site diversities in species *i*, and , and the rates of substitution per site per generation for silent and nonsynonymous mutations, and (second-order terms in small quantities are ignored):(1a)(1b)(1c)(1d)where is the mean equilibrium diversity at sites subject to purifying selection, and is the mean substitution rate at such sites.

#### Strong purifying selection:

These equations become greatly simplified if *N*_{e}*s* > 1 for all nonneutral nonsynonymous mutations in both species. The equilibrium diversity contributed by sites subject to purifying selection with selection coefficient *s* is then well approximated by the deterministic expression, 2*u*/*s*, as can be shown by numerical solutions of the general equations (McVean and Charlesworth 1999). We then have(2a)where θ_{i} = , and *s*_{h} is the harmonic mean of the selection coefficients for all mutations that are not effectively neutral (this is the same for both species, since we assume that nearly neutral mutations are absent).

For *N*_{e}*s* > 2, is negligibly small, so that we can replace Equation 1d by(2b)With just two species for which diversity data are available, Equations 1a and 2a lead to the following formula for *c*_{n}:(3a)Substituting this into Equation 2a and using Equation 1a, we obtain(3b)From Equations 1c and 2b, *c*_{a} is given by(4a)assuming that species 1 and 2 have the same mean divergences for silent and nonsynonymous sites from the third species, so that subscripts can be dropped. This is necessarily the case with strong selection, since species 1 and 2 are equally close to species 3, and only neutral nonsynonymous mutations can become fixed (there is no reason in principle why divergence between species 1 and 2 could not be used in the absence of data on a third species, but in the present case the level of divergence between *D. miranda* and *D. pseudoboscura* is so low that estimates based on this would be very unreliable).

The proportion of nonsynonymous substitutions that are caused by the fixation of advantageous mutations is(4b)Given the above assumptions, all the parameters of the model can be estimated by equating expectations to observed values.

#### Arbitrary purifying selection:

The validity of the assumption that *N*_{e}*s* > 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 , and so we represent it as 1/κ.

κ is the ratio of the forward and backward mutation rates for mutations creating deleterious amino acids, a measure similar to the mutational bias parameter used in models of codon usage bias (Li 1987; Bulmer 1991; McVean and Charlesworth 1999). At a site fixed for a deleterious amino acid, there will thus be a mutation rate *u*/κ back to the preferred amino acid; there will also be a neutral mutation rate *u*(κ − 1)/κ, if all deleterious amino acids at this site are selectively equivalent. These should be taken into account in the contribution to net diversity. Using the argument that led to Equations 6 and 7 of McVean and Charlesworth (1999), we then find that in Equation 1b is given by(5a)where *N _{i}* is the total number of breeding individuals in species

*i*;

*H*

_{0}(

*s*) and

*H*

_{1}(

*s*) are the expected total heterozygosities that are contributed during their sojourn in the population by new mutations to preferred and unpreferred amino acids, respectively, for a selection coefficient

*s*;

*m*

_{0}(

*s*) and

*m*

_{1}(

*s*) are the fractions of sites fixed for unpreferred and preferred amino acids, respectively, among sites with selection coefficient

*s*; and ϕ(

*s*) is the probability density function for the distribution of selection coefficients. Formulas for the

*H*and

*m*functions are given by McVean and Charlesworth (1999), Equations 5 and 10.

Similarly, in Equation 1d is given by(5b)where *U*_{0}(*s*) and *U*_{1}(*s*) are the fixation probabilities for new mutations to preferred and unpreferred amino acids, respectively, given a selection coefficient *s*, using the standard diffusion equation formula (Kimura 1962).

Even if the distribution ϕ(*s*) is described by only two parameters, such as the arithmetic mean and standard deviation, there are too few degrees of freedom in the data to estimate all the parameters of interest by equating observed and expected values of diversities and divergence, unless we are prepared to assume that there are no nonsynonymous sites with neutral effects (*c*_{n} = 0). These equations do, however, provide a means of evaluating the sensitivity of the results to our assumptions about *c*_{n} or the properties of the distribution, as is described below. Following convention (Piganeau and Eyre-Walker 2003), we use the gamma distribution for ϕ(*s*),(6)where α is the shape parameter, β is the scale parameter, Γ is the gamma function, αβ is the arithmetic mean, and αβ^{2} is the variance of *s* (R-Project 2005).

## 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 is, however, very wide, and infinite values were sometimes generated for the weighted data. The estimate of *N*_{e}*s*_{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*_{e}*s* < 1 or even 0.5. The results of McVean and Charlesworth (1999) suggest that, with *N*_{e}*s* < 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*_{e}*s* = 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*_{e}*s* = 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: , , and , 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*_{e}*s* < 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*_{e}*s* < 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^{−9}. The distributions for the two smaller values of *c*_{n} peaked at selection coefficients ∼10^{−4}, 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*_{e}*s*_{h} and again suggest that the magnitudes of the estimates of the main parameters of interest are fairly robust to details of the assumptions.

To simplify even further, one could assume that all polymorphic nonsynonymous sites are effectively under purifying selection; *i.e.*, *c*_{n} = 0. The mean frequency of such deleterious mutations is given by *q* = *u*/*s*_{h} (see Equation 2a), where *u* is the mutation rate per site per generation and *s*_{h} is the harmonic mean of the heterozygous selection coefficients. Since *q* is small, we have(8a)Using π_{S} = 4*N*_{e}*u*, we obtain(8b)

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*_{e}*s*_{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*_{e}*s*_{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*_{e}*s* 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*_{e}*s* 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*_{e}*s* 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 4*N*_{e}*u* for this species. Accordingly, *N*_{e}*s* 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^{−9} mutations per base per generation happen in *D. melanogaster*. Powell (1997, p.369–371) reported rates between 0.67 × 10^{−9} and 3.3 × 10^{−9}, 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^{−9} (95% C.I. = 1.0 × 10^{−9}–2.5 × 10^{−9}), 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^{−9} and 4.75 × 10^{−9}. Thus 0.5 × 10^{−9} and 8 × 10^{−9} 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 (). 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^{−4} and 10^{−5}, 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^{−3}, 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

We can extend the strong purifying selection model as follows, to lighten the assumptions involved. The aim is to place bounds on the estimates of the fraction of nearly neutral mutations and the harmonic mean of *s*, for the species with the lower effective population size, *i.e.*, species 1. We use the approximate formula for the equilibrium diversity at sites under selection given by Equation 15 of McVean and Charlesworth (1999). A simple extension to this, using the formulation that led to Equation 5a, yields the following expression for species 1,(A1)where γ_{1} = is the scaled measure of selection intensity for species 1, ψ(*s*) is the probability density of *s*, conditional on *s* falling outside the domain of effective neutrality (defined by the relation γ_{1} ≤ = 2), and *c*_{ne} is the fraction of neutral and effectively neutral mutations for species 1 (see Table 2). Numerical integrations have shown that Equation A1 typically predicts diversity patterns with a relative error of ∼5–10% when compared to our more accurate method.

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 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.

For our purposes, the term in braces is just a nuisance parameter, since we are interested only in , the harmonic mean of selection coefficients *s* > *s′* for species 1. Using the mean-value theorem, we can replace the integrals for the two species bywhere is a value of γ_{i} in the domain of integration above *s′*, and *s*_{h} is the harmonic mean of *s* with respect to ψ(*s*) over this domain, *i.e.*, the harmonic mean of *s* for amino acid mutations other than effectively neutral ones in our focal species, species 1.

*I _{i}* is an increasing function of for κ ≥ 0 and ≥ 0, so that a lower bound is obtained by setting to γ

_{i}′ = . The upper bound is 1; in the present case, this is very close to the actual value of

*I*

_{2}and is used in its place. By the same argument that led to Equations 2a and 3a, and using the lower bound of

*I*

_{1}and the upper bound of

*I*

_{2}, together with the fact that α < 1, after some algebra we obtain a lower bound to the estimate of

*c*

_{ne},(A2)where is the estimate of

*c*

_{n}from Equation 3a, and

*r*is the ratio of silent diversity for species 2 to that for species 1. This expression can in turn be used to yield a lower-bound estimate for by using Equations 1a and 1b:(A3)Similarly, approximate estimates of

*c*

_{a}and

*P*

_{a}are obtained by substituting (A2) into Equations 4.

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 γ_{1}′ = 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*_{e}*s* 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 γ_{1} = 2 and with κ = 2, *I*_{1} = 0.787, so we get ≥ 2.15 ( ≥ 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*_{e}*s*_{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^{−4}. 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*_{e}*s*_{h} for nonneutral variants becomes 2.42.

## Acknowledgments

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.

## Footnotes

↵2

*Present address:*Unidade de Xenética Evolutiva, Instituto de Medicina Legal, Facultade de Medicina, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain.Communicating editor: S. W. Schaeffer

- Received June 23, 2005.
- Accepted November 5, 2005.

- Copyright © 2006 by the Genetics Society of America