METHODS

Computer simulations: Throughout, computer simulations were used to test theoretical predictions and to examine the properties of the suggested statistical tests. In these simulations a haploid population of size N chromosomes was simulated under the Wright-Fisher model. To model the mutational process an infinite sites model was used (Watterson 1975). The mutational makeup of each individual chromosome and the generation in which each segregating mutation arose were stored in computer memory. In each generation the population of chromosomes underwent mutation, selection, and genetic drift. The probability of a new mutation occurring in a given generation on a chromosome is μ=μ_s +μ_n, where μ_s is the mutation rate of selected mutations and μ_n is the mutation rate of neutral mutations. Simple multiplicative fitness was employed; thus the relative fitness of the ith chromosome is given by wi=(1+s)ni∕∑j(1+s)nj, (1) where s is the selection coefficient of selected mutations and n_i is the number of selected mutations carried by the ith chromosome. For the sake of simplicity we assume that s does not vary in time or among selected mutations and that there are no epistatic interactions. Selection together with genetic drift is simulated by sampling multinomially between generations among the k segregating chromosomes according to the probabilities specified by the w_i, i = 1, 2,.., k given by (1). In general, samples were taken from the population every 10N generations after the first 100N generations. For each sample, some function of the sample was evaluated, such as the mean age of selected or neutral segregating sites or the P value under one of the suggested tests. Denote the obtained value of the evaluated function in the ith sample by f_i. Assuming that stationary distributions exist for these functions, the theorem of ergodicity ensures that 1n∑i=1nfi→E(f)asn→∞, where E(f) is the expectation of the function f and n is the number of samples. To examine the power of the McDonald/Kreitman test, we also scored the number of fixed differences between the sequences in each sample and a single sequence sampled t generations in the past. Population sizes of N = 1000 chromosomes were used in the simulations throughout. Except where otherwise noted, the parameter values used in the simulations were chosen to fit the estimated values based on the analyzed data sets (see Data analysis); i.e., the number of sequences in a sample was set to n = 25, θ = 10 (θ = 4Nμ), and the level of intraspecific divergence was set to μt = 180. The ratio of nonsynonymous to synonymous mutations was set to μ_n/μ_s = 1, corresponding to assuming that ∼¾ of all new nonsynonymous mutations are immediately lethal to the organism.

Test statistics: We explored several tests of the hypothesis that the mean ages of segregating sites in two distinct classes of sites (e.g., neutral and selected) in a sample are equal. The basics of these tests are the same: a score is calculated for each segregating site. Denote the value of this score for the ith site by ω_i. Then, we define the test statistic Ω=∑i(1−Ii)⋅∑iIiωi∑iIi⋅∑i(1−Ii)ωi, (2) where I_i is an indicator function that returns 0 if site i is a neutral site and 1 if site i is selected. Thus Ω is the ratio of the average of ω for selected sites to the average of ω for neutral sites. The null distribution of Ω under the hypothesis that the value of ω_i is independent of the class of site i can be approximated by repeatedly recalculating Ω after random relabeling of sites such that the total number of selected sites remains constant and the total number of neutral sites remains constant. A total of 100,000 replicates of this resampling scheme were used to establish critical values of Ω. If we denote the ith value of Ω by Ω_i and the observed value of Ω by Ω_obs, then a test based on k resampled samples rejects at the 5% level if (1/k)R^k_i=1I_{(Ωobs≥Ωi)} < 0.05, where I₍₎ is an indicator function that returns 1 if the Boolean condition is true and 0 otherwise. We refer to the statistic Ω as the age-of-mutation (AOM) statistic.

The power of such a test to detect differences in mean age between classes is obviously dependent on how strongly ω_i is correlated with the age of the mutation responsible for polymorphic site i. The power of the test using the different estimators for ω was evaluated as described in Computer simulations. One obvious candidate for ω is an estimator of the age of a mutation under a neutral equilibrium infinite sites model suggested by Griffiths and Tavaré (1999). This estimator uses E(T | X) (where T is the age of a mutation and X is the observed data) as an estimator of T. E(T | X) is calculated by integrating over all possible gene genealogies (branch lengths and topologies) while weighting each genealogy by its likelihood. The integral is evaluated by the Monte Carlo method described in Griffiths and Tavaré (1995, 1999). This estimator may have desirable statistical properties because all the information in the data regarding T is applied in the estimation of T. However, application of this estimator is computationally very time consuming. Furthermore, in cases where the data do not conform to the infinite sites model, it is necessary to remove conflicting sites or conflicting sequences. For this reason, considerable effort was devoted to developing additional estimators that could be calculated fast and did not require elimination of any data.

Given only b, the observed number of chromosomes bearing a mutation in a sample of n chromosomes, the expected age (a) of the mutation under strict neutrality is given by E(a∣b)=2(n−1b)−1∑j=2n(n−jb−1)n−j+1n(j−1) (3) (Griffiths and Tavaré 1998). Results such as this, based on the marginal frequency of a mutation in a single site, can obviously be applied to estimate ages of individual mutations. However, they were in general found to perform much worse than estimators employing more of the data. Also, estimators based on neutrality such as E(T | X) may not perform well when the assumption of selective neutrality is violated. Instead, a heuristic estimator was devised on the basis of estimation of the structure of the underlying genealogy. In this method, an estimate of the underlying genealogy was obtained by the unweighted pair group method using arithmetic averages (UPGMA; Sokal and Michener 1958). Mutations were subsequently assigned to branches (edges) using Fitch and Wagner parsimony (e.g., Fitch 1971), and the expected age of a mutation was estimated as the average age of the branch on which the mutation occurred, i.e., the average of the time of the start of the branch and the time of the end of the branch. In cases where two or more mutations were needed to account for the observed site pattern, the average of the estimates obtained for the different assignments of the mutation was used as an estimator. We denote this estimator by â. Note that under the infinite sites model, â is the expected age of a mutation conditional on the estimated genealogy.

This estimator (â) is applied in addition to E(T | X) because it can be calculated quickly, making it suitable for simulation purposes, and because it does not require the elimination of data when the assumptions of the infinite sites model are violated. However, we do not suggest that this heuristically derived estimator should be used for more general purposes. It was chosen here solely because it provided more power in the test than estimators based on the marginal site pattern such as estimators derived from Equation 3.

Power analysis: The power of the test statistic given in Equation 2 and of the McDonald/Kreitman test to detect the presence of deleterious mutations was evaluated by computer simulations as described in Computer simulations for different values of the selection coefficient. The proportion of simulated data sets that gave significant test statistic values was tabulated for a range of values of Ns for θ = 10 and for a range of values of θ for Ns = –3.

McDonald/Kreitman tests (McDonald and Kreitman 1991) of simulated data were performed as a G-test in a 2 × 2 contingency table using the Williams correction for small sample size (see, for example, Sokal and Rohlf 1981). The numbers entering the contingency table were the number of fixed selected and neutral mutations and the number of segregating selected and neutral mutations in the sample. Empty cell entries were ignored when calculating the test statistic. This mimics a common method for performing the McDonald/Kreitman test and should not inflate the chosen significance level markedly under biologically reasonable parameter values.

Data analysis: The null hypothesis that we wish to test is that the average age of nonsynonymous mutations is the same as the average age of synonymous mutations. To analyze the data, all polymorphisms were partitioned into nonsynonymous and synonymous polymorphisms. To do this, the following heuristic approach was applied: the most common sequence was chosen as a reference sequence representing a hypothesized ancestral state. This was done only to establish an unambiguous criterion to determine if a mutation was synonymous or nonsynonymous. The identity of a variable site in a codon was determined by comparison to the homologous codon in the reference sequence. In this manner, all polymorphisms could be scored unambiguously as nonsynonymous or synonymous even when there are multiple segregating sites in a single codon. If three nucleotides were segregating in the same site, two mutations were inferred and represented as two variable sites. In most cases only one mutation was observed in each codon, suggesting that the effect of this heuristic approach does not seriously bias the analysis. Also, because the approach was applied blindly with regard to the class of the mutation, it should not inflate the chosen significance level in the resampling test.

For each data set, two resampling tests were performed, one using E(T | X) and one using the heuristically derived â as estimators of the age of a mutation. In each case, the data were resampled 100,000 times to provide a P value. When E(T | X) was applied, the Watterson (1975) estimator of θ was used to drive the simulations, and 100,000 genealogies (runs of the Markov chain) were sampled to provide estimates of the age of each mutation (see Griffiths and Tavaré 1995). Calculation of E(T | X) requires that the data strictly conform to the infinite sites model. For example, Harding et al. (1997) chose to eliminate apparent recombinant sequences. In the present analysis we instead eliminated sites to obtain a subset of sites that are consistent with the infinite sites model. The following heuristic algorithm was implemented for this purpose. Two diallelic sites are in conflict if all of the patterns in Table 1 are observed, where the labeling of state (0, 1) is arbitrary. Let C_i be the number of sites that conflict with site i. Then site i is eliminated from the data set when C_i = max{C₁, C₂,..., C_h}, where h is the number of sites. This elimination scheme is repeated until C₁ = C₂ =... = C_h = 0. The test based on E(T | X) was performed on these reduced data sets, whereas the test based on â was performed on the full data sets. Data sets so small that a significant result was not theoretically possible were not included in the analysis. Four data sets were eliminated for this reason.

Amino acid saturation dynamics: To examine the relative degree of mutational saturation in nonsynonymous and synonymous substitutions in mtDNA, the number of nonsynonymous and synonymous nucleotide differences between pairs of species was plotted for 13 mtDNA genes (NADH1, NADH2, NADH3, NADH4, NADH4L, NADH5, NADH6, ATPase 6, ATPase 8, COI, COII, COIII, Cyt. b). The species included in this analysis were Homo sapiens (Andersonet al. 1981); Pan paniscus, Pongo pygmaeus (Horaiet al. 1995); Mus musculus (Bibbet al. 1981); Gallus gallus (Desjardins and Morais 1990); and Drosophila melanogaster (Garesse 1988). Sequences were aligned by eye and nonsynonymous and synonymous nucleotide differences were counted as described in Data analysis.