Detecting Ancient Admixture in Humans Using Sequence Polymorphism Data

Corresponding author: Jeffrey D. Wall, University of Chicago, 1101 E. 57th St., Chicago, IL 60637., jdwall{at}midway.uchicago.edu (E-mail)

Communicating editor: S. TAVARÉ

	ABSTRACT

TOP ABSTRACT MODEL AND TEST STATISTICS RESULTS DISCUSSION LITERATURE CITED

A debate of long-standing interest in human evolution centers around whether archaic human populations (such as the Neanderthals) have contributed to the modern gene pool. A model of ancient population structure with recent mixing is introduced, and it is determined how much information (i.e., sequence data from how many unlinked nuclear loci) would be necessary to distinguish between different demographic scenarios. It is found that ~50–100 loci are necessary if plausible parameter estimates are used. There are not enough data available at the present to support either the "single origin" or the "multiregional" model of modern human evolution. However, this information should be available in a few years.

THE question of modern human origins has fascinated physical anthropologists for decades. Recently, two main hypotheses have been proposed. One view, often called the single origin model, claims that modern humans evolved in a single location (probably in Africa) roughly 150,000 years ago and from there expanded and replaced the existing hominid populations around the world (e.g., HOWELLS 1976 ; STRINGER and ANDREWS 1988 ). In contrast, the multiregional model claims no single location for the origin of modern humans; instead, it posits that over the last 1–2 million years humans have evolved gradually from their Homo erectus ancestors throughout the Old World (e.g., WEIDENREICH 1943 ; WOLPOFF et al. 1984 ). Both theories claim support from the fossil record (e.g., WOLPOFF et al. 1984 ; STRINGER and ANDREWS 1988 ), and it is hard to sort out their conflicting claims. However, recent fossil studies tend to support the single origin hypothesis (e.g., LAHR 1994 , LAHR 1996 ; HUBLIN et al. 1995 ; but see DUARTE et al. 1999 ), as do recent genetic studies (e.g., VIGILANT et al. 1991 ; ARMOUR et al. 1996 ; TISHKOFF et al. 1996 ). Despite this trend, there are still vocal supporters of the multiregional model (e.g., WOLPOFF 1996 ).

One way of rephrasing the debate is to ask how much contribution regional "archaic" populations have had to current human morphological and genetic diversity. Under the single origin model, this contribution is thought to be very small or nonexistent, whereas the multiregional model predicts that it is quite large. Other intermediate hypotheses predict a range of contributions from relatively small (e.g., BRAUER 1984 ) to relatively large (e.g., SMITH 1992 ).

One recent genetical breakthrough involved the recovery and sequencing of a fragment of Neanderthal mtDNA. The inferred sequence was found to be quite different from the homologous sequences of all extant humans (KRINGS et al. 1997 ). This observation has been taken as evidence for the claim that modern humans and Neanderthals did not interbreed, and thus that Neanderthals left no descendants among extant humans (e.g., KAHN and GIBBONS 1997 ; LINDAHL 1997 ; WARD and STRINGER 1997 ). However, it is impossible to make any strong inferences from single locus data because of the randomness inherent in the underlying evolutionary process. In fact, the mitochondrial data are consistent with a high level of interbreeding between humans and Neanderthals (NORDBORG 1998 ). Nordborg calculates that even if Neanderthal mtDNA comprised 25% of the ancestry of extant human mtDNA roughly 68,000 years ago, there would still be a >50% chance that by the present time all Neanderthal mtDNA would have been lost because of genetic drift. To have any power, we must have data from multiple unlinked loci.

In this article, I explore how much power we have to exclude the possibility that Neanderthals (or other archaic humans) left descendants among extant humans given sequence data from multiple unlinked neutral nuclear loci. I assume that no further loci from Neanderthals will be sequenced; given the difficulty encountered with ancient mtDNA, which is present in hundreds of copies in each cell, it is highly unlikely that any single copy nuclear DNA can be recovered from fossils as old as Neanderthals. Instead, I focus on the effect of past demographic events on the patterns of current observed variation in human populations. This can be done by examining the change in the shape of the genealogy caused by alternative demographic scenarios. Both changes in population size and the presence of population subdivision produce characteristic changes in the shape of genealogies when they are compared with genealogies generated from a constant size panmictic model (e.g., TAKAHATA 1988 ; TAJIMA 1989A ; SLATKIN and HUDSON 1991 ; MARJORAM and DONNELLY 1994 ). For example, NORDBORG 2000 proposes a model that could characterize Neanderthal admixture. Consider two populations that diverged at time T₁ in the past and remain completely isolated until time T₀, at which point they merge into a single panmictic population in proportions c and 1 - c (see Fig 1). The genealogy in Fig 1 differs from a panmictic one in that the amount of time when there are two ancestors is much larger than its expectation under panmixia. NORDBORG 2000 shows that under some extreme admixture scenarios, panmixia can be rejected because this time is too long.

View larger version (15K): In this window In a new window Download PPT slide   Figure 1. A model of population subdivision. The population is panmictic from time 0 to time T0 (going backward in time). At time T0, lineages are separated into two isolated subpopulations, in proportions c and 1 - c. Then, at T1 and further back in time, the remaining lineages are randomly mating. Throughout, a coalescent model with recombination is employed (cf. HUDSON 1990 ).

Nordborg notes that his results are not conservative (i.e., the power to reject panmixia with actual data is less than in his simulations). This is because he assumes that the genealogy and branch lengths are known. In practice, they must be inferred from observed sequence data. This inference is especially difficult when there is intragenic recombination. I deal with the problem of intragenic recombination by examining patterns that are directly observable in sequence data. The model assumes that there are data from many unlinked neutral loci. Unlinked loci essentially provide independent runs of the evolutionary process, because demographic history is expected to equally affect the whole genome. Another advantage of considering multiple unlinked loci is that we can minimize the confounding effect of local selection.

I introduce ad hoc test statistics based on measures of inferred linkage disequilibrium or inferred genealogical tree shape, and I examine their powers to reject the null model (of panmixia) when the actual demography is as in Fig 1. More effective measures likely exist. However, my goal was to construct reasonable upper bounds for how much information (i.e., how many loci) would be necessary to distinguish between various demographic hypotheses. Full-likelihood methods were not considered. Algorithms that include intragenic recombination have only been implemented for panmictic models (e.g., GRIFFITHS and MARJORAM 1996 ), and even these are not computationally feasible for the recombination rates considered here (results not shown).

An outline of the remainder of the article is as follows. First, the model and test statistics are introduced. Then, the test statistics are compared with each other, and the best one is chosen for further simulations. These simulations are parameterized to model Neanderthal-modern human admixture in Europe ~30,000–45,000 years ago (e.g., MERCIER et al. 1991 ; HUBLIN et al. 1995 , HUBLIN et al. 1996 ) and more speculative admixture between modern humans and H. erectus remnants in Southeast Asia (ROBERTS et al. 1994 ; SWISHER et al. 1996 ). Finally, shortcomings of the model are discussed, as are general issues regarding demographic inference from current sequence data.

	MODEL AND TEST STATISTICS

TOP ABSTRACT MODEL AND TEST STATISTICS RESULTS DISCUSSION LITERATURE CITED

To model archaic admixture, I use an infinite-sites coalescent model with recombination (cf. HUDSON 1983 , HUDSON 1990 ). I assume there is no selection and a constant rate of recombination per base pair per generation. Consider a panmictic population with diploid effective population size N. Going backward in time, suppose that at time T₀ the remaining ancestors are placed randomly into one of two subpopulations with probabilities c and 1 - c. From time T₀ to time T₁, these subpopulations are assumed to be completely isolated from each other and panmictic with diploid effective population sizes of cN and N, respectively. [Replacing cN by N or N by (1 - c)N did not change any of the results.] Then, at time T₁, all remaining lineages are placed into a single panmictic population with diploid effective population size N. All further coalescences occur in this panmictic population. The null model is c = 0, which yields constant-sized panmixia throughout time (i.e., no contribution of archaic populations to the modern gene pool). T₀ and T₁ are scaled in units of 2N generations. This model is similar to that proposed by NORDBORG 2000 , except that it includes intragenic recombination. One possible problem with the model is that the mixing event at time T₀ is assumed to be instantaneous, whereas possible interbreeding between modern and archaic human populations might have been spread over thousands of years. However, the instantaneous approximation is reasonable because few coalescent events are expected to occur during the time of contact and because only small values of c are considered.

The data simulated using the above scheme are in the standard population genetic format of an ordered list of segregating sites. Call a pair of segregating sites congruent if the subset of the data consisting of the two sites contains only two different haplotypes. Congruence is a pairwise measure of linkage disequilibrium; it has been shown to be effective in detecting certain types of population structure (WALL 1999 ). If T₁ >> T₀, then at time T₀ there will be many sites that are fixed for alternative alleles in the two subpopulations. This will show up in data as a set of sites that is dispersed throughout a locus and is pairwise congruent. (A set of sites is pairwise congruent if all possible pairs in the set are congruent.) Over time these associations will decay, but they serve as the motivation for the following test statistics.

Label the segregating sites 1, 2, ... , S in order, and suppose 1 i j S. Define the indicator variable

Consider the subset of the data consisting of the sites i, i + 1, ... , j. If this subset has evidence of recombination (by containing a pair of sites with | D' | < 1; cf. HUDSON and KAPLAN 1985 ), then set I_r(i, j) = 1. Otherwise, define I_r(i, j) = 0. Now, let

g is the number of pairs of sites that are congruent, g_r is the number of pairs of sites that are both congruent and have evidence of recombination between them, and g_d is the maximal distance (measured in segregating sites) between any two congruent sites. For g_d, a distance measure based on the number of intervening base pairs should be no more effective; those areas containing Neanderthal ancestry will also have on average higher levels of polymorphism (due to longer genealogies) and thus longer distances when measured by the number of intervening segregating sites.

NORDBORG 2000 uses the total time when there are exactly two ancestors as a test statistic. This is similar to using the longest branch length in the unrooted tree as a test statistic, because the longest branch is often the oldest one as well. (This correspondence is less meaningful when recombination causes different sites to have different genealogies.) As mentioned before, tree topologies and branch lengths must be estimated from the data. One measure with similar motivation is

l_b is the maximum number of segregating sites that are all congruent to each other. If there were no recombination, l_b would be an estimate of the number of mutations on the longest branch. In practice, it is still useful even in the presence of intragenic recombination (see RESULTS). Another possible test statistic is the maximal number of nucleotide differences between two sampled individuals, denoted by k_max. This statistic has been used in estimating the time to the most recent common ancestor (T_MRCA) of a sample (e.g., TEMPLETON 1993 ; FU and LI 1997 ). The T_MRCA is expected to be substantially longer in scenarios such as depicted in Fig 1.

One final test statistic considered is TAJIMA's (1989b) D (denoted here by d_t), a commonly used measure of the skew in the observed frequency spectrum. When c is small, mutations on the oldest branch will tend to be at low frequency (assuming there is no outgroup), leading to a shift toward negative d_t values. Because the multilocus likelihood calculations described below require discrete variables, each d_t value was rounded down to the nearest tenth.

Suppose an alternative model Alt = (c, T₀, T₁, , C) is specified. (Here = 4Nµ is the population mutation parameter and C = 4Nr is the population recombination parameter.) It is compared with a null model Null = (c = 0, T₀, T₁, , C). The 's for the two models are chosen so that the average pairwise divergences (; cf. TAJIMA 1983 ) for the simulated data sets have roughly equal expectation. Furthermore, we assume that = C for each model. (This assumption is relaxed in the DISCUSSION.) Recent high-resolution linkage maps (e.g., BOUFFARD et al. 1997 ; NAGARAJA et al. 1997 ) and studies of human genetic variation (e.g., HARDING et al. 1997 ; ZIETKIEWICZ et al. 1997 ; NICKERSON et al. 1998 ; HARRIS and HEY 1999 ) suggest that this is a reasonable approximation. Note that the actual values of T₀ and T₁ are irrelevant for the null model.

We can estimate the distribution of the seven test statistics by simulation. A statistic will be effective when its distribution under the alternative model of archaic admixture is noticeably different from its distribution under the null model of equilibrium panmixia. For any value of k, set Pr_Null (g = k) as the probability that g = k if the null model is simulated. (Pr_Alt is defined similarly.) This probability can be estimated from the proportion of simulated trials where g = k. In general, 10⁶ trials were considered sufficient for each combination of parameter values. Now, suppose that data is available from x independent loci, each simulated according to either the null model or the alternative model. Define g(j) = the value of g from the jth locus, and

G is akin to a likelihood ratio and is used as a one-tailed test statistic. An empirical distribution for G is determined from 10⁵ simulations of the null model. Then, the power to reject the null model (at the 5% level) is determined from 10⁵ simulations under the alternative model. G values that are too low are rejected. G_r, G_q, etc., are all defined analogously.

	RESULTS

TOP ABSTRACT MODEL AND TEST STATISTICS RESULTS DISCUSSION LITERATURE CITED

Preliminary simulations were run to estimate the optimal sample size, locus size, and test statistic. It was assumed that there is a constant per base pair cost for sequencing; under the constraint of a fixed total cost, it was determined what combination of n, , number of loci, and test statistic led to maximal power to reject the null model. These simulations were meant to be exploratory, not exhaustive. They suggest that the optimal configuration has a relatively small sample size (n ~ 20) and a relatively large locus size ( ~ 10), and that power is maximized using G_r (results not shown). All further simulations fixed = 10 and n = 50. This corresponds to roughly 8–10 kb of sequence. The larger than optimal sample size was taken to conform more closely to the large sample sizes considered in most recent studies of human genetic variation (e.g., HARDING et al. 1997 ; ZIETKIEWICZ et al. 1997 ; NICKERSON et al. 1998 ). As an illustration, Fig 2 shows the power (in percentage) vs. the number of loci for all six test statistics. The other parameters are c = 0.1, T₀ = 0.1, and T₁ = 1.0. If we make the standard assumptions of N = 10⁴ and a generation time of 20 years, then T₀ is 40,000 years ago and T₁ is 400,000 years ago. G_r, L_b, and G are all reasonably powerful with ~40 or more loci, while the others remain weak even when a large number of loci are considered. G_r is slightly more powerful than G because it weights those congruent pairs that are "far" (as measured by evidence of recombination) over those that are "near." The former are much more unusual under the panmictic equilibrium null model.

View larger version (19K): In this window In a new window Download PPT slide   Figure 2. Power vs. the number of loci for all six test statistics. Here C = = 10, n = 50, c = 0.1, T0 = 0.1, and T1 = 1.0. () Gr, () Lb, () G, (X) Gd, (*) Kmax, (•) Dt.

The most powerful test statistic (G_r) was chosen for more extensive simulations that sought to determine how much information (i.e., how many loci) would be required to distinguish between plausible human demographic scenarios. The results are described first for possible Neanderthal admixture in Europe and second for possible H. erectus admixture in Southeast Asia.

Recent archaeological evidence suggests that modern humans first entered Europe roughly 45,000 years ago, and that remnant Neanderthal populations remained until ~25,000–30,000 years ago (e.g., MERCIER et al. 1991 ; HUBLIN et al. 1995 , HUBLIN et al. 1996 ; LAHR 1996 ; DUARTE et al. 1999 ). Making the same approximations of N = 10⁴ and a generation time of 20 years yields the estimate T₀ ~ 0.07–0.11. It is more difficult to make a reasonable guess for T₁. This is due to both the multiple interpretations of the fossil record and the fact that genetic isolation might be only loosely correlated with morphological divergence. FOLEY and LAHR 1997 suggest that the divergence may be relatively recent (i.e., roughly 250,000 years ago), while others have postulated that it is far older (e.g., CARBONELL et al. 1995 ; HUBLIN 1998 ). I simulated T₁ values from 0.6 to 2.0 (increment 0.2), which cover the whole range of possibilities. Also, T₀ values from 0.06 to 0.2 (increment 0.02) were considered. Results are shown in Fig 3 and Fig 4. Fig 3 has c = 0.1 while Fig 4 has c = 0.02. For both, the first graph shows the power to reject the null model when there are 60 loci, and the second shows the power when 250 loci are considered. Not surprisingly, more extreme subdivision (i.e., smaller T₀ values or larger T₁ values) leads to greater power. For the genealogy of a given site to look like Fig 1, at least one lineage must be present in both subpopulations. For the parameters considered here (i.e., c << 1 and many expected ancestors at time T₀), representation in the nonarchaic subpopulation is certain. So, a site's genealogy will appear as in Fig 1 if and only if it contains archaic ancestry. Higher values of c lead to a greater proportion of sites having archaic ancestry and hence greater power. These same qualitative results were noted by NORDBORG 2000 . The probability that a given site has archaic ancestry depends on both c and the number of ancestors present at time T₀. As T₀ increases, there are fewer expected ancestors at time T₀ and thus a smaller probability of a particular site having archaic ancestry and less power to reject the null model. Alternatively, as T₁ increases, those sites with archaic ancestry will have a longer oldest branch; this will lead to more segregating sites in strong linkage disequilibrium with each other and greater power to reject the null model. Fig 3A shows that 60 loci are sufficient to distinguish between high levels of Neanderthal admixture and panmixia in most situations, whereas Fig 4A suggests that the power to detect lower (possibly more realistic) levels of admixture depends much more on specific parameter assumptions. In particular, if Neanderthals were genetically isolated only recently (cf. FOLEY and LAHR 1997 ), then high power to reject the null model might require on the order of hundreds of loci.

View larger version (30K): In this window In a new window Download PPT slide   Figure 3. Power as a function of T0 and T1. The other parameters are C = = 10, n = 50, and c = 0.1. (a) Power with 60 independent loci; (b) power with 250 loci.

View larger version (40K): In this window In a new window Download PPT slide   Figure 4. Power as a function of T0 and T1. The other parameters are C = = 10, n = 50, and c = 0.02. (a) Power with 60 independent loci; (b) power with 250 loci.

With the recent claim that H. erectus survived on Java until just 25,000 years ago (SWISHER et al. 1996 ), there has arisen the intriguing possibility that H. erectus and modern humans might have been contemporaneous in Southeast Asia. Modern humans first colonized Australia ~60,000 years ago (ROBERTS et al. 1994 ), and due to geographical constraints they must have at least temporarily occupied Southeast Asia. This time of possible contact corresponds to T₀ values of 0.06–0.15 or more. H. erectus first colonized East Asia almost 2 mya (e.g., SWISHER et al. 1994 ; HUANG et al. 1995 ), and the marked morphological conservatism in the Southeast Asian H. erectus fossil record is suggestive of a high degree of genetic isolation. The 1.6–1.8 mya dates of SWISHER et al. 1994 suggest an upper bound for T₁ of ~4.0–4.5. Fig 5 shows the results of power simulations where T₀ varies from 0.06 to 0.20 (increment 0.02) and T₁ varies from 1.0 to 5.0 (increment 0.5). To reflect the knowledge that any contribution of H. erectus to the modern gene pool is likely to be minor, a low value of c (c = 0.005) was assumed. Fig 5A shows power with 60 loci, while Fig 5B shows power with 250 loci. If H. erectus were isolated for most of its existence in Southeast Asia (e.g., T₁ ~ 3 or more), then there would be reasonable power with 60 loci and very high power with 250 loci even when T₀ is large; those loci with archaic ancestry would appear extremely unusual, and this would more than counteract the low value of c. However, if there had been intermittent gene flow then lower values of T₁ would be appropriate; if so, there would be reduced power, because the vast majority of sites with no archaic ancestry would swamp out the pattern from other sites.

View larger version (44K): In this window In a new window Download PPT slide   Figure 5. Power as a function of T0 and T1. The other parameters are C = = 10, n = 50, and c = 0.005. (a) Power with 60 independent loci; (b) power with 250 loci.

	DISCUSSION

TOP ABSTRACT MODEL AND TEST STATISTICS RESULTS DISCUSSION LITERATURE CITED

Fig 3 Fig 4 Fig 5 show that the number of loci needed to have reasonable power to reject the null model is quite sensitive to assumptions about when and how long the archaic population was isolated and how much of a contribution it made to the modern gene pool. So, while as few as a dozen loci might be sufficient for extreme demographic scenarios, it seems that 60 or more loci would be necessary if parameter values are close to our best guesses for them. Furthermore, some demographic scenarios will be impossible to distinguish from equilibrium panmixia even if the entire genome is sequenced in every living human. If the contact between archaic and modern populations happened too far in the past (i.e., if T₀ is too large) and/or the archaic contribution was too small (i.e., if c is too small), then genetic drift in the intervening years could easily obscure whatever pattern there once was. One conclusion is that genetic sequence data can never exclusively support a single origin model with complete replacement, because very low levels of hybridization with the local archaic populations (thus very few sites with archaic ancestry) can never be ruled out. Another is that there is not enough human sequence information available at the present for us to make any real demographic inferences. As of July 1999, there are only four published nuclear data sets of reasonable size (i.e., n > 10 and > 5) for which haplotype information is available (Lpl, NICKERSON et al. 1998 ; Pdha1, HARRIS and HEY 1999 ; Xq13.3, KAESSMANN et al. 1999 ; ACE, RIEDER et al. 1999 ; I exclude MHC loci, which are known to be under strong diversifying selection). Sufficient data will probably be collected in the next 5–10 years, and until then it is probably most appropriate to take an agnostic attitude toward the debate on whether archaic populations have contributed to the modern human gene pool.

Even though very low values of c (e.g., c < 0.001) might be impossible to distinguish from c = 0 by the methods described here, they could still produce patterns in the data that would be highly suggestive of population structure. Lower c values imply a higher proportion of sites that conform to an equilibrium panmictic model because they have no archaic ancestors. However, conditional on a given site containing archaic ancestry, the expected length of sequence sharing this ancestry is quite large. For example, when T₀ = 0.1 and n = 50, the expected length is roughly the distance over which C ~ 30 (i.e., ~20–30 kb; results not shown). Also, the expected frequency of archaic variants (conditional on their existence) is roughly one over the expected number of ancestors at time T₀ (~7% when T₀ = 0.1), so is not vanishingly small (results not shown). If T₁ is large (e.g., T₁ > 3 in the case of Southeast Asian H. erectus ancestry), then areas of the genome with archaic ancestry would contain an unusually large number of reasonably rare segregating sites in high linkage disequilibrium with each other. These sites would not be tightly bunched, but would be dispersed over several kilobases of sequence.

I have focused on the power to reject c = 0 when the actual situation is c > 0. It would be just as easy to use multilocus data to come up with a point estimate of c (i.e., an estimate of what proportion of our genome came from archaic human populations). Suppose M(c) = (c, T₀, T₁, , C) is a class of models where all but the first parameter are fixed, and g_r(j) is the value of g_r from the j^th locus. Then one could estimate c by taking the value of c that maximizes ^l_j=1Pr_M(c) (g_r = g_r (j)). This method of using maximum likelihood on summary statistics has been used in other situations (e.g., FU and LI 1997 ; WEISS and VON HAESELER 1998 ). It will be informative to examine the likelihood curve for c when more real data are available.

The model considered in this article is somewhat simplistic; it assumes that modern humans (and most of their direct ancestors) form a constant-sized panmictic population, that the archaic population is formed and dissolved instantaneously, and that the archaic population is completely isolated for the duration of its existence. Some of these assumptions are known to be false. Below I discuss some possible shortcomings of the model assumptions and how more realistic models would change the general conclusions reached here.

Complete isolation of the archaic population:
It has been suggested that multiple dispersal events might be a common pattern in human evolution (LAHR and FOLEY 1994 ; LAHR 1996 ). Also, the multiregional model posits low levels of gene flow between all hominid populations (e.g., WOLPOFF et al. 1984 ). One way of incorporating this into the model is to include a (presumably low) migration rate between the two subpopulations. Because migration has a homogenizing effect, increasing the migration rate leads to a less extreme subdivision model and less power to reject the null model. Scenarios that include a reasonable amount of migration are difficult to distinguish from the null model because they are in actuality very similar to panmictic models.

Constant population size:
The recent population explosion highlights the fact that human population sizes have changed over time. Some researchers have postulated that human population sizes were small in the distant past and then started to grow rapidly during the upper Paleolithic (e.g., ROGERS 1995 ; STINER et al. 1999 ). This has been modeled by constant size followed by exponential growth (e.g., MARJORAM and DONNELLY 1994 ). One aspect of exponentially growing populations is that coalescent events tend to be bunched around the time of the onset of growth. If the onset of growth precedes T₀, then a sample will tend to have more ancestors present at time T₀ than it would in the constant-sized case, which leads to a greater proportion of sites with archaic ancestry and greater power to reject the null model. If, on the other hand, T₀ is further back in time than the exponential growth phase, then power depends on the scaled value of T₀, which in turn depends on assumptions about N. The power in this situation is very close to the power for the comparable constant-sized population case. Thus, the assumption of a constant population size is conservative with respect to recent population growth.

Panmixia:
The human population is not panmictic. Genetic studies have consistently found differences between regional populations (e.g., VIGILANT et al. 1991 ; CAVALLI-SFORZA et al. 1994 ; HARDING et al. 1997 ; HARRIS and HEY 1999 ). In general, any form of geographic subdivision tends to produce an excess of linkage disequilibrium between alleles and could lead to spurious rejection of the null model (i.e., the null assumption of panmixia in the modern human population is not conservative with respect to underlying population structure). This problem might be avoided if a suitable model of geographic subdivision could be used as the null model. For example, an island model could be considered with migration rates determined from observed F_st values between populations. Under an island model, the archaic contribution to the modern gene pool would likely be geographically localized. The actual population structure in humans is probably quite complex, and it remains to be shown whether simple models (e.g., island models) are appropriate substitutes.

Uncertainty in C:
This uncertainty can take two forms. First, the ratio C/ may be greater than one. As this ratio increases, there will be less linkage disequilibrium and less correlation in the genealogies of nearby segregating sites and thus less useful information for making demographic inferences. However, as long as C is known (so we can condition on how much linkage disequilibrium is expected), there is no loss of power. For example, if Fig 3A is rerun with C = 20 (instead of C = 10) everywhere, the power to reject the null model is essentially unchanged (results not shown).

A more serious concern is that regional variation in recombination rates might cause C to be different for each locus. If so, the test statistics described here could be used only if C could be measured accurately for each locus. This would require linkage maps with ~10-fold finer resolution than BOUFFARD et al. 1997 and NAGARAJA et al. 1997 , because current methods of estimating C from extant sequence data (e.g., HUDSON 1987 ; GRIFFITHS and MARJORAM 1996 ; HEY and WAKELEY 1997 ) are all biased and have large variances (results not shown). One possible solution would be to choose the most conservative recombination rate for each locus [i.e., for locus j the one that maximizes Pr_Null ]. Another option would be to construct a test statistic whose definition does not depend on an assumed recombination rate. The loss in power by either method is substantial. An example of the latter type of test statistic is

Here k is an arbitrary constant. Note, however, that determining the rejection power of this statistic still requires simulations over many possible values of C. If c = 0.1, T₀ = 0.1, T₁ = 1.0, and k = 7, then the power to reject the null model with f_r(k) and 60 loci is ~30%. In Fig 3A, the comparable power is 94.6%.

Nature of recombination:
The model used here includes crossing over but not gene conversion. The incorporation of gene conversion might lead to higher values of g_r; if a conversion tract is surrounded on both sides by regions with the same genealogy, then it is more likely that one observes two congruent sites with evidence of recombination between them. To avoid this problem, either (1) the null model could be modified to include intragenic recombination or (2) test statistics of similar power (e.g., l_b or g; see Fig 2) could be used that are likely to be less sensitive to assumptions regarding recombination.

Even though there are not yet enough genetic sequence data to answer conclusively questions surrounding modern human origins, it is important to start thinking about how multilocus data can be used to infer demographic history. Unlike selection, which tends to be localized, demography is expected to affect the whole genome. As methods become more sophisticated, we will have the power to consider more complicated (and hopefully more realistic) models of demographic history. The model of ancient structure with recent mixing presented here might be a reasonable approximation of human demographic history. If so, estimating c will give us some idea of what proportion of our genome came from Neanderthals or other archaic human populations. The results obtained here suggest that at least in some cases, tens of thousands of years of random mating are not sufficient to obliterate the patterns of population structure from the distant past.

	ACKNOWLEDGMENTS

I thank R. Hudson and M. Nordborg for helpful discussions and R. Hudson, M. S. McPeek, M. Przeworski, and two anonymous reviewers for comments on an earlier version of this manuscript. This work was motivated by discussions at the North Atlantic Treaty Organization Advanced Study Institute on Human Evolution at the Newton Institute in Cambridge.

Manuscript received April 18, 1999; Accepted for publication October 14, 1999.

	LITERATURE CITED

TOP ABSTRACT MODEL AND TEST STATISTICS RESULTS DISCUSSION LITERATURE CITED

ARMOUR, J. A. L., T. ANTTINEN, C. A. MAY, E. E. VEGA, and A. SAJANTILA et al., 1996  Minisatellite diversity supports a recent African origin for modern humans. Nat. Genet. 13:154-160[Medline].

BOUFFARD, G. G., J. R. IDOL, V. V. BRADEN, L. M. IYER, and A. F. CUNNINGHAM et al., 1997  A physical map of human chromosome 7: an integrated YAC contig map with average STS spacing of 79 kb. Genome Res. 7:673-692[Abstract/Free Full Text].

BRÄUER, G., 1984  The `Afro-European sapiens hypothesis' and hominid evolution in East Asia during the late Middle and Upper Pleistocene. Cour. Forschungsinst. Senckenb. 69:145-165.

CARBONELL, E., J. M. BERMUDEZ DE CASTRO, J. L. ARSUAGA, J. C. DIEZ, and A. ROSAS et al., 1995  Lower Pleistocene hominids and artifacts from Atapuerca-TD6 (Spain). Science 269:826-830[Abstract/Free Full Text].

CAVALLI-SFORZA, L. L., P. MENOZZI and A. PIAZZA, 1994 The History and Geography of Human Genes. Princeton University Press, Princeton, NJ.

DUARTE, C., J. MAURICIO, P. B. PETTITT, P. SOUTO, and E. TRINKAUS et al., 1999  The early Upper Paleolithic human skeleton from the Abrigo do Lagar Velho (Portugal) and modern human emergence in Iberia. Proc. Natl. Acad. Sci. USA 96:7604-7609[Abstract/Free Full Text].

FOLEY, R. and M. M. LAHR, 1997  Mode 3 technologies and the evolution of modern humans. Camb. Arch. J. 7:3-36.

FU, Y.-X. and W.-H. LI, 1997  Estimating the age of the common ancestor of a sample of DNA sequences. Mol. Biol. Evol. 14:195-199[Abstract].

GRIFFITHS, R. C. and P. MARJORAM, 1996  Ancestral inference from samples of DNA sequences with recombination. J. Comp. Biol. 3:479-502.

HARDING, R. M., S. M. FULLERTON, R. C. GRIFFITHS, J. BOND, and M. J. COX et al., 1997  Archaic African and Asian lineages in the genetic ancestry of modern humans. Am. J. Hum. Genet. 60:772-789[Medline].

HARRIS, E. E. and J. HEY, 1999  X chromosome evidence for ancient human histories. Proc. Natl. Acad. Sci. USA 96:3320-3324[Abstract/Free Full Text].

HEY, J. and J. WAKELEY, 1997  A coalescent estimator of the population recombination rate. Genetics 145:833-846[Abstract].

HOWELLS, W. W., 1976  Explaining modern man: evolutionists versus migrationists. J. Hum. Evol. 5:477-495.

HUANG, W., R. CIOCHON, Y. GU, R. LARICK, and Q. FANG et al., 1995  Early Homo and associated artefacts from Asia. Nature 378:275-278[Medline].

HUBLIN, J. J., 1998 Climatic changes, paleogeography, and the evolution of the neandertals, pp. 295–310 in Neandertals and Modern Humans in Western Asia, edited by T. AKAZAWA, K. AOKI and O. BAR-YOSEF. Plenum Press, New York.

HUBLIN, J. J., C. BARROSO RUIZ, P. MEDINA LARA, M. FONTUGNE, and J. L. REYSS, 1995  The Mousterian site of Zafarraya (Andalucia, Spain): dating and implications on the palaeolithic peopling processes of Western Europe. C. R. Acad. Sci. Paris 321:931-937.

HUBLIN, J. J., F. SPOOR, M. BRAUN, F. ZONNEVELD, and S. CONDEMI, 1996  A late Neanderthal associated with upper Paleolithic artefacts. Nature 381:224-226.

HUDSON, R. R., 1983  Properties of a neutral allele model with intragenic recombination. Theor. Popul. Biol. 23:183-201[Medline].

HUDSON, R. R., 1987  Estimating the recombination parameter of a finite population model without selection. Genet. Res. 50:245-250[Medline].

HUDSON, R. R., 1990 Gene genealogies and the coalescent process, pp. 1–44 in Oxford Surveys in Evolutionary Biology, Vol. 7, edited by P. H. HARVEY and L. PARTRIDGE. Oxford University Press, New York.

HUDSON, R. R. and N. L. KAPLAN, 1985  Statistical properties of the number of recombination events in the history of a sample of DNA sequences. Genetics 111:147-164[Abstract/Free Full Text].

KAESSMANN, H., F. HEI&SZLIG;IG, A. VON HAESELER, and S. PÄÄBO, 1999  DNA sequence variation in a non-coding region of low recombination on the human X chromosome. Nat. Genet. 22:78-81[Medline].

KAHN, P. and A. GIBBONS, 1997  DNA from an extinct human. Science 277:176-178[Free Full Text].

KRINGS, M., A. STONE, R. W. SCHMITZ, H. KRAINITZKI, and M. STONEKING et al., 1997  Neanderthal DNA sequences and the origin of modern humans. Cell 90:19-30[Medline].

LAHR, M. M., 1994  The multiregional model of modern human origins: a reassessment of its morphological basis. J. Hum. Evol. 26:23-56.

LAHR, M. M., 1996 The Evolution of Modern Human Diversity. Cambridge University Press, Cambridge, United Kingdom.

LAHR, M. M. and R. FOLEY, 1994  Multiple dispersals and modern human origins. Evol. Anthropol. 3:48-60.

LINDAHL, T., 1997  Facts and artifacts of ancient DNA. Cell 90:1-3[Medline].

MARJORAM, P. and P. DONNELLY, 1994  Pairwise comparisons of mitochondrial DNA sequences in subdivided populations and implications for early human evolution. Genetics 136:673-683[Abstract].

MERCIER, N., H. VALLADAS, J. L. JORON, J. L. REYSS, and F. LEVEQUE et al., 1991  Thermoluminescence dating of the late Neanderthal remains from Saint-Cesaire. Nature 351:737-739[Medline].

NAGARAJA, R., S. MACMILLAN, J. KERE, C. JONES, and S. GRIFFIN et al., 1997  X chromosome map at 75-kb STS resolution, revealing extremes of recombination and GC content. Genome Res. 7:210-222[Abstract/Free Full Text].

NICKERSON, D. A., S. L. TAYLOR, K. M. WEISS, A. G. CLARK, and R. G. HUTCHINSON et al., 1998  DNA sequence diversity in a 9.7-kb region of the human lipoprotein lipase gene. Nat. Genet. 19:233-240[Medline].

NORDBORG, M., 1998  On the probability of Neanderthal ancestry. Am. J. Hum. Genet. 63:1237-1240[Medline].

NORDBORG, M., 2000 On detecting ancient admixture, in Proceedings of the NATO ASI Workshop, "Genes, Fossils, and Behaviour: An Integrated Approach to Human Evolution" (NATO Science Series: Life Sciences, Vol. 310). IOS Press, Amsterdam.

RIEDER, M. J., S. L. TAYLOR, A. G. CLARK, and D. A. NICKERSON, 1999  Sequence variation in the human angiotensin converting enzyme. Nat. Genet. 22:59-62[Medline].

ROBERTS, R. G., R. JONES, N. A. SPOONER, M. J. HEAD, and A. S. MURRAY et al., 1994  The human colonisation of Australia: optical dates of 53,000 and 60,000 years bracket human arrival at Deaf Adder Gorge, Northern Territory. Q. Sci. Rev. 13:575-583.

ROGERS, A. R., 1995  Genetic evidence for a Pleistocene population explosion. Evolution 49:608-615.

SLATKIN, M. and R. R. HUDSON, 1991  Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics 129:555-562[Abstract].

SMITH, F. H., 1992 The role of continuity in modern human origins, pp. 145–156 in Continuity or Replacement? Controversies in Homo sapiens Evolution, edited by G. BRÄUER and F. H. SMITH. A. A. Balkema, Rotterdam, The Netherlands.

STINER, M. C., N. D. MUNRO, T. A. SUROVELL, E. TCHERNOV, and O. BAR-YOSEF, 1999  Paleolithic population growth pulses evidenced by small animal exploitation. Science 283:190-194[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[Abstract/Free Full Text].

SWISHER, C. C., G. H. CURTIS, T. JACOB, A. G. GETTY, and A. SUPRIJO et al., 1994  Age of the earliest known hominids in Java, Indonesia. Science 263:1118-1121[Abstract/Free Full Text].

SWISHER, C. C., W. J. RINK, S. C. ANTON, H. P. SCHWARCZ, and G. H. CURTIS et al., 1996  Latest Homo erectus of Java: potential contemporaneity with Homo sapiens in Southeast Asia. Science 274:1870-1874[Abstract/Free Full Text].

TAJIMA, F., 1983  Evolutionary relationships of DNA sequences in finite populations. Genetics 105:437-460[Abstract/Free Full Text].

TAJIMA, F., 1989a  The effect of change in population size on DNA polymorphism. Genetics 123:597-601[Abstract/Free Full Text].

TAJIMA, F., 1989b  Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123:585-595[Abstract/Free Full Text].

TAKAHATA, N., 1988  The coalescent in two partially isolated diffusion populations. Genet. Res. 52:213-222[Medline].

TEMPLETON, A. R., 1993  The "Eve" hypothesis: a genetic critique and reanalysis. Am. Anthropol. 95:51-72.

TISHKOFF, S. A., E. DIETZSCH, W. SPEED, A. J. PAKSTIS, and J. R. KIDD et al., 1996  Global patterns of linkage disequilibrium at the CD4 locus and modern human origins. Science 271:1380-1387[Abstract].

VIGILANT, L., M. STONEKING, H. HARPENDING, K. HAWKES, and A. C. WILSON, 1991  African populations and the evolution of human mitochondrial DNA. Science 253:1503-1507[Abstract/Free Full Text].

WALL, J., 1999  Recombination and the power of statistical tests of neutrality. Genet. Res. 74:65-79.

WARD, R. and C. STRINGER, 1997  A molecular handle on the Neanderthals. Nature 388:225-226[Medline].

WEIDENREICH, F., 1943  The skull of Sinanthropus pekinensis: a comparative study of a primitive hominid skull. Palaeontol. Sin. Ser. D 10:1-485.

WEISS, G. and A. VON HAESELER, 1998  Inference of population history using a likelihood approach. Genetics 149:1539-1546[Abstract/Free Full Text].

WOLPOFF, M. H., 1996  Interpretations of multiregional evolution. Science 274:704-706[Free Full Text].

WOLPOFF, M. H., X. WU and A. G. THORNE, 1984 Modern Homo sapiens origins; a general theory of hominid evolution involving the fossil evidence from East Asia, pp. 411–483 in The Origins of Modern Humans: A World Survey of the Fossil Evidence, edited by F. SMITH and F. SPENCER. Liss, New York.

ZIETKIEWICZ, E., V. YOTOVA, M. JARNIK, M. KORAB-LASKOWSKA, and K. K. KIDD et al., 1997  Nuclear DNA diversity in worldwide distributed human populations. Gene 205:161-171[Medline].

This article has been cited by other articles:

				K. E. Lohmueller, C. D. Bustamante, and A. G. Clark Methods for Human Demographic Inference Using Haplotype Patterns From Genomewide Single-Nucleotide Polymorphism Data Genetics, May 1, 2009; 182(1): 217 - 231. [Abstract] [Full Text] [PDF]

				G. Wlasiuk, S. Khan, W. M. Switzer, and M. W. Nachman A History of Recurrent Positive Selection at the Toll-Like Receptor 5 in Primates Mol. Biol. Evol., April 1, 2009; 26(4): 937 - 949. [Abstract] [Full Text] [PDF]

				M. Fumagalli, R. Cagliani, U. Pozzoli, S. Riva, G. P. Comi, G. Menozzi, N. Bresolin, and M. Sironi Widespread balancing selection and pathogen-driven selection at blood group antigen genes Genome Res., February 1, 2009; 19(2): 199 - 212. [Abstract] [Full Text] [PDF]

				M. Slatkin and J. L. Pollack Subdivision in an Ancestral Species Creates Asymmetry in Gene Trees Mol. Biol. Evol., October 1, 2008; 25(10): 2241 - 2246. [Abstract] [Full Text] [PDF]

				M. P. Cox, F. L. Mendez, T. M. Karafet, M. M. Pilkington, S. B. Kingan, G. Destro-Bisol, B. I. Strassmann, and M. F. Hammer Testing for Archaic Hominin Admixture on the X Chromosome: Model Likelihoods for the Modern Human RRM2P4 Region From Summaries of Genealogical Topology Under the Structured Coalescent Genetics, January 1, 2008; 178(1): 427 - 437. [Abstract] [Full Text] [PDF]

				C. Becquet and M. Przeworski A new approach to estimate parameters of speciation models with application to apes Genome Res., October 1, 2007; 17(10): 1505 - 1519. [Abstract] [Full Text] [PDF]

				P. D. Evans, N. Mekel-Bobrov, E. J. Vallender, R. R. Hudson, and B. T. Lahn Evidence that the adaptive allele of the brain size gene microcephalin introgressed into Homo sapiens from an archaic Homo lineage PNAS, November 28, 2006; 103(48): 18178 - 18183. [Abstract] [Full Text] [PDF]

				A. Geraldes, N. Ferrand, and M. W. Nachman Contrasting Patterns of Introgression at X-Linked Loci Across the Hybrid Zone Between Subspecies of the European Rabbit (Oryctolagus cuniculus) Genetics, June 1, 2006; 173(2): 919 - 933. [Abstract] [Full Text] [PDF]

				D. Garrigan, Z. Mobasher, S. B. Kingan, J. A. Wilder, and M. F. Hammer Deep Haplotype Divergence and Long-Range Linkage Disequilibrium at Xp21.1 Provide Evidence That Humans Descend From a Structured Ancestral Population Genetics, August 1, 2005; 170(4): 1849 - 1856. [Abstract] [Full Text] [PDF]

				K. J. Schmid, S. Ramos-Onsins, H. Ringys-Beckstein, B. Weisshaar, and T. Mitchell-Olds A Multilocus Sequence Survey in Arabidopsis thaliana Reveals a Genome-Wide Departure From a Neutral Model of DNA Sequence Polymorphism Genetics, March 1, 2005; 169(3): 1601 - 1615. [Abstract] [Full Text] [PDF]

				D. Garrigan, Z. Mobasher, T. Severson, J. A. Wilder, and M. F. Hammer Evidence for Archaic Asian Ancestry on the Human X Chromosome Mol. Biol. Evol., February 1, 2005; 22(2): 189 - 192. [Abstract] [Full Text] [PDF]

				E. Trinkaus, O. Moldovan, s. Milota, A. Bilgar, L. Sarcina, S. Athreya, S. E. Bailey, R. Rodrigo, G. Mircea, T. Higham, et al. An early modern human from the Pestera cu Oase, Romania PNAS, September 30, 2003; 100(20): 11231 - 11236. [Abstract] [Full Text] [PDF]

				J. L. Mountain, A. Knight, M. Jobin, C. Gignoux, A. Miller, A. A. Lin, and P. A. Underhill SNPSTRs: Empirically Derived, Rapidly Typed, Autosomal Haplotypes for Inference of Population History and Mutational Processes Genome Res., November 1, 2002; 12(11): 1766 - 1772. [Abstract] [Full Text] [PDF]

				D. J. Mathews, C. Kashuk, G. Brightwell, E. E. Eichler, and A. Chakravarti Sequence Variation Within the Fragile X Locus Genome Res., August 1, 2001; 11(8): 1382 - 1391. [Abstract] [Full Text] [PDF]

				Y. Koda, H. Tachida, H. Pang, Y. Liu, M. Soejima, A. A. Ghaderi, O. Takenaka, and H. Kimura Contrasting Patterns of Polymorphisms at the ABO-Secretor Gene (FUT2) and Plasma {{alpha}}(1,3)Fucosyltransferase Gene (FUT6) in Human Populations Genetics, June 1, 2001; 158(2): 747 - 756. [Abstract] [Full Text] [PDF]