## Abstract

We develop a Poisson random-field model of polymorphism and divergence that allows arbitrary dominance relations in a diploid context. This model provides a maximum-likelihood framework for estimating both selection and dominance parameters of new mutations using information on the frequency spectrum of sequence polymorphisms. This is the first DNA sequence-based estimator of the dominance parameter. Our model also leads to a likelihood-ratio test for distinguishing nongenic from genic selection; simulations indicate that this test is quite powerful when a large number of segregating sites are available. We also use simulations to explore the bias in selection parameter estimates caused by unacknowledged dominance relations. When inference is based on the frequency spectrum of polymorphisms, genic selection estimates of the selection parameter can be very strongly biased even for minor deviations from the genic selection model. Surprisingly, however, when inference is based on polymorphism and divergence (McDonald-Kreitman) data, genic selection estimates of the selection parameter are nearly unbiased, even for completely dominant or recessive mutations. Further, we find that weak overdominant selection can increase, rather than decrease, the substitution rate relative to levels of polymorphism. This nonintuitive result has major implications for the interpretation of several popular tests of neutrality.

CHARACTERIZING the various forces that shape patterns of genetic polymorphism within and between species is the central goal of population genetics (Lewontin 1974). To that end, statistical inference using Poisson random field (PRF) models (Sawyer and Hartl 1992; Hartl* et al.* 1994; Bustamante* et al.* 2001) provides powerful likelihood and Bayesian methods for quantifying some of these forces, such as mutation and directional selection. Because PRF models assume high levels of recombination between sites, they are particularly well suited to the analysis of polymorphism and divergence at multiple loci distributed across a genome. For example, using a set of sequences from 12 genes in Arabidopsis and 34 genes in Drosophila, Bustamante* et al.* (2002) demonstrated that amino acid substitutions in Drosophila tended to be more advantageous than amino acid substitutions in Arabidopsis, which they attributed to the very high rate of selfing in Arabidopsis. PRF models also provide extremely efficient methods for the simulation of polymorphism and divergence data under the assumption of free recombination, and it has been used in this respect to estimate the power of several statistical tests of neutrality (Akashi 1999).

In developing the original PRF model, Sawyer and Hartl (1992) made some fairly restrictive assumptions regarding population processes. They assumed equal mutational fitness effects, random mating, genic selection (*i.e.*, no dominance), independence among sites (*i.e.*, free recombination), and a stationary population size. Recently, an effort has been made to relax these assumptions one by one. Bustamante* et al.* (2003) relaxed the assumption of equal mutational effects by assuming that the fitness effects of different classes of new mutations are drawn from some underlying distribution. Using this method, it is possible to estimate the distribution of mutational effects using DNA sequence data; the exact form of this distribution is critical to a general understanding of microevolution. Wakeley (2003) relaxed the assumption of random mating by applying an island model of population structure (Wright 1931). He found that, depending on the sampling regime among demes, population structure can strongly bias estimates of mutation rates and divergence times obtained using the basic PRF model if sampling among demes is not accounted for. Surprisingly, however, he found that the estimate of the scaled selection parameter is largely unaffected by island-model population structure.

In this article, we relax the assumption of genic selection. We generalize the original PRF model to account for arbitrary diploid selection models. Our generalized model yields maximum-likelihood estimators (MLEs) for the scaled selection coefficient, γ, and the dominance parameter, *h*. This is the first DNA sequenced-based estimator of the dominance parameter, and simulations reveal that it performs very well when large numbers of segregating sites are available. We also use our generalized PRF model to investigate the impact of dominance on polymorphism and divergence data. Two surprising results emerge from this analysis. First, we find that dominance relations generally have very little impact on ratios of polymorphism to divergence. Therefore, statistical inferences based on these ratios (Sawyer and Hartl 1992) are robust to a violation of the assumption of genic selection. Second, we find that weak and moderate overdominant selection can sometimes decrease the ratio of polymorphism to divergence relative to a neutral standard. This result is contrary to intuition. One might expect that balancing selection (overdominance is a special case of balancing selection) would always increase, rather than decrease, the ratio of polymorphism to divergence because balancing selection actively maintains polymorphism.

A wealth of theory exists relating how different types of dominance relations affect patterns of polymorphism and divergence (*e.g.*, Cherry 2003, 2004; Griffiths 2003; Roze and Rousset 2003; Whitlock 2003). Noticeably absent, however, are any methods of statistical inference for distinguishing nongenic from genic selection, so the role of dominance has been relatively understudied by empirical population geneticists. The dominance parameter has important implications for a number of evolutionary phenomena, such as inbreeding depression (Lynch and Walsh 1998, Chap. 10), the maintenance of genetic variation by mutation-selection-drift balance (Charlesworth and Hughes 2000), and genetic load (Crow 1993). The methods we present here open up the possibility of using large DNA sequence and single-nucleotide polymorphism (SNP) data sets to investigate how dominance affects variation at the molecular level.

## THEORY

For a given site, consider the case of irreversible mutation from an ancestral nucleotide *A*_{1} to a derived nucleotide *A*_{2}, occurring at rate μ. When this mutation process is applied across many sites, it corresponds to the infinite-sites mutation model (Kimura 1968, 1971; Watterson 1975). Let 1, 1 + 2*sh*, and 1 + 2*s* be the relative fitnesses of the genotypes *A*_{1}*A*_{1}, *A*_{1}*A*_{2}, and *A*_{2}*A*_{2}, respectively. With random mating, this model is formally equivalent to models of frequency-dependent selection where fitness is linearly related to allele frequency (*e.g.*, Cherry 2004). Under the usual assumptions of the Wright-Fisher model (random mating, constant population size, nonoverlapping generations), Wright (1938) derived the quasi-stationary distribution of allele frequency for the above diploid selection scheme. Kimura (1964) later derived a more concise form, 1where *q* is the frequency of the derived nucleotide, and γ = 2*Ns*. Note that this expression is equivalent to Griffiths' (2003) expression (32) for the distribution of allele frequencies under stationarity, which was derived in a different manner. This can be demonstrated by substituting μ(*x*) = 4γ(*h* + (1 − 2*h*)*x*)*x*(1 − *x*) and σ^{2}(*x*) = *x*(1 − *x*) into expressions (12), (13), (14), and (32) in Griffiths (2003).

Expanding to multiple sites, if each site is independent, then the mutant allele frequency at each site is a random draw from the above distribution, with the relative density of the distribution proportional to the mutation rate (Sawyer and Hartl 1992). To estimate the parameters of the above distribution, consider polymorphism data in the form of a site frequency spectrum. If outgroup sequence data are available so that ancestral and derived nucleotides can be distinguished, then the site frequency spectrum is a vector, **x**, where each entry, *x _{i}*, is a count of the number of sites at which the derived nucleotide is represented

*i*times in a sample of size

*n*, for

*i*= 1, 2, … ,

*n*− 1. For a given allele frequency,

*q*, at a given site, the probability of choosing

*i*derived nucleotides in a sample of

*n*individuals is given by a binomial distribution with mean

*nq*. Thus, the expectation of each of the

*x*is θ

_{i}*F*(

*n, i*; γ,

*h*), where 2θ = 4

*N*μ, and μ is the per-generation mutation rate in the entire region sampled. Further, if a Poisson number of mutations enter the population each generation, then each of the

*x*will be Poisson distributed (Ewens 1974; Sawyer and Hartl 1992). With the full probability distribution of each entry of the site frequency spectrum, the model parameters θ,

_{i}*h*, and γ can be estimated using maximum-likelihood methods.

Kimura (1964) also derived the fixation rate under the dominance model. Measuring time in 2*N* generations, the instantaneous fixation rate at stationarity is 3

Note that *u*(γ, *h*) is equal to the stationary distribution (1), evaluated at *q* = 1. PRF theory predicts that the number of fixations over an arbitrary length of time, τ, will be Poisson distributed with mean 4

For example, if two species diverged *t*_{div} generations ago, and both species have the same population size *N*, then the expected number of fixed differences observed in a sample from the two species is 5where *n*_{1} and *n*_{2} are the sample sizes in the two populations.

The effects of various types of nonadditive allelic interactions on the stationary distribution of allele frequency are shown in Figure 1. Also shown are the expected site frequency spectra for different combinations of the selection (γ) and dominance (*h*) parameters. Note that *h* = 0.5 corresponds to genic selection. For negative selection (γ < 0), the density of the stationary distribution, and hence the expected number of segregating sites in the site frequency spectrum, is negatively related to *h*. This result is straightforward: the more recessive the deleterious nucleotide is, the more likely it is to drift to observable frequencies. For positive selection (γ > 0), the effect of dominance on the stationary distribution is more subtle. For high allele frequencies, the density of the stationary distribution is positively related to *h*. This is due to the fact that, once a dominant, advantageous mutation has attained high frequency (*e.g.*, *q* = 0.9), the sojourn time to fixation will be relatively long because the ancestral homozygote genotype will be rare. More surprisingly, at very low frequencies the dominance parameter does not strongly affect the density of the stationary distribution. One might expect to observe a relative excess of recessive, advantageous mutations at low frequencies because additive and dominant mutations attain high frequency much faster. However, this force is apparently counterbalanced by the greater probability of stochastic loss of the recessive mutations. In the case of heterozygote advantage (*h* > 1 for γ > 0, or *h* < 0 for γ < 0), there is sometimes an interior mode in the stationary distribution. This mode is centered on the deterministic prediction for the stable equilibrium allele frequency (Fisher 1922, 1930; Haldane 1926; Wright 1931), which, using our parameterization, occurs at *q̃* = (*h* − 1)/(2*h* − 1).

### Maximum-likelihood estimation conditioning on the number of segregating sites:

Because , and because each of the *x _{i}* are Poisson distributed, the joint probability of the

*x*conditional on

_{i}*S*is given by the multinomial distribution with

*n*− 1 frequency classes, and the probability of each class is 6

The denominator of the expression (5) is the expectation of *S*. Also note that the mutation parameter will cancel. The log-likelihood of a given site frequency spectrum is then 7

Maximum-likelihood estimates of γ and *h* can be obtained by finding the maximum of (7) using standard optimization techniques. We investigate the sampling properties of these MLEs using both asymptotic-likelihood theory and stochastic simulations.

Asymptotic-likelihood theory predicts that, for large sample sizes (in our case the number of segregating sites), the joint sampling distribution of our MLEs will be multinormally distributed (Kendall and Stuart 1973), with means given by the underlying true values of γ and *h* and a variance-covariance matrix given by the inverse of the Fisher information matrix * I*: 8

Thus the direct way to predict the sampling variance and covariance is to calculate the second derivatives of the log-likelihood function and evaluate their expectations. Let *F _{i}* be shorthand for

*F*(

*n, i*; γ,

*h*). The first derivatives of the log-likelihood function in γ and

*h*are 9a9bexchanging the order of integration and the derivative, we have 10a10bwhere the first derivatives of

*f*are given in the appendix. The second derivatives of the log-likelihood function are 11a11b11cwhere 12a12b12cand 13a13b13c

The expressions for the second derivatives of *f* are also given in the appendix. To arrive at the entries of the Fisher information matrix, we evaluate the expectations of second derivatives of ℓ, noting that the *T* terms do not depend on *x _{i}*, 14a14b14cwhere . The inverse of

*is 15*

**I**Therefore, the sampling variances and covariances of the MLEs for the selection parameter, γ, and the dominance parameter, *h*, are 16a16b16c

When investigating the sampling properties of the MLEs, we utilize assumed underlying “true” values of γ and *h* in evaluating *E*[*x _{i}*] and the various

*T*

_{- -}terms. In practice, one may approximate the sampling variances and covariance of the parameter estimates by substituting the MLEs into expressions (16).

### Sampling properties of γ̂ and *ĥ*:

The asymptotic joint sampling distribution of γ̂ and *ĥ* is shown in Figure 2 for a reasonable sample size (*n* = 25), a large number of segregating sites (*S* = 10,000), and several underlying true values of γ and *h*. Also shown are the simulated joint sampling distributions (see below for details of the simulations). In general, the estimation procedure seems to perform exceptionally well when data are gathered from a large number of segregating sites: The entries of the variance-covariance matrix are small, and simulated MLEs cluster tightly around true values. Also, the agreement between the asymptotic prediction for the joint sampling distribution of our MLEs and the simulated distribution is generally quite good. These results indicate that it is possible to estimate dominance parameters from DNA polymorphism data alone. The main exception to this result is the case of strongly deleterious (γ < −5) and at least partially dominant (*h* > 0.5) mutations. In this region of the parameter space, the sampling variances and covariance of γ̂ and *ĥ* become extremely large. For instance, for γ = −20 and *h* = 0.8, the asymptotic prediction for the 95% confidence interval on *ĥ* is 0.8 ± 4.11, and the 95% confidence interval on γ̂ is −10 ± 98.71. In this situation, one will have virtually no power to make inferences or to reject null hypotheses such as γ = 0 (neutrality) or *h* = 0.5 (genic selection). Fortunately, though, quantitative genetic analyses (*e.g.*, Simmons and Crow 1977; Crow and Simmons 1983; Willis 1999; Kelly 2003) suggest that strongly deleterious mutations tend to be recessive. Also, biochemical models (Kacser and Burns 1981) suggest that mutations of large effect will tend to be recessive. To summarize, the one situation where our estimator performs poorly is thought to occur rarely in natural populations.

Thus far we have reported results for a large number of segregating sites, with the presumption that the method would be applied at the genomic level to large SNP data sets, and we have found that the estimation procedure works surprisingly well in this situation for most parameter combinations. For smaller data sets (*S* < 100), our ability to simultaneously estimate the selection and dominance parameters is greatly diminished. Noting that *E*[*x _{i}*] ∝

*S*, we see that the sampling variance and covariance terms in expressions (15) are inversely proportional to

*S*. This hyperbolic dependence on

*S*is reflected in Figure 3, which shows the standard deviation of the MLEs as a function of

*S*.

### Power to detect nonadditive allelic interactions:

The two-parameter model with selection and dominance can be compared to the one-parameter basic PRF model (Sawyer and Hartl 1992) by employing a likelihood-ratio test (LRT). Here, genic selection (*h* = 0.5) is the null hypothesis and the LRT statistic, 2(ℓ(γ̂, *ĥ*|**x**) − ℓ(γ̂_{o}, ^{1}/_{2}|**x**)), is expected to be chi-square distributed with 1 d.f. under the asymptotic assumption (γ̂_{o} denotes the maximum-likelihood estimate of the selection parameter under the assumption of genic selection). We conducted stochastic simulations to determine the statistical power (probability of rejecting genic selection) this test has in detecting deviations from strictly genic selection. The simulation procedure is straightforward. First, to simulate data, we conducted *S* pseudo-random draws from a multinomial distribution with *n* − 1 classes, where the probability of each class is given by expression (6). Using this simulated data, the likelihood function (7) was evaluated using the extended midpoint numerical integration algorithm (Press* et al.* 1988, Chap. 4) and then maximized using the Fletcher-Reeves-Polak-Ribiere optimization routine (Press* et al.* 1988, Chap. 10). The likelihoods from the two different models were then compared using the LRT statistic. This procedure was repeated 1000 times for each parameter combination, using the 190-node Computational Biology Service Unit cluster at the Cornell Theory Center (www.tc.cornell.edu).

Statistical power, *i.e*., the proportion of tests that reject the null hypothesis, is shown in Figure 4 as a function of the dominance parameter *h*. Results are not reported for dominant, strongly deleterious mutations (γ = 20 and *h* > 0.5) because the large sampling variances of the MLEs in this region of the parameter space made it difficult to optimize the likelihood function. In general, given enough segregating sites, the LRT is extremely powerful in detecting deviations from the genic selection model, even for very weak selection (|γ| < 5) and incomplete dominance (0 < *h* < 1, *h* ≠ 0.5). This indicates that the “signature” of nongenic selection is evident from patterns of DNA sequence polymorphism, and it is indeed possible to identify nongenic selection and estimate dominance parameters. Figure 4 also shows that, above *n* = 25, the sample size makes very little difference in one's ability to detect departures from genic selection. However, an increase in the observed number of segregating sites can make a substantial difference in statistical power. Despite the positive relationship between power and the number of segregating sites, we do observe appreciable power even for a small number of segregating sites in some situations. For example, if advantageous mutations are completely dominant or completely recessive, statistical power is high even for *S* = 100.

In estimating statistical power, we have made the asymptotic assumption that our LRT statistic is chi-square distributed. To assess the validity of this assumption, we conducted several simulations under the null hypothesis of genic selection. If the asymptotic assumption is appropriate, then the 95% quantile of our simulated LRT statistics should closely approximate the 95% (*P* = 0.05) critical value of the chi-square distribution with 1 d.f. The 95% quantiles of our simulated LRTs are shown in Table 1 for several combinations of the selection parameter, sample size, and number of segregating sites. In general, the simulated critical values are very close to the chi-square critical value. Also, the simulated distribution of the LRT statistic closely conforms to the chi-square distribution with 1 d.f. Some example distributions are shown in Figure 5.

### Bias in estimating selection parameters:

To assess the bias caused by estimating selection in the presence of unacknowledged dominance relations, we simulated data for various degrees of dominance and then estimated the selection parameter assuming the one-parameter genic-selection model. The relative bias due to unacknowledged dominance relations is shown in Figure 6. Dominance can have a major impact on the genic-selection estimate of the selection parameter, especially for the case of strong positive selection. If new mutations are advantageous but recessive, then the genic-selection model substantially underestimates the selection parameter, to the point that one sometimes incorrectly identifies negative, rather than positive, selection. If new mutations are advantageous and dominant, then the genic-selection model yields large to extreme overestimates of the selection parameter, even for slightly dominant mutations (0.5 < *h* < 0.7). Our genic-selection estimates of the selection parameter were sometimes pushed to the upper limit (γ = 100) of the range of possible selection parameters that we allowed in our simulations. In general, the genic-selection model for the site-frequency spectrum does not reliably characterize nongenic, positive selection. For the case of negative selection, the bias introduced by dominance relations is less extreme, but still substantial. In general, if new mutations are deleterious and at least partially recessive, then one tends to overestimate the selection parameter. The opposite pattern is observed if new mutations are deleterious and partially dominant; *i.e.*, one tends to underestimate the selection parameter.

Given the sometimes extreme bias in the genic-selection estimate of the selection parameter, we require some criteria to determine when it is appropriate to apply the genic-selection model to data in the form of a site-frequency spectrum. A simple chi-square goodness-of-fit test is appropriate for this purpose. Under genic selection, the following statistic should be chi-square distributed with *n* − 3 d.f., 17where 18

In our simulations, this goodness-of-fit test leads us to convincingly reject the genic-selection model in cases where the genic-selection estimate of the selection parameter is strongly biased. For example, for γ = 20 and *h* = 0.7, the average MLE of γ under the genic-selection model is γ̂_{o} = 66.9. In this case, the average goodness-of-fit statistic over 1000 simulations has a *P*-value <10^{−100}.

### McDonald-Kreitman polymorphism and divergence data:

The McDonald-Kreitman test of neutrality (McDonald and Kreitman 1991) contrasts the ratios of polymorphism to divergence across different classes of mutations. If one of the classes (*e.g.*, synonymous sites) is thought to evolve neutrally for *a priori* reasons, this class can be used as a “neutral standard,” and the ratio of polymorphism to divergence from the potentially selected class (*e.g.*, nonsynonymous sites) can be compared to this standard to detect the action of selection. Because the McDonald-Kreitman test is based on an observed standard for neutral evolution in the particular population in question, it is fairly robust to demographic deviations from the equilibrium neutral model, such as population subdivision or fluctuating population size (McDonald and Kreitman 1991; Akashi 1999; Nielsen 2001). Therefore, results of a McDonald-Kreitman test are often easier to interpret than those of many other neutrality tests, such as Tajima's (1989) *D*-test and Fu and Li's (1993) series of tests, which are sensitive to demographic, as well as selective, forces (*e.g.*, Golding 1997).

The original McDonald-Kreitman test (McDonald and Kreitman 1991) was devised as a 2 × 2 contingency table analysis. As an alternative, Sawyer and Hartl (1992) developed a maximum-likelihood framework for analyzing polymorphism and divergence, which allows the estimation of selection, mutation, and divergence time parameters, as well as hypothesis testing. Let *S*_{n}, *S*_{s}, *D*_{n}, and *D*_{s} be the observed numbers of nonsynonymous polymorphisms, synonymous polymorphisms, nonsynonymous fixed differences, and synonymous fixed differences, respectively (hereafter referred to as McDonald-Kreitman data). With the usual assumptions of the PRF model (see above), one can obtain parameter estimates by setting the observed values to their expectations under the PRF model and then solving for the parameters. These estimates, derived using the method of moments, are equivalent to the maximum-likelihood estimates (Sawyer and Hartl 1992). Note that it is not possible to estimate additional parameters, such as the dominance parameter, from these types of data because, under genic selection, there are four equations and four unknown parameters: the nonsynonymous and synonymous nonlethal mutation rates, the divergence time, and the selection parameter. Therefore, without prior knowledge of some of the parameters, McDonald-Kreitman data contain no additional information regarding dominance.

Even though it is not possible to estimate dominance parameters from McDonald-Kreitman data, we can still investigate how dominance relations affect parameter estimates obtained using McDonald-Kreitman data under the assumption of genic selection. This is important because, even though the McDonald-Kreitman test and Sawyer and Hartl's parametric methods are thought to be robust to deviations from an ideal population model (*e.g.*, nonrandom mating, nonstationary population size), it is not known how sensitive this approach is to deviations from the assumed form of selection. Using fixed values for the mutation rate and divergence time (the divergence time is usually estimated from synonymous sites, and the mutation rate does not need to be estimated to estimate the selection parameter), we simulated nonsynonymous polymorphism and divergence data for several values of the selection and dominance parameters and then estimated the selection parameter under genic selection. For each iteration, the number of nonsynonymous segregating sites was drawn from a Poisson distribution with mean 19and the number of nonsynonymous fixed differences was drawn from a Poisson distribution with mean 20

Given these simulated data, we estimated the selection parameter, γ, under genic selection by numerically solving Equation 22 in Sawyer and Hartl (1992).

The bias in Sawyer and Hartl's (1992) estimate of the selection parameter is shown in Figure 7 as a function of the degree of dominance. Surprisingly, dominance relations, including the case of weak overdominance, have very little impact on estimates of γ obtained assuming genic selection. This is due to the fact that dominance has a similar effect on both polymorphism and divergence. Conditional on the scaled divergence time, τ, which is generally estimated from synonymous sites, Sawyer and Hartl's estimate of the selection parameter depends solely on the ratio of nonsynonymous polymorphism to divergence. Therefore, if dominance has a similar effect on polymorphism and divergence, then it will not appreciably affect estimates of the selection parameter. For instance, in the case of positive selection, an increase in the dominance parameter causes both an increase in the level of polymorphism and an increase in the fixation rate (Figure 8). If these effects roughly cancel out in a ratio of polymorphism to divergence, then Sawyer and Hartl's method would be insensitive to nongenic selection. The ratio of the expected number of polymorphisms to the expected number of fixed differences is shown in Figure 9 for several different values of dominance parameter, *h*. Dominance relations have very little impact on the ratio of polymorphism to divergence, which explains why we observe so little bias in the genic-selection estimate. Wakeley (2003) recently demonstrated that, when applied to McDonald-Kreitman data, Sawyer and Hartl's method is robust to the assumption of random mating by applying an island model of population structure. Using simulations, Weinreich and Rand (2000) demonstrated that McDonald-Kreitman ratios—and, consequently, Sawyer and Hartl's estimate of the selection parameter—are not sensitive to dominance relations for a limited range of the dominance parameter (0 ≤ *h* ≤ 1 in our notation, 0 ≤ *h* ≤ 2 in their notation). Our results indicate that this is also true of weakly overdominant and underdominant mutations.

The effect of heterozygote advantage (*h* > 1 for γ > 0, or *h* < 0 for γ < 0) on McDonald-Kreitman data deserves special attention. Heterozygote advantage is a special case of balancing selection, and its effect on polymorphism and divergence approximates the effect of several other types of balancing selection (Wright and Dobzhansky 1948; Denniston and Crow 1990; Takahata and Nei 1990). One might expect that balancing selection should cause an increase in the ratio of nonsynonymous polymorphism to divergence because balancing selection actively maintains polymorphisms. In our results from the case of heterozygote advantage, we have assumed the best-case scenario for detecting balancing selection: *every* new mutation is independent and subject to heterozygote advantage, yet we see that the ratio of polymorphism to divergence is actually sometimes less than the neutral standard. This occurs when the derived homozygote is more fit than the ancestral homozygote and the degree of heterozygote advantage is fairly weak (*h* = 1.3, γ > 0 in Figure 9). As the strength of heterozygote advantage increases (*i.e*., γ and *h* increase) and allele frequencies are more tightly maintained at intermediate frequencies, this trend reverses because the level of polymorphism increases faster than the fixation rate (*h* = 2, γ > 6 in Figure 9). However, we submit that McDonald-Kreitman tables provide information regarding balancing selection that is ambiguous at best—if heterozygote advantage is weak or moderate, then significant results could be interpreted as positive directional selection rather than balancing selection.

## DISCUSSION

Dominance plays a very important role in a number of evolutionary phenomena at the heart of population genetics. For instance, a central controversy in evolutionary genetics has been the dispute over whether a balance between deleterious mutations and purifying selection can account for the bulk of genetic variation in fitness-related traits (*e.g.*, Lewontin 1974). Alternative explanations suggest that selection actively maintains genetic variation via frequency-dependent selection, variation in selection intensity over time and space, heterozygote advantage, or other higher-order processes. Theoretical models of deleterious mutation/purifying selection balance depend primarily on the deleterious mutation rate and the product of the selection and dominance parameters (*e.g.*, Charlesworth and Hughes 2000). Therefore, assessing the validity of the deleterious mutation hypothesis will require characterizing both selection and dominance at the level of the entire genome.

In the fields of molecular population genetics and molecular evolution, there has been a growing consensus that weak negative selection plays an important role in evolution (*e.g.*, Ohta 1992). Virtually all of the arguments in support of this finding are based on the assumption of genic selection. However, if moderately or strongly deleterious mutations tend to be recessive, they could “appear” to be weakly selected on the basis of their frequency profile. For instance, we have shown that nongenic selection can strongly bias estimates of the selection parameter that are based on the site-frequency spectrum. To address whether the apparent signal of weak purifying selection is actually due to more deleterious, but recessive, mutations, it will be necessary to quantify both selection and dominance parameters. The analyses presented here are a first step toward that goal.

At the level of the phenotype, two independent lines of evidence suggest that dominance plays an important role in shaping polymorphism and divergence in fitness-related traits. The first line of evidence is the simple observation of inbreeding depression, *i.e.*, the reduction in offspring fitness due to inbreeding. The two leading explanations of inbreeding depression are: (1) Fitness is reduced because inbreeding “unmasks” deleterious recessives, and (2) fitness is reduced because fewer overdominant loci occur in the heterozygous state (see Lynch and Walsh 1998, Chap. 10 for an overview of the two hypotheses). Both explanations are based on dominance relations among segregating, nonneutral alleles at many loci. The second line of evidence comes from speciation genetics: Growing consensus supports the dominance theory (Muller 1942; Turelli and Orr 1995) of Haldane's rule (Haldane 1922; Orr 1997). Haldane's rule is an observation from studies of postzygotic isolation between species with chromosomal sex determination. It states that, in hybrid crosses between species, if only one sex is inviable or infertile, it will be the heterogametic sex. Recent analyses (Turelli and Begun 1997; Presgraves and Orr 1998) strongly support the hypothesis that Haldane's rule is due to dominance relations in sex-chromosome-linked “speciation genes,” *i.e.*, genes at which fixed differences between lineages contribute to postzygotic isolation. Given the importance of dominance in shaping both polymorphism and divergence at the phenotypic level, and given the implications of dominance for evolutionary theory, there is an urgent need to quantify patterns of nongenic selection at the molecular level.

In this article we have made the somewhat artificial assumption that all new mutations have the same fitness effect. An obvious extension to this method is to allow variable fitness and dominance effects by assuming that the effects of each new mutation are drawn from some bivariate distribution and then integrating over all possible fitness and dominance effects to arrive at predictions for the stationary distribution of allele frequency. This approach is computationally challenging because it requires numerically evaluating complicated four-dimensional integrals during each evaluation of the likelihood function. We have also assumed high levels of recombination between sites. This assumption may be appropriate for SNPs distributed across a genome, but it does not apply to small regions with limited or no recombination such as single protein-coding genes or animal mitochondrial genomes. It may be possible to address this problem using a composite-likelihood approach (Hudson 2001).

## APPENDIX

To evaluate the first and second derivatives of the stationary distribution, *f*, we apply the product rule from calculus several times and exchange the order of integration and derivation when necessary. The first derivatives in γ and *h* are then A1aA1bwhere A2aA2b

And the second derivatives are

A3a

A3b

A3c

## Acknowledgments

We thank J. F. Crow, J. K. Kelly, and R. Nielsen for very helpful discussions and suggestions. Two anonymous reviewers greatly improved the manuscript. This research was conducted using the resources of the Cornell Theory Center, which receives funding from Cornell University, New York State, federal agencies, foundations, and corporate partners. C.D.B. is supported by a U.S. Department of Agriculture-HATCH grant (NYC-151411) and National Science Foundation grant 0319553.

## Footnotes

Communicating editor: J. B. Walsh

- Received November 19, 2003.
- Accepted May 15, 2004.

- Genetics Society of America