Frequentist Estimation of Coalescence Times From Nucleotide Sequence Data Using a Tree-Based Partition
Hua Tang, David O. Siegmund, Peidong Shen, Peter J. Oefner, Marcus W. Feldman

Abstract

This article proposes a method of estimating the time to the most recent common ancestor (TMRCA) of a sample of DNA sequences. The method is based on the molecular clock hypothesis, but avoids assumptions about population structure. Simulations show that in a wide range of situations, the point estimate has small bias and the confidence interval has at least the nominal coverage probability. We discuss conditions that can lead to biased estimates. Performance of this estimator is compared with existing methods based on the coalescence theory. The method is applied to sequences of Y chromosomes and mtDNAs to estimate the coalescent times of human male and female populations.

FOR a sample of DNA sequences that exhibits variation, a quantity of interest is the time until all of the sequences coalesce. This quantity is termed time to the most recent common ancestor (TMRCA). The age of the common ancestor indicates the degree to which the sample sequences relate to one another. Strictly speaking, we aim to estimate the TMRCA of the sample, which may not coincide with that of the underlying population. This issue is discussed further in the next section. Because unlinked loci have different genealogies, this ancestor is defined for a specific genetic locus. In studies of human populations, two loci have received great attention because of their implications for the demographic histories of the male and female human populations: the nonrecombining region of the Y chromosome and mitochondrial DNA (mtDNA), respectively. In addition, we may imagine a situation in which the sample consists of sequences sharing some unique event polymorphism (UEP). In this case, the TMRCA represents a lower bound for the age of the mutation that defines this polymorphism. Other applications of TMRCA are illustrated in Takahata et al. (2001) and Ruvolo (1996).

It is, therefore, desirable to have an estimator of TMRCA whose statistical properties are well understood. In the next section, we review some of the methods that have been developed recently. Most of these methods make assumptions about the population structure, e.g., population size and rate of growth, amount of immigration, etc., as well as the mutation process. The robustness of these methods against model misspecification has not been characterized. A recent article by Stumpf and Goldstein (2001) advocates a model-free approach in population studies. Unfortunately, their “model-free” approach actually assumes a very special model, namely a star-shaped genealogy. Our goal in this article is to find an estimator of TMRCA that does not assume that the population parameters are known and whose statistical properties can be assessed easily. We propose an inferential approach to TMRCA that explicitly models the mutation process but not the population history. Simulations demonstrate that the estimate is relatively unaffected by some important components of population structure.

METHOD

Interpretations of the estimator: Before reviewing the existing methods and introducing the new one, we wish to clarify the genetic and statistical interpretations of the quantity under investigation.

We restrict attention to the estimation of the sample TMRCA. Whether the TMRCA of a sample can be interpreted as an approximation of the TMRCA of a population depends on the history of the population as well as the sampling strategy. In the simplest case, for a random sample of size n from a panmictic population, the probability that the common ancestor of a sample coincides with the common ancestor of the population is ∼(n - 1)/(n + 1) (Saunderset al. 1984). In practice, the panmictic assumption is seldom satisfied, and the sampling scheme is rarely completely random, which lead to uncertainties in the interpretation of the sample TMRCA (Templetonet al. 2000).

Statistically speaking, because of the stochastic nature of the evolutionary process, TMRCA can be regarded as the (unobserved) value of a random variable. More precisely, the DNA at each locus evolves as one realization of a stochastic process. Further, the TMRCA of one sample may differ from that of another sample. However, for a fixed sample and at each locus, there is an unknown true TMRCA, denoted by T, at which all individuals in the sample share a common ancestor. We may imagine that an immortal person has kept a genealogy of each locus since the beginning of human history. If we had the genealogy, we would know T. Since we do not have the genealogy, we estimate T using the genetic variation of the sample at that locus. In statistical terms, we estimate T conditional on its true value. Thus, in spite of our view that T is a random variable, in our analysis it becomes an unknown parameter that is estimated on the basis of a sample of current DNA.

Review of existing methods: The foundation of the majority of existing methods of estimating TMRCA is Kingman’s (1982) coalescence theory. It determines a probability distribution of TMRCA solely on the basis of the assumed population model, i.e., in the absence of any sequence data. For example, under the standard model of a panmictic population of constant size N (N → ∞), coalescence theory stipulates that the waiting times between consecutive coalescent events are independent exponential random variables with known parameters. When there is also mutation, an important parameter is θ= 4Nμ1 (or θ= 2Nμ1 for a haploid population), where μ1 is the mutation rate per genetic locus. When time is measured in units of N generations, θ is simply twice the mutation rate per unit time. Most existing methods share a common structure with three components:

  1. Choose a population model and parameters, which may include population size, growth rate, and migration rate.

  2. Formulate a prior distribution of TMRCA according to the coalescent.

  3. Update the prior with a likelihood function or summary statistics derived from the data.

In practice, the selection of a population model involves somewhat arbitrary choices because of our ignorance of demographic history; it is also constrained by the computational tools currently available. Thus, these methods vary in their choice of population parameters and the method for evaluating the likelihood function.

One group of methods places prior distributions on one or more of the model parameters and reduces sequence data to one or a few summary statistics, such as the number of segregating sites, heterozygosity, or maximum pairwise number of nucleotide differences. With values of the parameters chosen from their prior distributions, proposal genealogies are simulated using the coalescent. The likelihood is the probability of observing data with identical (or similar) summary statistics to those observed in the sample. These methods often implement rejection sampling, where each proposed genealogy is accepted with a probability proportional to the likelihood. The TMRCA distribution of the accepted genealogies is taken as the posterior distribution. If the model parameters are chosen from prior distributions, the same procedure also produces posterior distributions of population parameters, such as population size, growth rate, etc. Examples of these methods can be found in Tavaré et al. (1997) and Pritchard et al. (1999). We call these methods Bayesian approaches because they view T as a random variable with the prior distribution determined by the coalescent and summarize conclusions in terms of the posterior distribution of T given the data (or summary statistics of the data).

In contrast, a second group of methods estimates the model parameters on the basis of the likelihood of the complete data. The posterior distribution of TMRCA is computed at the estimated parameter values without regard to the variability due to estimation of unknown parameters. Evaluation of the full likelihood can be accomplished by methods such as importance sampling or Markov chain Monte Carlo. We refer to these methods as “empirical Bayes approaches,” to indicate that the posterior distributions of the TMRCA are computed at the maximum-likelihood estimate of population parameters. Examples of empirical Bayes methods include algorithms implemented in genetree (Griffiths and Tavaré 1994; available from http://www.stats.ox.ac.uk/mathgen/software.htm) and the larmarc package (available from http://evolution.genetics.washington.edu/lamarc.html), which includes programs fluctuate (Kuhneret al. 1998), recombine (Kuhneret al. 2000), and migrate (Beerli and Felsenstein 1999). These programs have the advantage of producing maximum-likelihood estimates of model parameters as well as the posterior distribution of the sample TMRCA. At present, however, the convergence behavior of the Markov chain Monte Carlo (MCMC) algorithms is still not completely understood in the context of ancestral processes. We find that these MCMC procedures require considerable experience in the “art” of tuning the MCMC parameters and patience for the substantial computational costs incurred.

An estimator that is similar in certain respects to ours was suggested by Templeton (1993), who exploited expressions for the conditional moments of the coalescent time T2 between two randomly sampled sequences given k, the observed number of nucleotide differences between those sequences. The expressions, derived by Tajima (1983), are E(T2k)=θ(1+k)2μ1(1+θ) (1) and var(T2k)=θ2(1+k)4μ12(1+θ)2, (2) where θ= 2Nμ1 (haploid model as above). For a sample of n sequences, Templeton evaluated D, the number of pairwise nucleotide differences averaged over all pairs whose most recent common ancestor is the root of a reconstructed genealogy. Substituting D for k, he then applied (1) and (2). This method has been criticized on two grounds (Tavaréet al. 1997). First, because the derivation of (1) and (2) uses coalescence theory, Tajima’s results are applicable only for two randomly sampled sequences, while Templeton’s method uses pairs selected to be more different than average. Second, since D is derived from many dependent pairs of sequences, it has different statistical properties from k, the number of nucleotide differences between two observed sequences.

It is important to recognize that all of these estimators of TMRCA assume knowledge in the form of prior distributions or precise values of population parameters such as effective population size, growth rates, or migration rates, which are difficult to estimate accurately. Taking into consideration the uncertainty of the estimated parameter values will increase the sampling variance of the TMRCA estimators. This is demonstrated by the Y chromosome microsatellite analysis of Pritchard et al. (1999), which concludes that the TMRCA for the worldwide Y population may range between 16,000 and 126,000 years.

Worse yet, most populations have gone through periods of expansion and bottleneck. The existing parametric models are quite restrictive and may fail to describe the population history adequately. In particular, genetree, fluctuate, migrate, and micsat employ only a few greatly simplified models. For example, both genetree and fluctuate assume a constant and deterministic rate of exponential growth; migrate assumes constant population size and constant migration rate with respect to time. A large number of simulations would be required to investigate the bias and variance associated with those estimators of TMRCA when the underlying population model deviates from those allowed in these programs. We are not aware of any such study.

We now describe a new method that is explicit about the mutation model but not the population model.

Model and assumptions: We make the following assumptions on the mutation process:

  1. Mutations occur as a Poisson process with constant rate along all branches (i.e., constant molecular clock).

  2. The mutation rate is the same known constant at all sites.

  3. Mutation is independent among all sites.

  4. The mutation rate is sufficiently low that the infinitely many sites model is a reasonable approximation (i.e., there are no recurrent mutations).

    In addition, we assume the following:

  5. The sample consists of fully sequenced regions of DNA. There is no recombination.

  6. Generations are nonoverlapping and of constant length.

Most of these assumptions are made in other methods as well. As explained in a later section, assumption 4 is not necessary, but is made to simplify the exposition. As long as 4 holds, 2 is not strictly required, provided the total mutation rate for the locus is known. Since we take this rate from phylogenetic estimates, where it is given as a per site rate, and since we want to apply our methods to mtDNA data, where hypothesis 4 is untenable, we are in practice required to make the assumption 2 in some applications.

Point estimate of TMRCA: Our method is motivated by the following observation. As a consequence of the constant molecular clock hypothesis and the infinitely many sites mutation model, the number of nucleotide differences dij between any two sequences i and j (not necessarily randomly sampled) follows a Poisson distribution with mean 2μ2113;T2, where μ is the mutation rate per site, 2113; is the length of the sequenced region, and T2 = T2(i, j) is the coalescent time of the sequences from their common ancestor measured in the same time units as the mutation rate. In particular, the expected number of nucleotide differences between pairs of sequences in a sample that diverged from a common ancestor is proportional to the TMRCA of the sample. This suggests a frequentist estimator: T^2=dij2μ. (3) If T2 is the true coalescence time for a pair of sequences, by basic properties of the Poisson distribution, we have E(T^2T2)=T2 and var(T^2T2)=T22μ. Thus, 2 is an unbiased estimator of T2. A plug-in estimator for its variance is simply dij/(4μ22113;2).

To generalize 2 to a sample of n sequences, we must solve two problems. First, not all pairs of sequences in the sample share the same TMRCA. Second, because the sequences are related through a genealogy, distances between pairs of sequences have complicated joint distributions.

To solve the first problem, we follow Templeton (1993) and partition the sample of sequences into two subsets, corresponding to the two clades that appear to have diverged from the common ancestor of the sample. Denote the two clades as L and R ; clade L includes sequences l1, l2,..., lp, and clade R includes sequences r1, r2,..., rq, where p + q = n. Next, we compute the average pairwise nucleotide differences between L and R . Let D=1pqiLjRdij. (4) The estimator for the TMRCA we propose is T^=D2μ. (5) Keeping in mind that μ1 in (1) and (2) equals μ2113; here, we find it instructive to rewrite Templeton’s estimator as T^Temp=θ1+θ(12μ+D2μ)=θ1+θ(12μ+T^), which illustrates that Templeton’s estimator is a linear transformation of our estimator. The shrinkage by a factor of θ2113;/(1 +θ2113;) and the offset of 1/(2μ2113;) arise from Templeton’s use of the assumed constant population size model as a prior.

Partitioning algorithms: In practice, to select L and R , one can use any tree-constructing algorithm and take the deepest branch (Felsenstein 1993). The root of the tree may be established by comparison with an outgroup, for example, using chimpanzee for human data. Although the topology of the entire genealogy is seldom known, and different criteria may lead to more than one plausible tree, the two clades that diverge most from the ancestor can often be constructed with less ambiguity. Templeton observed that all of the plausible genealogies in his example agreed on the branching pattern at the root. Our simulation results suggest the same unless the true genealogy is extremely star-shaped.

To analyze the simulated data, which do not include an outgroup sequence, we take advantage of two heuristic ideas to estimate the two clades directly. First, when the underlying genealogy is not star-shaped, the pair of sequences that differ most from each other are most likely to represent the two clades; each of the remaining sequences can be “classified” into one of the clades by its relative similarity to the existing sequences in each clade. Second, in the case of a perfectly star-shaped genealogy (that is, all sequences have evolved independently from the root of the genealogy), the above algorithm may lead to a completely incorrect partition and a biased point estimate. We prefer a random partition that produces a point estimate with less bias than would result from a deterministic partition. To accomplish this, we suggest the following “two-step iterative” partition algorithm, which seems to combine the best features of both these ideas.

  1. Compute the Hamming distance matrix K for the sample, where Kij is the number of nucleotide differences between sequences i and j.

  2. Find a pair (i*, j*), such that Ki*j*Ki,j for all pairs (i, j). Assign i* and j* to the opposite clades.

  3. Add each of the remaining sequences, in random order, to the clade closer to it, as measured by the mean pairwise Hamming distance of the candidate sequence to the current members of the clades.

  4. Compute the average pairwise nucleotide difference, D0, by Equation 4.

  5. Compute the upper triangular matrix, Mij = Λ(Kij, D0) (i, j = 1, 2,..., n, ij), where Λ(x, λ) is the Poisson probability density function with mean parameter λ, evaluated at x.

  6. Compute the likelihood ratio matrix, Lij = Mij/maxk,l(Mkl) (ij). If Lij is less than a prespecified threshold (which we set at 0.2), set it to 0.

  7. Choose a pair (i, j) according to the multinomial probability proportional to Lij, and assign i and j to the opposite clades. This pair of sequences, i and j, takes the role of i* and j* in step 2.

  8. Repeat step 3 to obtain a partition.

For our simulation, we used this algorithm to obtain 20 partitions and took the averaged point estimates and variance estimates over all the partitions.

Sampling variance of T: The second problem in generalizing 2 in (3) is the correlation between pairwise distances, dij. Although the correlation does not affect the point estimate, it must be accounted for in the computation of the sampling variance of : σ2=var(T^T)=var(DT)4μ22. (6)

To estimate var(D|T) from data, we ignore the variability due to the uncertainty of the clade partition. Then var(DT)=1p2q2(iLjRvar(dijT)+2×i,kL,ikjRcov(dij,dkjT)+2×iLj,lR,jlcov(dij,dilT)+2×i,kL,ikj,lR,jlcov(dij,dklT)). (7) Since each dij is a Poisson random variable, its variance is equal to its mean. Let ik be the most recent common ancestor of sequences i and k, and let Tik be the time between i and k (twice their coalescence time). Given Tik, the mutation processes from the two sequences, i and k, to their most recent common ancestor, ik are independent; i.e., cov(dik,i,dik,kTik,T)=0 . Hence, for example, for i,kL,jR cov(dij,dkjT)=E[cov(dij,dkjTik,T)T]=E[E(ik,jTik,T)T]=E(dik,jT), so var(dijT)=E(dijT),iL,jR (8) cov(dij,dkjT)=E(dik,jT),i,kL,jR (9) cov(dij,dklT)=E(dik,jlT),i,kL,j,lR. (10) A plug-in estimator for (8) is dij. A plug-in estimator for (10) is dik,jl , and this can be approximated by counting the number of bases that are the same in each i and k and in j and l but are different between, say, i and j. Likewise, a plug-in estimator for (9) is the number of sites that are the same in sequences i and k, but are different in sequence j. Thus, we can estimate var(D|T) by var^(DT)=1p2q2(iLjRdij+2×ik,i,kLjRdik,j+2×iLjl,j,lRdi,jl+2×ik,i,kLjl,j,lRdik,jl). (11)

Substituting (11) into (6) gives the sampling variance of . The amount of computation required for the variance calculation is on the order of p2q2, which should not pose any difficulty unless the sample is very large and the tree is balanced (i.e., p = q = n/2). Further, the variance can be calculated sequentially. We observe very fast convergence of the variance estimate by adding sequences from alternating clades.

Confidence interval for T: Because is a weighted sum of Poisson random variables, its sampling distribution is skewed to the right. Our simulations indicate that because the distribution of is skewed, confidence intervals based on a normal approximation can have incorrect coverage probability in samples with a small number of segregating sites. To correct for the skewness, we utilize a square root transformation. We propose to approximate the sampling distribution of T^ by the normal distribution; i.e., T^N(T,σ24T), where σ2 is defined in (6) and estimated in (11). We observe that the MSE of can be reduced slightly by an adjustment based on the square root transformation. The modified estimator, b, is T^b=T^σ^24T^, (12) where σ^2 is the sampling variance estimated by (11). Tables 1, 2, 3, 4 show that the confidence intervals constructed this way have coverage probabilities close to or higher than the nominal levels, and the bias is often negligible.

Correction for recurrent mutation: The exposition so far has assumed the infinitely many sites mutation model. Our simulation program generates sequences according to a model involving two states per site. The pairwise distance, dij, is the number of nucleotide differences between sequences i and j. To account for recurrent mutations, we replace all dij with a corrected distance measure dij , where dij=2log(12dij) (13) (Nei 1987). If a sample contains sites with three and more alleles, Kimura’s two parameter model can be used to correct for recurrent mutations (Nei 1987).

Figure 1.

—Scatter plot of true TMRCA vs. Tˆ (dots) under a constant, panmictic population model, with 95% confidence intervals calculated using (6) and (11) (vertical bars). Of 1000 simulated genealogies, 953 true TMRCAs are covered by the 95% confidence interval.

Computational details: For the numerical studies reported in the next section, sequence data are simulated on the basis of coalescence. The estimation procedure described above is implemented in MATLAB and is applied to both the simulated data and the real Y chromosome and mtDNA polymorphism. For genetree estimates of TMRCA in the simulations, we input the assumed values of θ and use the mean TMRCA of 500,000 sample genealogies. The publicly available version of the program fluctuate does not estimate the sample TMRCA. However, by modifying the source code, we record the TMRCA of all the sample genealogies in the last long chain and use the mean of this distribution as an estimate. The parameters for the MCMC sampling used to estimate TMRCA for the Y and mtDNA sample are 10 short chains of 100,000 steps each and 5 long chains with 106 steps each.

RESULTS FROM SIMULATED DATA

Simulations of a panmictic population: Figures 1 and 2 display results from two simulation experiments. In the first experiment, 1000 genealogies from a sample of 20 haploid individuals are simulated under a model of a panmictic population whose effective size has remained constant at N = 5000. Mutation events are superimposed on each genealogy. The mutation model is quite simple: In a region of 15 kb, each site mutates between two states at a rate of μ= 2.5 × 10-7/generation. These figures are comparable to what is known about mtDNA. We adopted a two-state mutation model rather than a four-state model for computational simplicity. Recurrent mutations do occur in these simulations and are corrected by (13).

Figure 2.

—Scatter plot of true TMRCA vs. Tˆ (dots) under an exponential growth, panmictic population model, with 95% confidence intervals calculated using (6) and (11) (vertical bars). Of 1000 simulated genealogies, 960 true TMRCAs are covered by the 95% confidence interval.

The second experiment is done in a similar fashion: Each genealogy of 20 sequences is simulated assuming that the panmictic population has grown exponentially at a rate of g = 1 × 10-3/generation with a present effective population size of N0 = 4 × 104. The mutation rate and sequence length correspond to what is known about the Y chromosome: 2113; = 50 kb and μ= 2.5 × 10-8/generation (Thomsonet al. 2000). The vertical bars in Figures 1 and 2 are the 95% confidence intervals.

In the same fashion, we have simulated genealogies under a panmictic population model with different values of the mutation rate and sequence length (Table 1), effective population size (Table 2), growth rate (Table 3), and the sample size (Table 4). Under each parameter set, we evaluate the performance of the estimator with two summary statistics: bias and coverage probability. For a measure of (relative) bias, we compute Bias=1Ii=1IT^ibTiTi, (14) where I is the number of genealogies simulated in an experiment. Coverage probability is defined as the proportion of runs in which the confidence interval covers the true T. Each of Tables 1, 2, 3, 4 varies one aspect of the mutation model, the population model, or the number of samples, while keeping all other aspects fixed.

View this table:
TABLE 1

Coverage probability of the confidence intervals and average bias of the point estimates as functions of mutation rate and sequence length

View this table:
TABLE 2

Coverage probability of the confidence intervals and average bias of the point estimates as functions of the present effective population size, N0

Several trends become evident from the simulation results. First, the bias of the estimate becomes nonnegligible in two situations. In one, the mutation rate is low or the sequenced region is short (Table 1). In the other, the present population has resulted from a small founding population that has gone through rapid growth (Table 3). In this case, the genealogy tends to be short and star-shaped, and nearly all branches evolve independently. In both situations, the resulting DNA samples contain only a few segregating sites for most of which there is a single representation of the rare type. Thus, the observed genetic variation contains little information regarding the relationship among the sequences. Consequently, distance-based partition algorithms often introduce bias. The two-step algorithm described above reduces the bias by averaging over 20 partitions. Parsimony- or likelihood-based algorithms are less likely to introduce bias in these situations. Second, as seen in Table 4, a larger sample size will reduce the mean square error (MSE) and the width of the confidence interval, but the improvement tends to be negligible, especially when the sample size is already moderately large. Finally, the total heterozygosity (summed over all sites), H, increases with the mutation rate, the length of the sequenced region (Table 1), and the population size near the root of the genealogy (Table 2); it decreases as the growth rate increases (Table 3), but is roughly constant with changing sample size (Table 4). Hence, H can be used to predict the performance of the estimator on a particular data set. Roughly speaking, an H value <5 provides a warning sign that the data do not contain enough information, in which case special care should be taken in the partition step by using non-distance-based tree-building algorithms or by considering more than one partition.

View this table:
TABLE 3

Coverage probability of the confidence intervals and average bias of the point estimates as functions of the population growth rate, g

Next, we compare our estimator, b, to two existing methods: the Bayesian method based on summary statistics and the empirical Bayes approach as implemented by genetree, using complete sequence data.

View this table:
TABLE 4

Coverage probability of the confidence intervals and average bias of the point estimates as functions of number of sequences included in the sample

Comparison with methods based on summary statistics: As an example, we compare our approach with the rejection method of Tavaré et al. (1997). A data set is simulated under a chosen population model, and its TMRCA is estimated by (5). For the rejection method, we assume perfect knowledge of both the population model and parameter values. Each proposed genealogy is accepted with Prob(k segregating sites|the length of the proposal genealogy). Figure 3 shows that the 95% confidence interval constructed by our frequentist approach is narrower than the 95% credibility interval from the rejection method. This is hardly surprising, since our frequentist estimator uses the complete sequence information, while the rejection method is based on summary statistics.

Comparison with methods based on the complete sequence data: It is important to compare the relative efficiencies of our estimator and a coalescent sampling estimator that uses complete data instead of summary statistics. We compare the MSE of the two approaches for the model of constant population size. Recall that the mutation rate parameter for the entire sequenced region is θ= 2Nμ2113;.

First, consider an idealized situation, in which we not only observe DNA sequences, but we also know the positions of each mutation relative to the coalescent events. A schematic genealogy is shown in Figure 4. The true times between two consecutive coalescence events are Wj, indexed by the number of lineages existing during that time interval. We do not observe Wj, but we do know the number of mutations that occurred during time interval Wj. By coalescence theory under the constant population model, estimation of TMRCA amounts to estimating each Wj (j = 2, 3,..., n - 1) independently, on the basis of Sj, the number of mutations occurring in the corresponding time interval. Consider a Bayesian estimator and assume the population size, N, is known without error. The prior distribution of Wj, based on coalescent theory (Tavaréet al. 1997), is Prior(Wj)exp(j(j1)2). Modeling the mutation as a Poisson process, we write the likelihood of observing Sj = sj mutations during time Wj as P(Sj=sjWj=wj)=e(12)θjwj((12)θjwj)sjsj!. The posterior distribution of Wj can be derived by Bayes theorem, P(WjSj)=Prior(Wj)P(SjWj)wPrior(w)P(Sjw)dw, and is a gamma distribution: L(WjSj)=Γ(1+Sj,2j(j1+θ)). Let j be a Bayes estimator taken as the mean of the posterior distribution of Wj. We have W˘j=2(1+Sj)j(j1+θ). The MSE associated with j is MSE(W˘j)=ws=0(w˘w)2P(sw)×Prior(w)dw=4j2(j1+θ)(j1). (15)

Figure 3.

—Comparison of our approach to the posterior distribution obtained by rejection methods. The histogram is the posterior distribution obtained by the rejection method. The curve is the probability density function of the square of a Formula variate. The data are simulated under a constant size model, N = 5000, g = 0, μ= 5 × 10-8, 2113; = 50 kb.

Figure 4.

—A schematic genealogy with four samples.

Now we consider a frequentist estimator, j, which is based only on the number of mutations in each time interval of duration Wj and is defined by Wj=Sjθ2j. The MSE associated with this estimator is MSE(Wj)=ws=0(ww)2P(sw)×Prior(w)dw=4j2θ(j1). (16) For the entire genealogy with n sequences, let T˘=Σj=2nW˘j and T=Σj=2nWj . Then MSE(T˘)=j=2nMSE(W˘j),MSE(T)=j=2nMSE(Wj) since all cross terms are 0.

Comparing (16) with (15), we see that MSE(j) < MSE(j) for each j and all values of θ. This is expected, as a Bayes estimator minimizes MSE under the assumed model. The difference between the two MSE terms, however, is Δ(j)=4(j1+θ)θj20for largeθ. (17) The computation so far establishes lower bounds on the MSE of a Bayesian and a frequentist estimator, respectively. Both MSE terms converge to a finite limit as n → ∞, while the difference converges to 0 as θ → ∞.

Our next task is to examine how close to the lower bound one gets by using the existing Bayesian or our frequentist estimator. We simulated genealogies and DNA sequences under a specific population model. For the Bayesian estimate, we used genetree with assumed values of θ, and the estimate is based on 500,000 sampled genealogies. Table 5 indicates that both the Bayesian and frequentist estimators approach their respective theoretical lower bounds for large θ. The discrepancy between the MSE of our estimator, , and the theoretical lower bound is due to our imperfect knowledge of the position of each mutation relative to the coalescence events. The discrepancy between the MSE of the gene tree estimator and the Bayesian theoretical lower bound is due to two factors: imperfect knowledge of the position of each mutation relative to the coalescence events and the sampling variance due to the finite number of trees sampled by the program. The simulation study in Stephens (2001) shows that millions of genealogies are required to reduce this latter source of variation. Finally, the value of θ= 2Nμ2113; is seldom known. In fact, a primary goal of genetree is to estimate θ empirically. This additional source of uncertainty would increase the MSE associated with TMRCA from the genetree estimator.

View this table:
TABLE 5

Formula associated with Bayesian and frequentist estimators of TMRCA

APPLICATIONS TO Y-CHROMOSOMAL DNA AND mtDNA

As examples, we apply our method to the worldwide samples of Y chromosome and mtDNA we have recently sequenced.

Y chromosome TMRCA: The Y chromosome data are from n = 108 individuals sampled worldwide. A region of 69,000 bp has been screened, and 114 single nucleotide polymorphisms (SNPs) have been found. The assumed mutation rate is derived from the divergence between chimpanzees and humans in the corresponding genomic region: μ=substitutions between chimp and human2Tdiv. Assuming a divergence time of 5 million years between chimpanzees and humans and a constant generation length of 25 years, the mutation rate is ∼3 × 10-8/(site × generation) (Thomsonet al. 2000). Under a model with constant population size and panmixia, the program fluctuate produces an estimate of 10,300 for the effective population size and 117,000 years for the worldwide TMRCA. If we assume an exponential growth model, however, the same program estimates the rate of growth to be 1.3 × 10-3, the current effective population size to be 22,400, and the worldwide TMRCA to be 75,000 years (H. Tang, R. Thomson, L. L. Cavalli-Sforza, P. Shen, P. J. Oefner and M. W. Feldman, unpublished results). Clearly, a Bayesian approach depends heavily on the assumed population model. In the application of our approach, the clades are formed on the basis of the haplogroup definitions in Underhill et al. (2000). In particular, for the worldwide TMRCA, one clade is represented by the five San and Pygmy samples in haplogroup I, while haplogroups II-X form the other clade. The sample heterozygosity is ∼7.49. TMRCA estimates for the world and for the five continents are given in Table 6.

View this table:
TABLE 6

Estimated Y chromosome TMRCA for each continent and for the world

Although the sample heterozygosity is relatively small, the simulations reported above suggest that the confidence interval is still reasonable. In view of the care with which the tree has been constructed (Underhillet al. 2000), we believe the bias due to the partition choice in the point estimators is relatively insignificant.

The coalescent time of the African Y chromosomal sequences is quite similar to that of the whole world, but much older than that of other continents, providing support for the out of Africa hypothesis. That the values for the estimated TMRCA of all non-African continents are similar does not necessarily suggest simultaneous settlement on these continents; rather, it may be a consequence of population bottlenecks. In fact, the TMRCA represents only a lower bound of the founding date and can be substantially more recent than the time of the settlement.

mtDNA TMRCA: The mtDNA data are from 179 individuals sampled worldwide. The entire mitochondrial genome, except the hypermutating D-loop, has been sequenced, and 971 SNPs have been found. The mutation rate has been derived from the divergence between chimpanzees and humans in the corresponding genomic region. Because of the relatively high mutation rate, we used Kimura’s two-parameter model to correct for recurrent mutations (Nei 1987). The average mutation rate is ∼2.43 × 10-7/(site × generation). The sample heterozygosity is ∼39. For the tree-partition step, we used dnapenny. The estimated TMRCA for individual continents and for the world are given in Table 7.

View this table:
TABLE 7

Estimated mtDNA TMRCA for each continent and for the world

As with the estimates for the Y chromosome TMRCA, the worldwide mitochondrial TMRCA is comparable to that of the African TMRCA, while mtDNA from populations on all other continents coalesce much more recently. It is also interesting to note that the TMRCA for European mtDNA populations is estimated to be smaller than that for Asian, Oceanic, and American populations, which may indicate a small female founding population in Europe, combined with a high population growth rate in that continent. More on this is in the discussion.

DISCUSSION

In the construction of the variance of , we assume that the two clades are defined without error. In our simulation, we used a rather naive partitioning algorithm to reduce computation and to facilitate automation. One might have expected that the confidence intervals based on (7) would be too narrow due to the uncertainties in the clade-defining step. On the contrary, for our simulated data, the estimate and the confidence interval appear to be relatively robust. We reason that a branch that would be placed in the wrong clade is likely to have diverged from the rest of the sequences near the root; hence the point estimate and the confidence interval are not greatly affected. One can make a bias-variance tradeoff argument: A slightly mistaken but more balanced partition produces a small bias, but reduces the variability of the estimate. We do not pursue an optimal partitioning algorithm in this study.

Effects of population subdivision: So far, our simulation has dealt with panmictic populations only. A natural question is how our method performs when the sampled population is subdivided. We argue that, unless the generation length is subpopulation dependent, the constant molecular clock assumption is independent of population subdivision; therefore, this frequentist method should be valid in the presence of population subdivision. Our simulations in one-population (panmictic) situations indicate that the performance of our estimator depends on the shape of the underlying genealogy and the number of segregating sites, but does not depend on the specific mechanisms that generate such genealogies. In fact, since subdivided populations tend to generate genealogies with very long branches near the root, the sample partition tends to be more accurate. Hence, we hypothesize that the estimates of TMRCA should have smaller biases in these situations. We performed a few simulations under a model that incorporates population subdivision and migration. The goal of these simulations was not to prove the validity of our method in a subdivided population; rather, it is to provide some evidence supporting our intuitions.

A population model with subdivision and migration: We performed simulations in which sample genealogies are generated as follows: An ancient population of size N0(T1 + T2) began to grow T1 + T2 generations ago. At time T2 generations ago, this population, now of size N0(T2), was split into two subpopulations, of size N1(T2) = f1N0(T2) and N2(T2) = f2N0(T2). The two subpopulations each underwent exponential growth, at rates g1 and g2, respectively. At the end of each generation, a fraction m1 of subpopulation 1 migrates to subpopulation 2, while a fraction m2 of subpopulation 2 migrates to subpopulation 1. The ancestral population remains at the constant size, N0(T1 + T2). Figure 5 plots the true TMRCA vs. the estimated TMRCA. Specifically, T1 = 500, T2 = 2000, m1 = 0.0005, m2 = 0.0001, g0 = 0.005, g1 = 0.001, g2 = 0.003, f1 = 0.75, f2 = 0.25, N0(T1 + T2) = 300, L = 15,000 bp, and μ= 2.5 × 10-7/site × generation. The mean simulated TMRCA is 2505 generations; on average each simulated data set contains 120 segregating sites. The bias defined in (14) over 1000 genealogies is 2% and the coverage probability is 96.5%. Thus, the method seems to perform well even in the presence of population structure.

Figure 5.

—Scatter plot of true TMRCA vs. Tˆ (dots) under a subdivided population model, with 95% confidence intervals calculated using (6) and (11) (vertical bars). Of 1000 simulated genealogies, 965 true TMRCAs are covered by the 95% confidence interval.

The robustness of the estimator does break down in some situations. For a diagnostic test, we suggest using the sample heterozygosity. The coverage of the confidence interval is reasonably accurate when the sample heterozygosity is at least 15 and exceeds the nominal level otherwise. The relative bias of the point estimator is also reasonably small and becomes very small when the sample heterozygosity exceeds 5. Smaller values of H indicate a lack of information in the data. In these cases the point estimator of the TMRCA is likely to be positively biased. To reduce bias, we recommend using a true tree-building algorithm instead of our naive algorithm to partition the sample into left and right clades. At least for human populations, since rapid population growth is likely to be a recent phenomenon and since worldwide populations have experienced repeated isolation and admixture, we believe that the genealogy underlying a real sample is unlikely to have the star shape, where our method shows its poorest behavior, unless the sampled individuals are closely related.

Limitations of the frequentist estimator: As one may suspect, the chief weakness of this frequentist estimator is its sensitivity to misspecification of the mutation model. Unless the rate of recurrent mutation is very high, our frequentist estimate of TMRCA varies linearly with the estimate of the mutation rate. The uniform independent mutation model assumed here is clearly not realistic. It is possible to generalize our estimator to allow for different transition-transversion rates or to model a higher mutation rate at third codon positions. In estimating the divergence time of two species, even the validity of the constant molecular clock hypothesis can be challenged. The real problem is, however, that our current knowledge of the mutation process is limited. One could argue that, in cases where researchers have more faith in their knowledge of population histories than they do in the estimates of mutation rate, Bayesian estimators can be more robust. Nonetheless, we believe that, with the availability of comparative genomes, our understanding of the mutation process will advance at a faster rate than our understanding of the demographic histories of the studied populations.

The simulation results in Tables 1 and 4 allow us to consider a question related to the sampling scheme. To acquire more information regarding a population TMRCA, we may elect to sequence a longer region for each individual already in the sample, or we may fix the size of the sequenced region but sample more individuals. Which scheme is more efficient? The answer seems to be the former. Equations 15 and 16 indicate that the MSE decreases as θ increases, and θ is proportional to the sequence length. Coalescent theory shows that if our sample is already of reasonable size, additional individuals are likely to be closely related to individuals already sequenced and thus would add little information regarding the TMRCA. These observations are also supported by the simulation results shown in Tables 1 and 4. However, in reality our methods estimate the sample TMRCA, not the population TMRCA [although in panmictic populations these are the same with probability close to one (Saunderset al. 1984)]. A small convenience sample may not fully represent the population diversity. When a population is subdivided, to increase the confidence that the sample TMRCA is close to the population TMRCA, we should sample additional sequences from underrepresented subpopulations.

Implications of the estimated TMRCA for human populations: Comparison between the estimates in Table 6 and those in Table 7 reveals that the Y chromosome TMRCA, of the world and on each continent, are much younger than their mtDNA counterparts. H. Tang, R. Thomson, L. L. Cavalli-Sforza, P. Shen, P. J. Oefner and M. W. Feldman (unpublished results) provide an explanation that is based on the disparity between the male and female effective population sizes. The estimates of TMRCA for the Y chromosome and mtDNA populations are puzzling in another way: How could the worldwide Y chromosome TMRCA be younger than that of the mtDNA of all continents: Africa, Asia, Oceanic, and America? At first sight, this seems to imply that the most recent mtDNA common ancestors for all continents (except for Europe) resided in Africa. One could certainly claim that all the confidence intervals overlap, and therefore the worldwide Y chromosome TMRCA could still be older than the non-African mtDNA TMRCA. Alternatively, as explained before, these estimates are sensitive to the errors in the estimated mutation rate. A small error in mutation rate will amplify to a much larger discrepancy in the estimated TMRCA. Thus, a slightly underestimated mtDNA mutation rate can account for the observation. However, we favor a simple explanation based on the unequal generation length between the human male and female populations. Because of the physiological constraints and social customs, it has been observed that the average generation length can be longer for the male than the female population (Cavalli-Sforzaet al. 1964; Tremblay and Vezina 2000). In our analyses, we used a generation length of 25 years for both the male and the female populations. Our paradox can be reconciled if we assume a 25-year generation length for the females, but 30 years per generation for the males. As seen in Table 6, the adjusted worldwide Y population TMRCA then becomes 104,000 years; similar adjustment results in estimates of 40,000-∼60,000 years for the non-African continents.

Finally, it is tempting to associate the TMRCA of a sample collected on a single continent with the date of settlement of that continent. This can be misleading for two reasons. On the one hand, many of the founding human populations are likely to have been small and to have gone through expansion; hence the population TMRCA of a continent can be much more recent than the date of settlement. Second, the MRCA of the present population may be much more ancient than the founding population due to subsequent migration events. For example, in the analysis of the European Y chromosome samples, we detected two sequences that represent haplotypes that are found almost exclusively in Africa. These were not included in our analysis of the European sequences. The mtDNA of both of these individuals resemble those of other European mtDNA sequences. We suspect that there have been migrations from Africa in the male lineages of these individuals. Because of the lack of recombination on the Y chromosome, the whole African Y chromosome has been preserved. Inclusion of these two Y chromosome sequences would have increased the European Y chromosome TMRCA considerably. We detected these two sequences using the jack-knife method, repeating the estimation procedure leaving one sequence out each time. However, when the sample is small, such as in the case of the American samples, it is often impossible to detect “outliers.”

CONCLUSION

We have introduced a new method of estimating TMRCA on the basis of a sample of DNA sequences. Unlike many existing methods, this estimator does not require assumptions on the demographic model and associated parameters. Construction of point estimates and confidence intervals is simple and fast. Simulation studies have been performed under simple population models. The MSE calculation and (17) suggest that, with a sufficient amount of data, our estimator, which requires no knowledge of the population model, achieves similar predictive power to a Bayesian approach that uses exact knowledge of the demographic parameters.

Acknowledgments

We gratefully acknowledge Dr. Mary Kuhner’s advice on obtaining TMRCA estimates with fluctuate. We thank Luca Cavalli-Sforza, Neil Risch, Noah Rosenberg, and Russell Thomson for helpful discussions. Two anonymous reviewers provided useful comments and criticisms. Hua Tang is supported by a Howard Hughes predoctoral fellowship. Research was supported in part by the National Institutes of Health grants GM28016 and GM28428 to M.W.F. and by the National Science Foundation grant DMS-0072523 to D.O.S.

Footnotes

  • Communicating editor: J. B. Walsh

  • Received July 19, 2001.
  • Accepted February 26, 2002.

LITERATURE CITED

View Abstract