Estimating Divergence Time with the Use of Microsatellite Genetic Distances: Impacts of Population Growth and Gene Flow

Lev A. Zhivotovsky

N. I. Vavilov Institute of General Genetics, Russian Academy of Sciences, Moscow, Russia


    Abstract
 TOP
 Abstract
 Introduction
 Methods and Results
 Discussion
 Appendix
 Acknowledgements
 literature cited
 
Genetic distances play an important role in estimating divergence time of bifurcated populations. However, they can be greatly affected by demographic processes, such as migration and population dynamics, which complicate their interpretation. For example, the widely used distance for microsatellite loci, ({delta}µ), assumes constant population size, no gene flow, and mutation-drift equilibrium. It is shown here that ({delta}µ) strongly underestimates divergence time if populations are growing and/or connected by gene flow. In recent publications, the average estimate of divergence time between African and non-African populations obtained by using ({delta}µ) is about 34,000 years, although archaeological data show a much earlier presence of modern humans out of Africa. I introduce a different estimator of population separation time based on microsatellite statistics, TD, that does not assume mutation-drift equilibrium, is independent of population dynamics in the absence of gene flow, and is robust to weak migration flow for growing populations. However, it requires a knowledge of the variance in the number of repeats at the beginning of population separation, V0. One way to overcome this problem is to find minimal and maximal bounds for the variance and thus obtain the earliest and latest bounds for divergence time (this is not a confidence interval, and it simply reflects an uncertainty about the value of V0 in an ancestral population). Another way to avoid the uncertainty is to choose from among present populations a reference whose variation is presumably close to what it might have been in an ancestral population. A different approach for using TD is to estimate the time difference between adjacent nodes on a phylogenetic population tree. Using data on variation at autosomal short tandem repeat loci with di-, tri-, and tetranucleotide repeats in worldwide populations, TD gives an estimate of 57,000 years for the separation of the out-of-Africa branch of modern humans from Africans based on the value of V0 in the Southern American Indian populations; the earliest bound for this event has been estimated to be about 135,000 years. The data also suggest that the Asian and European populations diverged from each other about 20,000 years, after the occurrence of the out-of-Africa branch.

Key Words: modern human divergence • microsatellite genetic distance • separation time • population growth • migration flow

Key Words: modern human divergence • microsatellite genetic distance • separation time • population growth • migration flow


    Introduction
 TOP
 Abstract
 Introduction
 Methods and Results
 Discussion
 Appendix
 Acknowledgements
 literature cited
 
Advances in DNA technology have facilitated the study of population history, and especially the study of the history of modern humans. Most of these studies strongly support the out-of-Africa theory (Andrews 1986Citation ; Stringer and Andrews 1988Citation ), which posits that modern humans emerged recently from Africa, split into major continental groups, and replaced existing hominid populations throughout the Old World with no significant genetic admixture (Cann, Stoneking, and Wilson 1987Citation ). However, there is the problem of how to estimate the time of separation of the branches on a population tree using genetic variation.

Two different methodologies have been developed for dating ancient evolutionary events, one based on a concept of genetic distances, and the other using haplotype data. The former is a "population" approach in which a certain demographic model of population bifurcation is assumed from archaeological and other findings, genetic distances are computed between the branches on the phylogenetic population tree, and, finally, the time since the divergence event is estimated using knowledge of how the genetic distance depends on mutation rate and other parameters (Cavalli-Sforza 1998Citation ). The latter methodology constructs a haplotype tree and estimates the time of the most recent common ancestor of the corresponding DNA region at which genetic variation is represented by the haplotypes. No recombination is assumed, so this method is relevant to mitochondrial DNA, the Y chromosome, and short regions in the autosomes and the X chromosome. Hey (1998)Citation argues that the population methodology is not able to estimate the time since bifurcation in the presence of gene flow, and instead suggests studying haplotypes as a key to genetic history. It might be added that the population approach is also complicated by population dynamics: genetic distances may deviate from their theoretically predicted values if population size is changing (Reich, Feldman, and Goldstein 1999Citation ).

The purpose of this paper is to propose an estimation procedure for divergence time that includes the possibilities of population size change and gene flow. We focus on microsatellite variation because it has been widely used in population studies. It is shown that the genetic distance ({delta}µ)2 (Goldstein et al. 1995bCitation ) greatly underestimates the time of divergence if populations are growing in size or are connected by weak gene flow. I introduce a different estimator of divergence time, TD, that does not depend on population dynamics and is sufficiently robust to weak gene flow. Analysis of data on di-, tri-, and tetranucleotide microsatellite loci results in higher values of divergence times between major human populations than previous estimates based on ({delta}µ)2.


    Methods and Results
 TOP
 Abstract
 Introduction
 Methods and Results
 Discussion
 Appendix
 Acknowledgements
 literature cited
 
The Influence of Demographic Processes on ({delta}µ)2 and TD
TD–Estimator of Divergence Time
Variation at microsatellite loci is amenable to the use of population statistics that treat allele repeat scores (the number of repeats) as a quantitative trait (Goldstein et al. 1995a, 1995bCitation ; Slatkin 1995Citation ; Zhivotovsky and Feldman 1995Citation ). Among them is the genetic distance ({delta}µ)2 (Goldstein et al. 1995bCitation ), the squared difference in the number of repeats, which is linear with time since divergence (Goldstein et al. 1995bCitation ; Nei 1995Citation ; Zhivotovsky and Feldman 1995Citation ; Cooper et al. 1999Citation ). That is, the expectation of ({delta}µ)2 is 2w{tau} if the populations are at mutation-drift equilibrium and have constant effective size. Here, {tau} is the number of generations since divergence, and w is the effective mutation rate, that is, the product of the mutation rate and the second moment of mutational change in the number of repeats (eq. 5 , in the appendix). Thus, ({delta}µ)2/2w is an estimator of divergence time {tau} since a population bifurcated into two.

I suggest a different estimator of separation time:


(1)
where 1 is the average over loci of the average squared difference between pairs of alleles sampled, one each from the two populations (Goldstein et al. 1995aCitation ), and 0 is the average over loci of the within-population variance in the number of repeats in the ancestral population prior to its subdivision; w is the effective mutation rate. TD estimates a divergence time in terms of generations. The statistical procedure for computing TD and its evolutionary (genetic) error are given in the appendix.

In order to investigate the extent to which population dynamics and gene flow affect both distances, a set of numerical computations was carried out using recursions (7) and (8), in the appendix.

Population Growth and ({delta}µ)2
If populations are growing in size, ({delta}µ)2/2w underestimates the actual divergence time (fig. 1a ). The faster the growth, the stronger the deviation from actual divergence time (fig. 1 ). A decline in population size results in overestimation of actual divergence time (fig.1a, triangles) although the effects are not pronounced.



View larger version (53K):
[in this window]
[in a new window]
 
Fig. 1.—Expectations of divergence time between populations with no gene flow based on (a) ({delta}µ)2 and (b) TD. Numerical iterations are based on equations (7) , with w = 0.001. One time unit is 1,000 generations. Population growth follows a logistic law, equation (9) , with the initial size N0 = 1,000, the rate of growth a = 0.001, and the ultimate I-fold increases in size (plus signs) I = 1 (no growth), (diamonds) I = 10, (boxes) I = 100, and (triangles) I = 0.1 (10-fold decline). TD does not depend of population dynamics in the absence of gene flow, and thus the curves on b coincide

 
Gene Flow and ({delta}µ)2
Gene flow between two bifurcated growing populations may influence ({delta}µ)2 even more dramatically, up to its nonmonotonic dependence on the time since bifurcation (figs. 2a and 3a ). Note that in the presence of gene flow, ({delta}µ)2/2w may eventually underestimate divergence times even if both populations decline in size: compare figures 2a and 3a (triangles) with figure 1a (triangles).



View larger version (51K):
[in this window]
[in a new window]
 
Fig. 2.—Expectations of divergence time between populations with gene flow, m = 0.0001, based on (a) ({delta}µ)2 and (b) TD. The parameters and notation are as in figure 1 , and additionally, x = a reference line for the actual divergence time

 


View larger version (49K):
[in this window]
[in a new window]
 
Fig. 3.—Expectations of divergence time between populations with gene flow, m = 0.001, based on (a) ({delta}µ)2 and (b) TD. The notation is as in figures 1 and 2

 
Population Growth and TD
It is proved in the appendix (eq. 8 ) that, unlike ({delta}µ)2, TD gives an unbiased estimate of the time of separation regardless of change in population size (see fig. 1b ).

Gene Flow and TD
Unlike ({delta}µ)2, TD only slightly underestimates divergence time if growing populations are connected by weak gene flow. Moreover, the faster the populations are growing, the smaller the deviation of the predicted divergence time from its actual value (figs. 2b and 3b ). The error in TD is relatively small even if populations are expanding slowly in size. This is illustrated by figures 2b and 3b (diamonds) in the case of an average increase in size by less than one individual per generation (from 1,000 to 10,000 over 10,000 generations). In the case of constant or decreased population size, deviation of the divergence time predicted by TD from the actual one can be great, although much less than that for ({delta}µ)2.

Estimating TD
The statistical procedure for estimating TD and its error is given in the appendix. As equation (1) shows, computing TD requires two different steps: estimating 1 and estimating 0. The first one is a purely statistical problem and uses data on current populations. The second one is less certain in most cases, because an ancestral population will not be available, and thus 0 cannot be directly estimated from data. The following approaches are suggested to bound or approximate it.

Earlier and Later Bounds for the Divergence Time
One approach is to find a bound for 0, and so for TD. An earlier bound for divergence time can be established by setting 0 to a certain minimum, V0min:


(2)
As an extreme, V0min = 0 can be used. In this case, the value of 1/2w gives the earliest bound in terms of generations. A later bound can be constructed in a similar way, by setting 0 to its possible maximal value, V0max:


(3)
It follows from equation (1) and equation (5) (in the appendix) that ({delta}µ)2/2w can be viewed as TD with 0 taken as the average of the within-population variances in the two populations under comparison. The divergence time estimated with ({delta}µ)2 is less than that estimated with TD if within-population variance in the two current populations exceeds its value in an ancestral population. This is why ({delta}µ)2 underestimates divergence time for growing populations—their variation increases with time.

Using Current Values of V
If there are data on several populations that are supposed to have a common ancestor, available archaeological and demographic information might suggest that some of these populations can serve as a reference, that is, they might have the same level of variation as the ancestor had. This approach can give a later bound for divergence time if the reference populations descended from a large growing population with high genetic variation and have not yet achieved the ancestral value of variance.

Difference in Divergence Time Between Two Adjacent Nodes
One might estimate a time difference between two subsequent, adjacent branch nodes on a population tree, assuming equal values of the within-population variance at the nodes. That is, if '1 and ''1 correspond to the two adjacent nodes separated at unknown times T' and T'', respectively, for which within-population variances in the number of repeats at microsatellite loci are assumed to be equal, the time difference between the nodes, {Delta}T = T'' - T', is simply


(4)
This estimate of {Delta}T can be used if the within-population variances at T' and T'' differ little from each other, for example, if the adjacent nodes are not separated by a long period. Otherwise, as follows from equation (1) , this procedure underestimates the time difference between adjacent nodes if the average within-population variance at the later node is greater than that at the preceding node, and vice versa.

Analysis of Human Data
In order to compare the performance on data of TD with that of ({delta}µ)2, we analyzed five microsatellite data sets. The first set was composed of 28 loci with dinucleotide repeats of the initial 30 loci collected and analyzed by Bowcock et al. (1994)Citation and used with ({delta}µ)2 by Goldstein et al. (1995b)Citation . One of the 30 loci, FES, was a tetranucleotide (Bowcock et al. 1994Citation ). Variation at another locus, D13S133, was so high that it constituted a large part of total variation, and the values of within-population variances in repeat numbers occurred as outliers in the total distribution of variances (Zhivotovsky et al. 2000Citation ). These two loci had been removed in our previous analyses (Zhivotovsky et al. 2000Citation ), as well as in that of Jin et al. (2000)Citation . Therefore, we also removed them from the current analysis. The second and the third data sets contained 22 tri- and 21 tetranucleotide loci, respectively (A. Bowcock and L. Bennett, personal communication). A list of loci and more detailed information on these three data sets are given in Zhivotovsky et al. (2000)Citation . The fourth data set was originally presented by Jorde et al. (1997)Citation and includes allele frequencies at 60 tetranucleotide loci. The fifth data set comprised the Y chromosome loci with tetranucleotide repeat motifs (Seielstad et al. 1999Citation ). For the latter data, the initial set of eight short tandem repeats (STRs) was reduced by removing DYS388 and DYS392 (with trinucleotide repeats) and complex loci DYS389a and DYS389a because they are highly unstable (Forster et al. 2000Citation ). The "evolutionary" mutation rates at the remaining four loci (DYS19, DYS390, DYS391, and DYS393) have been recently estimated (Forster et al. 2000Citation ) and were used in our study.

The five data sets contained samples from African, Asian, and European populations. Data on autosomal loci with di- (Bowcock et al. 1994Citation ), tri-, and tetranucleotide repeats (A. Bowcock and L. Bennett, personal communication) included African populations (Mbuti pygmies from Zaire, Biaka pygmies from the Central African Republic, and a sample from Lisongo), European populations (Italians and Northern Europeans), and Asian populations (Chinese, Japanese, and Cambodians) (see the map in Barbujani et al. 1997Citation ). Data on autosomal loci with tetranucleotide motifs (Jorde et al. 1997Citation ) included African populations (Mbuti and Biaka pygmies, Sotho-Tswana, Tsonga, Nguni, San), European populations (Northern Europeans, French, Poles, Finns), and Asian populations (Chinese, Japanese, Cambodians, Malay, Vietnamese). Data on Y-chromosome loci with tetranucleotide repeats (Seielstad et al. 1999Citation ) included African populations (including pygmies from Zaire and the Central African Republic and samples from Lissongo and other sites of Eastern and Central Africa), European populations (Basques, Catalans, and a combined sample from Italy and Germany), and Asian populations (Chinese, Japanese, Cambodians).

Three kinds of estimates were obtained with these data to estimate divergence time between Africans and non-Africans: (1) the earliest bound, TD with 0 = 0, which overestimates an actual divergence time if within-population variation at microsatellite loci in an African ancestor was not sufficiently small compared with variation in the current populations; (2) ({delta}µ)2, which underestimates an actual divergence time for growing populations; and (3) TD with the value of 0 in equation (1) based on genetic variation in Southern American Indian populations (samples from Karitiana and Surui; see below more discussion on choosing these populations as a reference). In addition to the three statistics, I computed {Delta}T, the time difference between two major events in modern human evolution: separation of the non-African branch and its division into Asian and European populations.

As mentioned above, the Southern American Indian populations (Karitiana and Surui) can be considered a reference for the ancient value of variance: the average values of the within-population variances taken across di-, tri-, and tetranucleotide loci in each of the first three data sets in the two populations were 4.416, 2.522, and 2.591, respectively, which can be taken as 0 in formula (1) for these three data sets. Note that the values were 0.661, 0.671, and 0.652 times the average variances for pygmy populations (Biaka and Mbuti), 6.676, 3.759, and 3.972, respectively, in the three data sets. The fourth data set, with tetranucleotide loci, did contain samples from the same African populations, Biaka and Mbuti pygmies (with a within-population variance of 4.973), but did not include samples from South Amerindians. To make a reference for the ancient variation for the fourth data set, I extrapolated to it the average of the above ratios, (0.661 + 0.671 + 0.652)/3 = 0.661, and computed the corresponding 0 as simply 4.973 x 0.661. The fifth data set included the Southern American Indian populations, and thus 0 in equation (1) could be obtained directly.

For each locus, each of the statistics, TD, {Delta}T, and ({delta}µ)2, was averaged over all possible pairs of populations under comparison. For each intercontinent comparison, the mean value of each statistic and its statistical error were computed over loci (see appendix). The effective mutation rate, w, was taken as 1.52 x 10-3, 0.85 x 10-3, and 0.93 x 10-3 for autosomal di-, tri-, and tetranucleotide loci, respectively (see table 1 of Zhivotovsky et al. 2000Citation ), with a mutation rate of 0.26 x 10-3 estimated for the Y-chromosome tetranucleotides (Forster et al. 2000Citation ). To convert the number of generations into years, 25 years was taken as a generation time for modern humans.


View this table:
[in this window]
[in a new window]
 
Table 1 Estimates of Divergence Time Between Africans and Non-Africans and Between Asians and Europeans (in thousands of years)

 
Table 1 shows that for the autosome microsatellite data sets, TD revealed similar estimates for divergence times between Africans and non-Africans and indicated that the out-of-Africa branch of modern humans bifurcated into Asians and Europeans about 20,000 years after the out-of-Africa branch had separated from the African ancestor (fig. 4 ). The Y-chromosome data had large statistical errors because of a small number of loci and did not show significant deviation from the estimates suggested by autosomal markers.



View larger version (16K):
[in this window]
[in a new window]
 
Fig. 4.—TD estimates for two major events in the history of modern humans: the occurrence of the out-of-Africa branch and its division into Asian and European populations

 

    Discussion
 TOP
 Abstract
 Introduction
 Methods and Results
 Discussion
 Appendix
 Acknowledgements
 literature cited
 
The Performance of ({delta}µ)2
The genetic distance ({delta}µ)2 is unbiased only if both populations are (1) at mutation-drift equilibrium, (2) of identical constant sizes during their evolution, and (3) not subject to gene flow. This is probably not the case for most species, including human populations whose total size has been growing and which might have been connected by gene flow (Cavalli-Sforza, Menozzi, Piazza 1996Citation , pp. 105–107; Hey 1998Citation ). The two demographic factors, population growth and gene flow, cause strong underestimation of divergence times by ({delta}µ)2 (figs. 1a, 2a, and 3a ). Moreover, phylogenetic trees based on ({delta}µ)2 may have incorrect topologies because the distance is a nonmonotonic function of time for a certain degree of gene flow and population growth (figs. 2a and 3a ).

Table 1 shows the earliest bound for divergence time between African and non-African populations estimated using 0 = 0 with autosomal STRs: about 135,000 years. Therefore, the actual divergence time value should be smaller. Probably, the estimate based on 0 for Southern Amerindians, 57,000 years (table 1 ), gives its later bound. Indeed, among the current human populations, South American aboriginals can be considered a reference for microsatellite variation in an ancient African ancestor because their population size is low and might be compared with that estimated for an African ancestor, from one to a few thousand gametes (Rogers and Harpending 1992Citation ; Rogers 1995Citation ; Rogers and Jorde 1995Citation ; Zhivotovsky et al. 2000Citation ), and they have maintained their style of life probably since they arrived in this area (L. L. Cavalli-Sforza, personal communication). On the other hand, the Southern American Indian populations have descended from Asians, whose variation is very high, and thus the time since human colonization of South America might not have been sufficient to reduce microsatellite variation to the level of that in the ancestral population in Africa. Also, gene flow between even small villages can greatly increase within-population variance at microsatellite loci (Feldman, Kumm, and Pritchard 1999Citation ). Nevertheless, there is no objective proof of whether or not variation at the studied loci in South Amerindians is close to an African ancestral value. Therefore, the corresponding estimates of TD should be taken cautiously.

It should be kept in mind that the earliest and latest bounds reflect our uncertainty in the ancient value of the variance V0 and cannot be regarded as a confidence interval because they refer to some extreme values of 0, whereas confidence intervals highlight sampling errors. It is well known that in estimating procedures that involve more than one locus, there are two sources of statistical errors: the between-locus evolutionary (or genetic) error due to sampling loci, and the within-locus error caused by multinomial sampling of individuals (Weir 1996Citation ). Sampling error is relatively negligible if the sample size, N, and the number of loci, L, are sufficiently large, because it is based on LN observations. In contrast, the evolutionary error is based on L observations and thus is the major source of statistical errors (its algorithm is given in the appendix). Table 1 shows that the evolutionary errors for our estimates are still large, even for combined data with 131 loci. This corresponds well to our conclusion (Zhivotovsky and Feldman 1995Citation ; Zhivotovsky et al. 2000Citation ) that hundreds of microsatellite loci might be needed to obtain statistically stable estimates.

({delta}µ)2 suggests extremely recent subdivision of Africans and non-Africans, about 29,000 years for autosomal microsatellite loci, which is at least half of the TD estimate based on the reference Southern American Indian populations, 57,000 years (table 1 ). Since archaeological data, as well as mtDNA, show earlier divergence of modern humans more than 50,000, or even more than 100,000, years ago (Cann, Stoneking, and Wilson 1987Citation ; Roberts, Jones, and Smith 1990Citation ; Quintana-Murci et al. 1999Citation ), such a low estimate by ({delta}µ)2 is probably a symptom of statistical bias for ({delta}µ)2 caused by population growth and gene flow. Indeed, the theory (fig. 1 ) predicts that ({delta}µ)2 gives about half the actual divergence time in 2,000 generations after divergence (or 50,000 years for humans if 25 years is used as the generation time length) even if the populations were slowly increasing in size; ({delta}µ)2 shows even greater deviation down from the actual divergence time if, additionally, migration presents (figs. 2 and 3 ).

Recently published results also give low values of divergence time inferred using of ({delta}µ)2. Although Goldstein et al. (1995b)Citation estimated the age of non-Africans as ~156,000 years, assuming one-step mutations for dinucleotide loci with a mutation rate of 0.00065, it can be argued that a multistep mutation model might be more reasonable for microsatellite loci with dinucleotide repeats. (It is worth noting that a highly variable locus, D13S133, had also contributed to the high value of this figure, 156,000 years; see Analysis of Human Data, above). Using the same dinucleotide data set, but without D13S133 and with multistep mutation parameters inferred from experiments on more than 5,000 dinucleotide loci (Dib et al. 1996Citation ), a much more recent estimate of divergence time between Africans and non-Africans, 39,000 years, was obtained (Jin et al. 2000Citation ); the latter differs from our estimate, 36,000 years (table 1 ) because the former used additional data on populations from Oceania and South America. For the same samples, Jin et al. (2000)Citation used a larger set of dinucleotide loci and found an even smaller value, 23,000 years since the separation Africans and non-Africans, and Feldman, Kumm, and Pritchard (1999)Citation estimated the divergence time as 39,000 years based on both di- and tetranucleotide loci. Cooper et al. (1999) found about 43,000 and 54,000 years with two different estimation procedures. Recent estimates of the time of splitting into Africans and non-Africans were obtained for Y-chromosome microsatellites: 12,000 and 25,000 years, depending on assumptions for mutation rate (Seielstad et al. 1999Citation ). The average over these values is 34,000 years, close to our lowest bounds inferred by using ({delta}µ)2, 28,700 and 31,900 years, obtained with the use of the autosomal and the Y-chromosome markers (table 1 ).

The Advantages and Disadvantages of the TD Estimator
The advantages of TD are that it (1) does not assume mutation-drift equilibrium, (2) does not depend on population dynamics (fig. 1 ), and (3) is sufficiently robust to gene flow for growing populations (Figs. 2b and 3b ). Therefore, TD can be useful for revealing bifurcation time from population genetic data even if populations have been changing in size and/or connected by gene flow since their divergence.

A major problem with application of TD is that the estimator heavily relies on knowledge of the initial variance in the number of repeats at the time of separation. This is a crucial point in the application of TD because, in most cases, 0 will not be known a priori. Nevertheless, a few workable approaches have been suggested above to overcome this problem: use of the earlier and later bounds for divergence time, use of a reference population that models variation in an ancestor, and estimation of the time separating adjacent nodes on a population tree. One more way to evaluate 0 might be to estimate effective population size in the initial generation, N0, by using models of growing populations (Rogers 1995Citation ; Rogers and Jorde 1995Citation ; Kimmel et al. 1998Citation ; Reich and Goldstein 1998Citation ; Pritchard et al. 1999Citation ; Zhivotovsky et al. 2000Citation ). However, to compute 0, one has to assume that the ancestral population was at mutation-drift equilibrium prior its subdivision, V0 {approx} 2wN0. It is worth noting, however, that it is the initial variance, not the initial population size, that is a prime term of TD in equation (1) . Therefore, it is necessary to develop an estimating procedure for ancient values of the within-population variance in allele repeat scores.

Another problem with application of TD is the necessity of knowledge of the effective mutation rate, w. That mutation rate may strongly depend on a locus has been demonstrated by Forster et al. (2000)Citation , whose estimate for the Y-chromosome loci used in our study was 10-fold less than the estimates involving other loci. There are few experimental data that permit direct estimation of mutation parameters for microsatellites. Dinucleotide markers in humans are now most appropriate for this because there already exist data on more than 5,000 dinucleotide polymorphisms (Dib et al. 1996Citation ). Using these and population data, Feldman, Kumm, and Pritchard (1999)Citation and Zhivotovsky et al. (2000)Citation estimated w for human microsatellites. It is worth noting that mutation parameters may vary across microsatellite loci, so the average effective mutation rates for a given set of loci may deviate from the total value. Suppose that the effective mutation rate varies across loci with mean w and variance {sigma}2w. Since w is the denominator in TD, using L microsatellite loci brings a bias with the factor 1 + C2w/L, where C2w = {sigma}2w/w2 is the squared coefficient of variation of w (see the appendix). However, very little is known about the amount of variation in effective mutation rate. As follows from table 1 in Di Rienzo et al. (1998)Citation , the coefficients of variation in the frequencies of somatic mutations at di-, tri-, and tetranucleotide loci are <1: 0.31, 0.41, and 0.39, respectively. Using population data, Zhivotovsky, Goldstein, and Feldman (unpublished data) estimated Cw for human microsatellites with di- and tetranucleotide repeats to be <1. With such variation among microsatellites, several dozens of loci make the relative bias, C2w/L, negligible. As for the evolutionary error of divergence time, which is computed as a statistical error due to sampling among loci, it may substantially depend on variation in mutation rates (see the appendix). Therefore, the evolutionary errors in table 1 must be even greater if variation in mutation rates is high. It should be emphasized that such an influence of variation in mutation rate on estimates of divergence time and their evolutionary errors is not a property of only TD. Other genetic distances are also affected in the same way if calibrated by mutation rate.

In the case of unknown w, a possible way for estimating divergence time would be cancelling w by using the ratio statistic Bt/Vt (Chakraborty and Nei 1982Citation ; Kimmel et al. 1996Citation ). This statistic is defined as the between-population variance in the number of repeats subdivided by the within-population variance and thus is formally identical to the analogous index TR (Slatkin 1995Citation ), the ratio of the between-population and the total variances. The expected value of Bt /Vt is linear with time and equals {tau}/N, where {tau} is the number of generations since the population split and N is the population size. However, this statistic can be used for estimating divergence time only with the following assumptions: (1) complete reproductive isolation, (2) constant population size, and (3) mutation-drift equilibrium for their entire history. Moreover, although w need not be known, N has to be known for computing {tau} with Bt/Vt. Generally, estimates of divergence time inferred from the ratio statistics have a tendency to substantially deviate from their theoretical prediction (Hedrick 1999Citation ).

One more problem with using TD is the assumption that gene flow between diverged populations is weak. High migration rates may complicate the behavior of TD (data not shown). Therefore, the estimator has to be carefully applied to intracontinental subdivisions of modern humans which might be fitted by an island model better than a bifurcation scheme.

Also, it should be kept in mind that constant population size or decrease in size in the presence of gene flow may cause underestimation of divergence time by TD, although not as extensive as the underestimation by ({delta}µ)2 (figs. 2 and 3 , plus signs and triangles). In this regard, it is important to note that even a slow decline in size of one of two diverged connected populations overrides the effect of the increase of the other one (data not shown). Therefore, one should be careful in interpreting the TD values obtained for populations in which a decrease in size is suspected.

It is worth noting that ({delta}µ)2 is not the only distance measure influenced by demographic processes. For the infinite-alleles model, Nei's (1987Citation , pp. 232–238) distance is more appropriate, and Cavalli-Sforza and Edwards' (1967)Citation distance better fits the process of population separation by genetic drift alone. However, like ({delta}µ)2, these distances are also biased if there is gene flow and change in population size. Therefore, to what extent the actual divergence times deviate from those predicted by these models is also of great concern. The major advantage of TD is that it is a nearly unbiased estimator of divergence time robust to population growth and gene flow and thus allows for the concept of genetic distances in dating ancient population divisions.

Despite the above-mentioned difficulties, TD gives satisfactory estimates of divergence time for major divisions in the history of modern humans (fig. 4 ). Using 131 autosomal STRs, it bounds the separation of African and Eurasian populations as early as 135,000 years ago (table 1 ). The separation time should exceed the values obtained with the use of ({delta}µ)2, 29,000 years, because the major human populations have generally been growing in size, which strongly decreases estimates based on ({delta}µ)2 (see fig. 1 ); the possibility of exchange by migrants can decrease these estimates even more strongly (figs. 2 and 3 ). The estimate of 57,000 years revealed by TD with use of microsatellite variation in the Southern American Indian populations as a reference might be considered a possible approximation for divergence time between the African and the non-African branches.

Although the uncertainty caused by using 0 in TD seems to be large, it is usual for estimates of evolutionary times. Moreover, the range from 57,000–135,000 years for divergence time between Africans and non-Africans revealed by TD for autosomal STR markers lies within the wider range of the existing estimates. Namely, the TD estimate of divergence time between Africans and non-Africans based on the variance in the number of repeats in Southern American Indian populations is larger than estimates based on ({delta}µ)2, from 12,000 years (Seielstad et al. 1999Citation ) to 54,000 years (Cooper et al. 1999Citation ), whereas the earliest bound, 134,600 years (table 1 ), is lower than the estimate based on the X-linked PDHA1 locus (Harris and Hey 1999Citation ), 200,000 years, and close to the 150,000 and 137,000 years inferred from data on 62 polymorphic loci (Nei 1987Citation , table 9.4) and Alu insertions (Stoneking et al. 1997Citation ), respectively.

It is important to emphasize that the use of TD allowed dating of the time of division of the out-of-Africa branch into Asian and European populations (fig. 4 ): the adjacent-nodes approach shows that the Asian and European populations had separated from each other about 20,000 years after the occurrence of the out-of-Africa branch of modern humans.


    Appendix
 TOP
 Abstract
 Introduction
 Methods and Results
 Discussion
 Appendix
 Acknowledgements
 literature cited
 
Mutation Scheme
Here, I consider a general multistep mutation scheme with constant bias {theta} that can be zero (no bias), positive, or negative. Let vc be the probability that an allele mutates to another allele whose number of repeats is larger by c (if c is negative, the mutation has a smaller size); v0 is the probability of having no mutation. The total mutation rate, µ, is {Sigma}s,s!=0 vs, or 1 - v0. Introduce the first four central moments in mutation changes


(1)
Neglecting terms of order µ2 and smaller, we have


(5)
where {eta}(2)m, {eta}(3)m, and {eta}(4)m are noncentral moments of mutational changes in allele size (Di Rienzo et al. 1998Citation ). In a special case of no mutation bias, {eta}(2)m is {sigma}2m, the variance in mutation changes (Slatkin 1995Citation ).

Dynamic Equations for ({delta}µ)2 and D1
Consider a diploid population of size N. With the above notation, equations (3)–(8) by Zhivotovsky and Feldman (1995)Citation remain valid (with the obvious replacement of N by the number of gametes, 2N) (see also Feldman, Kumm, and Pritchard 1999Citation ). In particular,


(6)
where r and V are the population mean and variance, respectively; a prime sign denotes the progeny generation. It can be shown that if two populations A and B are connected by gene flow at rate m, these relationships imply the following recursions:


(7)
where {tau} is time in generations. This system generalizes equations (8.14) and (8.16) of Feldman, Kumm, and Pritchard (1999) to a multistep mutation scheme for the case of two populations of unequal size.

Since D1 is defined as ({delta}µ)2 + VA + VB (Goldstein et al. 1995bCitation ), summing up equations (7) gives the following recursion for D1 in the case of no migration: D'1 = D1 + 2w. This implies that in the absence of gene flow,

(8)
which is an extension of equation (A8) of Goldstein et al. (1995a)Citation to the case of multistep mutation and population dynamics. Equation (8) underlies the definition of TD in equation (1) and explains why TD does not depend on population dynamics in the absence of gene flow (fig. 1b ). In the numerical calculations, we assume logistic population growth of the form

(9)
where a is the growth rate and I is the factor by which the population size ultimately will exceed the initial size N0. Replacing the differences V(t+1)A - V(t)A, etc. in equation (7) with the appropriate derivatives, we numerically solve the corresponding system of differential equations with Mathematica (Wolfram 1996Citation , pp. 886–894).

A Bound for the Initial Variance and the Initial Effective Population Size
Note that equation (1) can also be used to obtain a bound for the initial variance, V0, prior to population subdivision. Indeed, if we set TD in equation (1) to a definite value that can be chosen as a later or an earlier bound for divergence time (in terms of generations) assumed to be known from some sources, then the corresponding upper or lower bound for the initial variance in allele repeat scores is

(10)
If one assumes mutation-drift equilibrium prior to subdivision, then the corresponding bound for the initial population size can be estimated as


(11)

Sampling Estimate of TD
Let nAi and nBi be the counts of the allele with i repeats at a certain locus in samples from populations A and B, respectively. The allelic sample sizes in the populations are nA = {Sigma}i nAi and nB = {Sigma}i nBi, respectively; note that nA = 2NA and nB = 2NB. Then, we compute the average numbers of repeats in the samples, A and B, and the biased sampling estimates of the within-population variances, A and B (the tilde indicates sampling estimates):


(12)
An unbiased estimate of 1 at each locus is then computed as

(13)
For each locus, we compute tD = 1 - 2V0, where V0 is the value of within-population variance in the number of repeats in the ancestral population prior to its subdivision. Next, we compute across loci both the arithmetic mean of these single-locus estimates, tD, and its standard statistical error, stD; the latter is known to be a bootstrap error for an arithmetic mean (Efron and Tibshirani 1993Citation , pp. 42–43). Finally, TD = tD/2w, and its evolutionary (or genetic; Weir 1996Citation ) error is sTD = stD/2w. Note that the given estimating procedures for TD and sTD assume no variation in mutation rates among loci. If effective mutation rates vary across loci with mean w and variance {sigma}2w, we need to consider the expectation and the variance of tD/2; here, the bar indicates averaging across L loci involved in the study. Using the delta method (e.g., Weir 1996Citation ), it can be obtained that the expectation of TD is D/2w x (1 + C2w/L), where C2w = {sigma}2w/w2 is the squared coefficient of variation of w, and D is the expectation of tD. For the evolutionary error sTD, the same method gives the square of evolutionary error for TD as s2TD + T2DC2w, where s2TD is the above given estimate of evolutionary error assuming no variation in w.


    Acknowledgements
 TOP
 Abstract
 Introduction
 Methods and Results
 Discussion
 Appendix
 Acknowledgements
 literature cited
 
I am grateful to Dr. Luca Cavalli-Sforza, Dr. Marc Feldman, Dr. David Goldstein, Dr. Jeffrey Long, and two anonymous reviewers for helpful comments and constructive suggestions, and to Dr. Anne Bowcock, Dr. Lynda Bennett, and Dr. Lyn Jorde for providing me with their data on variation at microsatellite loci in worldwide human populations. This research was supported in part by the Russian Foundation of Basic Research (grant 98-04-49292), the National Institutes of Health (grants 1 R03 TW00491-01, GM 28016, and GM 28428), the N. I. Vavilov Institute of General Genetics, the Morrison Institute for Population and Research Studies at Stanford University, and the Centre for Human Genetics at Edith Cowan University.


    Footnotes
 
Jeffrey Long, Reviewing Editor

2 Address for correspondence and reprints: Lev A. Zhivotovsky, N. I. Vavilov Institute of General Genetics, Russian Academy of Sciences, 3 Gubkin Street, Moscow 117809, Russia. lev{at}vigg.ru Back


    literature cited
 TOP
 Abstract
 Introduction
 Methods and Results
 Discussion
 Appendix
 Acknowledgements
 literature cited
 

    Andrews, P. 1986. Fossil evidence on human origins and dispersal. Cold Spring Harb. Symp. Quant. Biol. 51:419–428.[ISI][Medline]

    Barbujani, G., A. Magagni, E. Minch, and L. L. Cavalli-Sforza. 1997. An appointment of human DNA diversity. Proc. Natl. Acad. Sci. USA 94:4516–4519.

    Bowcock, A. M., A. Ruiz-Linares, J. Tomfohrde, E. Minch, J. R. Kidd, and L. L. Cavalli-Sforza. 1994. High resolution of human evolutionary trees with polymorphic microsatellites. Nature 368:455–457.

    Cann, R. L., M. Stoneking, and A. Wilson. 1987. Mitochondrial DNA and human evolution. Nature 325:31–36.

    Cavalli-Sforza, L. L. 1998. The DNA revolution in population genetics. Trends Genet. 14:60–65.[ISI][Medline]

    Cavalli-Sforza, L. L., and A. W. F. Edwards. 1967. Phylogenetic analysis: models and estimation procedures. Am. J. Hum. Genet. 19:223–257.

    Cavalli-Sforza, L. L., P. Menozzi, and A. Piazza. 1996. The history and geography of human genes. Princeton University Press, Princeton, N.J.

    Chakraborty, R., and M. Nei. 1982. Genetic differentiation of quantitative characters between populations of species: 1. Mutation and random genetic drift. Genet. Res. Camb. 39:303–314.

    Cooper, G., W. Amos, R. Bellamy, M. R. Siddiqui, A. Frodsham, A. V. S. Hill, and D. C. Rubinsztein. 1999. An empirical exploration of the ({delta}µ)2 genetic distance for 213 human microsatellite markers. Am. J. Hum. Genet. 65:1125–1133.[ISI][Medline]

    Di Rienzo, A., P. Donnelly, C. Toomajian, B. Sisk, A. Hill, M. L, Petzl-Erler, G. K. Haines, and D. H. Barch. 1998. Heterogeneity of microsatellite mutations within and between loci, and implications for human demographic histories. Genetics 148:1269–1281.

    Dib, C., S. Faure, C. Fizames et al. (14 co-authors). 1996. A comprehensive genetic map of the human genome based on 5,264 microsatellites. Nature 380:152–154.

    Efron, B., and R. J. Tibshirani. 1993. An introduction to the bootstrap. Chapman and Hall, N.Y.

    Feldman, M. W., J. Kumm, and J. K. Pritchard. 1999. Mutation and migration in models of microsatellite evolution. Pp. 98–115 in D. G. Goldstein and C. Schlotterer, eds. Microsatellites: evolution and applications. Oxford University Press, Oxford, England.

    Forster, P., A. Rohl, P. L. Lunnermann, C. Brinkmann, T. Zerjal, C. Tyler-Smith, and B. Brinkmann. 2000. A short tandem repeat-based phylogeny for the human Y chromosome. Am. J. Hum. Genet. 67:182–196.[ISI][Medline]

    Goldstein, D. B., A. R. Linares, L. L. Cavalli-Sforza, and M. W. Feldman. 1995a. An evaluation of genetic distances for use with microsatellite loci. Genetics 139:463–471.

    ———. 1995b. Genetic absolute dating based on microsatellites and the origin of modern humans. Proc. Natl. Acad. Sci. USA 92:6723–6727.

    Harris, E. E., and J. Hey. 1999. X chromosome evidence for ancient human histories. Proc. Natl. Acad. Sci. USA 96:3320–3324.

    Hedrick, P. H. 1999. Perspective: highly variable loci and their interpretation in evolution and conservation. Evolution 53:313–318.

    Hey, J. 1998. Population genetics and human origins—haplotypes are key! Trends Genet. 14:303–304.

    Jin, L., M. L. Baskett, L. L. Cavalli-Sforza, L. A. Zhivotovsky, M. W. Feldman, and N. A. Rosenberg. 2000. Microsatellite evolution in modern humans: a comparison of two data sets from the same populations. Ann. Hum. Genet. 64:117–134.[ISI][Medline]

    Jorde, L. B., A. R. Rogers, M. Bamshad, W. S. Watkins, P. Krakowiak, S. Sung, J. Kere, and H. Harpending. 1997. Microsatellite diversity and the demographic history of modern humans. Proc. Natl. Acad. Sci. USA 94:3100–3103.

    Kimmel, M., R. Chakraborty, J. P. King, M. Bamshad, W. S. Watkins, and L. B. Jorde. 1998. Signatures of population expansion in microsatellite repeat data. Genetics 148:1921–1930.

    Kimmel, M., R. Chakraborty, D. N. Stivers, and R. Deka. 1996. Dynamics of repeat polymorphisms under a forward-backward mutation model: within- and between-population variability at microsatellite loci. Genetics 143:549–555.

    Nei, M. 1987. Molecular evolutionary genetics. Columbia University Press, N.Y.

    Nei, M. 1995. Genetic support for the out-of-Africa theory of human evolution. Proc. Natl. Acad. Sci. USA 92:6720–6722.

    Pritchard, J. K., M. T. Seielstad, A. Perez-Lezaun, andM. W. Feldman. 1999. population growth of human Y chromosomes: a study of Y chromosome microsatellites. Mol. Biol. Evol. 16:1791–1798.[Abstract/Free Full Text]

    Quintana-Murci, L., O. Semino, H.-J. Bandelt, G. Passarino, K. McElreavey, and A. S. Santachiara-Benerecetti. 1999. Genetic evidence of an early exit of Homo sapiens sapiens from Africa through eastern Africa. Nat. Genet. 23:437–441.[ISI][Medline]

    Reich, D. E., M. W. Feldman, and D. B. Goldstein. 1999. Statistical properties of two tests that use multilocus data sets to detect population expansions. Mol. Biol. Evol. 16:453–466.[Free Full Text]

    Reich, D. E., and D. B. Goldstein. 1998. Genetic evidence for a Paleolithic human population expansion in Africa. Proc. Natl. Acad. Sci. USA 95:8119–8123.

    Roberts, R. G., R. Jones, and M. A. Smith. 1990. Thermoluminescence dating of a 50,000-year-old human occupation site in northern Australia. Nature 345:153–156.

    Rogers, A. R. 1995. Genetic evidence for a Pleistocene population expansion. Evolution 49:608–615.

    Rogers, A. R., and H. C. Harpending. 1992. Population growth makes waves in the distribution of pairwise differences. Mol. Biol. Evol. 9:552–569.[Abstract]

    Rogers, A. R., and L. B. Jorde. 1995. Genetic evidence on modern human origins. Hum. Biol. 67:1–36.[ISI][Medline]

    Seielstad, M., E. Bekele, M. Ibrahim, A. Toure, and M. Traore. 1999. A view of modern human origins from Y chromosome microsatellite variation. Genome Res. 9:558–567.[Abstract/Free Full Text]

    Slatkin, M. 1995. A measure of population subdivision based on microsatellite allele frequencies. Genetics 139:457–462.

    Stoneking, M., J. J. Fontius, S. L. Clifford et al. (11 co-authors). 1997. Alu insertion polymorphisms and human evolution: evidence for a larger population size in Africa. Genome Res. 7:1061–1071.[Abstract/Free Full Text]

    Stringer, C. B., and P. Andrews. 1988. Genetic and fossil evidence for the origin of modern humans. Science 239:1263–1268.

    Weir, B. 1996. Genetic data analysis II. Sinauer, Sunderland, Mass.

    Wolfram, S. 1996. Mathematica. A system for doing mathematics by computer. 3rd edition. Wolfram Media/Cambridge University Press, N.Y.

    Zhivotovsky, L. A., L. Bennett, A. M. Bowcock, and M. W. Feldman. 2000. Human population expansion and microsatellite variation. Mol. Biol. Evol. 7:757–767.

    Zhivotovsky, L. A., and M. W. Feldman. 1995. Microsatellite variability and genetic distances. Proc. Natl. Acad. Sci. USA 17:11549–11552.

Accepted for publication December 21, 2000.