Keywords: population genetics, admixture, single nucleotide polymorphism (SNP) array

The three-population test and introduction off-statistics

We begin with a description of the three-population test. First we give some theory. Consider the tree ofFigure 1A. We see that the path fromC toA and the path fromC toB just share the edge fromC toX. Leta′,b′,c′ be allele frequencies in the populationsA,B,C, respectively, at a single polymorphism. Define

We, similarly, in an obvious notation define

Choice of the allele does not affect any ofF₂,F₃,F₄ as choosing the alternate allele simply flips the sign of both terms in the product. We refer toF₂(A,B) as thebranch length between populationsA andB. We use these branch lengths in admixture graph fitting for graph edges.

OurF-values should be viewed as population parameters, but we note that they depend both on the demography and choice of SNPs. In Appendix A we give formulae that use sample frequencies and that yield unbiased estimates of the correspondingF parameters. The unbiased estimates ofF computed using these formulae at each marker are then averaged over many markers to form ourf-statistics.

The results that follow hold rigorously if we identify the polymorphisms we are studying in an outgroup (that is, we select SNPs based on patterns of genetic variation in populations that all have the same genetic relationship to populationsA,B,C). Since only markers with variation inA,B,C are relevant to the analysis, then by ascertaining in an outgroup we ensure that our markers are polymorphic in the root population ofA,B,C. Later on, we discuss how other strategies for ascertaining polymorphisms would be expected to affect our results. In general, our tests for admixture and estimates of admixture proportion are strikingly robust to the ascertainment processes that are typical for human SNP array data, as we verify both by simulations and by empirical analysis.

Suppose the allele frequency of a SNP isr at the root. In the tree ofFigure 1A, leta′,b′,c′,x′,r′ be allele frequencies inA,B,C,X,R. Condition onr′. Then

sinceE[a′|x′] =x′, andE[x′ −b′] =E[r′ −b′ − (r′ −x′)] = 0. If the phylogeny hasC as an outgroup (switchingB,C inFigure 1A), then a similar argument shows that

There is an intuitive way to think about the expected values off-statistics, which relies on tracing the overlap of genetic drift paths between the first and second terms in the quadratic expression, as illustrated inFigure 2 and discussed further in Appendix B. For example,E[(c′ −a′)(c′ −b′)] can be negative only if populationC has ancestry from populations related to bothA andB. Only in this case are there paths betweenC andA andC andB that also take opposite drift directions through the tree (Figure 1C andFigure 2), which contributes to a negative expectation for the statistics. The observation of a significantly negative value off₃(C;A,B) is thus evidence of complex phylogeny inC. We prove this formally in Appendix C (Theorem 1). In Appendix D, we also relax our assumptions about the ascertainment process, showing thatF₃ is guaranteed to be positive ifC is unadmixed under quite general conditions, for example, polymorphic in the rootR and in addition ascertained as polymorphic in any ofA,B,C. It is important to recognize, however, that a history of admixture does not always result in a negativef₃(C;A,B)-statistic. If populationC has experienced a high degree of population-specific drift (perhaps due to founder events after admixture), it can mask the signal so thatf₃(C;A,B) might not be negative.

An important feature of this test is that it definitively shows that the history of mixture occurred in populationC; a complex history forA orB cannot produce negativeF₃(C;A,B). To explain why this is so, we recapitulate material fromReichet al. (2009). If populationA is admixed then if we pick an allele ofA, it must have originated in one of the admixing populations. Pick allelesα,β from populationsA andB andγ₁,γ₂ independently fromC, coding 1 for a reference allele, 0 for a variant, etc. Thus,F₃(C;A,B) =E[(γ₁ −α)](γ₂ −β)]. Suppose populationA is admixed;B andC are not admixed. The alleleα sampled from populationA can take more than one path through the ancestral populations.F₃(C;A,B) can then be computed as a weighted average over the possible phylogenies, in all of which the quantity has a positive expectation becauseA andB are now unadmixed (Appendix B andFigure 2). In conclusion, the diagram makes it visually evident that ifF₃(C;A,B) < 0 then populationC itself must have a complex history.

Additivity ofF₂ along a tree branch

In this article we consider generalizations of phylogenetic trees where graph edges indicate that one population is a descendant of another. Consider the phylogenetic tree inFigure 1B, and a marker polymorphic at the root. Drift on a given edge is a random variable with mean 0. For ifA →B is a graph edge, with corresponding allele frequenciesa′,b′,

This is themartingale property of allele frequency diffusion. Drifts on two distinct edges of a tree are orthogonal, where orthogonality of random variablesX,Y simply means thatE[XY] = 0. In our context this means that the drifts on distinct edges have mean0 and are uncorrelated.

A valuable feature of ourF-statistics definition is that branch lengths on the tree (as defined byF₂) are additive. We illustrate this with an example from human history (Figure 1B). (We note that all examples in this article refer to human history, although the methods should apply equally well to other species.) In this example,A, andC are present-day populations that split from an ancestral populationX.B is an ancestral population toC. For instance,A might be modern Yoruba,C a European population, andB an ancient population, perhaps a sample from archeological material of a population that existed thousands of years ago. We assume here that we ascertain in an outgroup (implying polymorphism at the root) and again assume neutrality and that we can ignore recurrent or backmutations. Then we mean by additivity that

for

but the last term is 0 since the change in allele frequencies (drifts)X →A,X →B,B →C are all uncorrelated.

We remark that ourF₂-distance resembles the familiarF_st, but is not the same. In particular, parts of a graph that are far from the root (in genetic drift distance) haveF₂ reduced. Some insight into this effect is given by considering the simple graph

whereτ₁,τ₂ are drift times on the standard diffusion timescale (two random alleles ofB have probability$e^{-{\tau }_2}$ that they have not coalesced in the ancestral populationA).

Ifr′,a′,b′ are allele frequencies inR,A,B, respectively, thenF₂(A,B) =E[(a′ −b′)²]. WriteE_r′,E_a′ for expectations conditional on population allele frequenciesr′,a′. Then$E_{a\prime }\left[{\left(a\prime -b\prime \right)}^2\right]=a\prime \left(1-a\prime \right)\left(1-e^{-{\tau }_2}\right)$ (Nei 1987, Chap. 13). Moreover$E_{r\prime }\left[a\prime \left(1-a\prime \right)\right]=r\prime \left(1-r\prime \right)e^{-{\tau }_1}$. Hence

Informally the drift fromR →A shrinksF₂(A,B) by a factor$e^{-{\tau }_1}$. Thus expected drift isadditive,

but the drift does depend on ascertainment. For a given edge, the more distant the root, the smaller the drift. A loose analogy is projecting a curved surface, such as part of the globe, into a plane. Locally all is well, but any projection will cause distortion in the large. Additivity inF₂ distances is all we require in what follows. We note that there is no assumption here that population sizes are constant along a branch edge, and so we arenot assuming linearity of branch lengths in time.

Expected values of ourf-statistics

We can calculate expected values for ourf-statistics, at least for simple demographic histories that involve population splits and admixture events. We assume that genetic drift events on distinct edges are uncorrelated, which as mentioned before will be true if we ascertain in an outgroup, and our alleles are neutral.

We give an illustration forf₃-statistics. Consider the demography shown inFigure 1C. PopulationsE,F split from a root populationR.G then was formed by admixture in proportionsα:β (β = 1 −α). Modern populationsA,B,C are then formed by drift fromE,F,G. We want to calculate the expected value off₃(C;A,B). Assume that our ascertainment is such that drifts on distinct edges are orthogonal, which will hold true if we ascertained the markers in an outgroup.

We recapitulate some material from (Reichet al. 2009, Supplementary S2, Sect. 2.2). As before leta′,b′,c′ be population allele frequencies inA,B,C, and letg′ be the allele frequency inG and so on:

We see by orthogonality of drifts that

which we write as

Now, label alleles at a marker 0, 1. Then picking chromosomes from our populations independently we can write

wherea₁,b₁ are alleles chosen randomly in populationsA,B andg₁,g₂ are alleles chosen randomly and independently in populationG. Similarly, we definee₁,e₂,f₁, andf₂. However,g₁ originated fromE with probabilityα and so on. Thus

wherea₁,a₂ are independently picked fromE andb₁,b₂ fromF. The first three terms vanish. Further

This shows that under our assumptions of orthogonal drift on distinct edges,

It might appear thatFigure 1C is too restricted, as it assumes that the admixing populationsE,F are ancestral toA,B and that we should consider the more general graph shown inFigure 1D. But it turns out that using ourf-statistics alone (and not the more general allelic spectrum) that even ifα,β are known, we can obtain information only about

Thus in fitting admixture graphs tof-statistics, we can, without loss of generality, fit all the genetic drift specific to the admixed population on the lineage directly ancestral to the admixed population (the lineage leading fromC toG inFigure 1C).

The outgroup case

Care though is needed in interpretation. ConsiderFigure 1E. Here a similar calculation to the one just given shows (again assuming orthogonality of drift on each edge) that

Note thatY has little to do with the admixture intoC and we obtain the sameF₃ value forany populationY that splits off fromA more anciently thanX.

We call this case, where we have apparent admixture betweenA andY, theoutgroup case, and it needs to be carefully considered when recovering population relationships.

Estimates of mixing proportions

We want to estimate, or at least bound, the mixing proportions that have resulted in the ancestral population ofC. With further strong assumptions on the phylogeny we can get quite precise estimates even without accurate surrogates for the ancestral populations (seeReichet al. 2009 and theF₄-ratio estimation that we describe below, for examples). Also if we have data from populations that are accurate surrogates for the ancestral admixing population (and we can ignore the drift post admixture), the problem is much easier. For instance, inPattersonet al. (2010) we give an estimator that works well even when the sample sizes of the relevant populations are small, and we have multiple admixing populations whose deep phylogenetic relationships we may not understand. Here we show a method that obtains useful bounds, without requiring full knowledge of the phylogeny, although the bounds are not very precise. Note that although our three-population test remains valid even if the populationsA,B are admixed, the mixing proportions we calculate are not meaningful unless the assumed phylogeny is at least roughly correct. Indeed even discussing mixing from an ancestral population ofA hardly makes sense ifA is admixed itself subsequent to the admixing event inC. This is discussed further when we present data from Human Genome Diversity Panel (HGDP) populations.

In much of the work in this article, we analyze some populationsA,B,C and need an outgroup, which split off from the ancestral population ofA,B,C before the population split ofA,B. For example, inFigure 1E,Y is such an outgroup. Usually, when studying a group of populations within a species, a plausible outgroup can be proposed. The outgroup assumption can then be checked using the methods of this article, by adding an individual from a more distantly related population, which can be treated as a second outgroup. For instance, with human populations from Eurasia, Yoruba or San Bushmen from sub-Saharan Africa will often be plausible outgroups.¹ Our second outgroup here is simply being used to check a phylogenetic assumption in our primary analysis, and we donot require polymorphism at the root for this narrow purpose. Chimpanzee is always a good second outgroup for studies of humans.

Consider the phylogeny ofFigure 1F. Hereα,β are mixing parameters (α +β = 1) and we show drift distances along the graph edges. Note that here we usea,b,…, as branch lengths (F₂ distances), not sample or population allele frequencies as we do elsewhere in this article. Thus, for example,F₂(O,X) =u. Now we can obtain estimates of

We also have estimates of

SetY_i =Z_i −Z₀,i = 0…5, which eliminatesu. This shows that any populationO which is a true outgroup should (up to statistical noise) give similar estimates forY_i (Figure 1F). We have three inequalities:

Usingαa =Y₁,βb =Y₂ we can rewrite these as

giving lower and upper bounds onα, which we write asα_L,α_U in the tables of results that follow. These bounds can be computed by a programqpBound in the ADMIXTOOLS software package that we make available with this article.

Although these bounds will be nearly invariant to choices of the outgroupO, choices for the source populationsA,B may make a substantial difference. We give an example in a discussion of the relationship of Siberian populations to Europeans. In principle we can give standard errors for the bounds, but these are not easily interpretable, and we think that in most cases systematic errors (for instance, that our phylogeny is not exactly correct) are likely to dominate.

We observe that in some cases the lower bound exceeds the upper, even when theZ-score for admixture of populationC is highly significant. We interpret this as suggesting that our simple model for the relationships of the three populations is wrong. A negativeZ-score indeed implies thatC has a complex history, but ifA orB also has a complex history, then a recovered mixing coefficientα has no real meaning.

Estimation and normalization

With all ourf-statistics it is critical that we can compute unbiased estimates of the populationF-parameter for a single SNP, with finite sample sizes. Without that, our estimates will be biased, even if we average over many unlinked SNPs. The explicit formulae forf₂,f₃,f₄ that we present in Appendix A (previously given inReichet al. 2009, Supplementary Material) are in fact minimum variance unbiased estimates of the correspondingF-parameters, at least for a single marker.

The expected (absolute) values of anf-statistic, such asf₃, strongly depends on the distribution of the derived allele frequencies of the SNPs examined; for example, if many SNPs are present that have a low average allele frequency across the populations being examined, then the magnitude off₃ will be reduced. To see this, suppose that we are computingf₃(C;A,B), and as beforea′,b′,c′ are population frequencies of an allele inA,B,C. If the allele frequencies are small, then it is obvious that the expected value off₃(C;A,B) will be small in absolute magnitude as well. Importantly, however, the sign of anf-statistic is not dependent on the absolute magnitudes of the allele frequencies (all that it depends on is the relative magnitudes across the populations being compared). Thus, a significant deviation of anf-statistic from 0 can serve as a statistically valid test for admixture, regardless of the ascertainment of the SNPs that are analyzed. However, to reduce the dependence of the value of thef₃-statistic on allele frequencies for some of our practical computations, in all of the empirical analyses we report below, we normalize using an estimate for each SNP of the heterozygosity of the target populationC. Specifically, for each SNPi, we compute unbiased estimates${\widehat{T}}_i$,${\widehat{B}}_i$ of both

Now we normalize ourf₃-statistic computing

This greatly reduces the numerical dependence off₃ on the allelic spectrum of the SNPs examined, without making much difference to statistical significance measures such as aZ-score. We note that we usef₃ and$f_3^{\mathrm{\star }}$ interchangeably in many places in this article. Both of these statistics give qualitatively similar results and thus if the goal is only to test iff₃ has negative expected value then the inference should be unaffected.

D-statistics

TheD-statistic test was first introduced inGreenet al. (2010) where it was used to evaluate formally whether modern humans have some Neandertal ancestry. Further theory and applications ofD-statistics can be found inReichet al. (2010) andDurandet al. (2011). A very similar statisticf₄ was used to provide evidence of admixture in India (Reichet al. 2009), where we called it a four-population test. TheD-statistic was also recently used as a convenient statistic for studying locus-specific introgression of genetic material controlling coloration in Heliconius butterflies (Dasmahapatraet al. 2012).

LetW,X,Y,Z be four populations, with a phylogeny that corresponds to the unrooted tree ofFigure 3A. For SNPi suppose variant population allele frequencies arew′,x′,y′,z′, respectively. Choose an allele at random from each of the four populations. Then we define a “BABA” event to mean that theW andY alleles agree, and theX andZ alleles agree, while theW andX alleles are distinct. We define an “ABBA” event similarly, now with theW andZ alleles in agreement. LetNum_i andDen_i be the numerator and denominator of the statistic:

For SNP data these values can be computed using either population or sample allele frequencies.Durandet al. (2011) showed that replacing population allele frequencies (w′,y′, etc.) by the sample allele frequencies yields unbiased estimates ofNum_i,Den_i. Thus ifw,x,y,z are sample allele frequencies we define

and, in a similar spirit to our normalizedf₃-statistic$f_3^{\mathrm{\star }}$ we define theD-statisticD(W,X;Y,Z) as

summing both the numerator and denominator over many SNPs and only then taking the ratio. If we ascertain in an outgroup, then if (W,X) and (Y,Z) are clades in the population tree, it is easy to see thatE[Num_i] = 0. We can compute a standard error forD using the weighted block jackknife (Businget al. 1999). The number of standard errors that this quantity is from zero forms aZ-score, which is approximately normally distributed and thus yields a formal test for whether (W,X) indeed forms a clade.

More generally, if the relationship of the analyzed populations is as shown inFigure 3B orFigure 3C and we ascertain in an outgroup or in {W,X} thenD should be zero up to statistical noise. The reason is that ifU is the ancestral population toY,Z andu′,y′,z′ are population allele frequencies inU,Y,Z, thenE[y′ −z′|u′] =E[y′|u′] −E[z′|u′] = 0. Here there is no need to assume polymorphism at the root of the tree, as for a SNP to make a nonzero contribution toD we must have polymorphism at both {Y,Z} and {W,X}. If the tree assumption is correct, drift betweenY,Z and betweenW,X are independent so thatE[Num_i] = 0. Thus testing whetherD is consistent with zero constitutes a test for whether (W,X) and (Y,Z) are clades in the population tree.

As mentioned earlier,D-statistics are very similar to the four-population test statistics introduced inReichet al. (2009). The primary difference is in the computation of the denominator ofD. For statistical estimation, and testing for “treeness,” theD-statistics are preferable, as the denominator ofD, the total number of ABBA and BABA events, is uninformative for whether a tree phylogeny is supported by the data, whileD has a natural interpretation: the extent of the deviation on a normalized scale from −1 to 1.

As an example, let us assume that two human Eurasian populationsA,B are a clade with respect to West Africans (Yoruba). Assume the phylogeny shown inFigure 3D and that we ascertain in an outgroup toA,B. Then

F₄-ratio estimation

F₄-ratio estimation, previously referred to asf₄-ancestry estimation inReichet al. (2009), is a method for estimating ancestry proportions in an admixed population, under the assumption that we have a correct historical model.

Consider the phylogeny ofFigure 4. The populationX is an admixture of populationsB′ andC′ (possibly with subsequent drift). We have genetic data from populationsA,B,X,C,O.

SinceF₄(A,O;C′,C) = 0 it follows that

Thus an estimate ofα is obtained as

where the estimates in both numerator and denominator are obtained by summing over many SNPs.

As we can obtain unbiasedF₄-statistics by sampling a single allele from each population, we can apply this test to sequence data, where we pick a single allele, from a high-quality read, for all relevant populations at each polymorphic site. In practice this must be done with care as both sequencing error that is correlated between samples and systematic misalignment of reads to a reference sequence can distort the statistics.

Examples ofF₄-ratio estimation

Reichet al. (2009) provide evidence that most human South Asian populations can be modeled as a mixture of Ancestral North Indians (ANI) and Ancestral South Indians (ASI) and that if we set, using the labeling above,

• Label Population

• A Adygei

• B CEU (HapMap European Americans)

• X Indian (many populations)

• C Onge (indigenous Andamanese)

• O Papuan (Dai and HapMap YRI West Africans also work)

we get estimates of the mixing coefficients that are robust, have quite small standard errors, and are in conformity with other estimation methods. See Reichet al. (2009, Supplementary S5) for further details.

As another example, inReichet al. (2010) andGreenet al. (2010) evidence was given that there was gene flow (introgression) from Neandertals into non-Africans. Further, a sister group to Neandertals, “Denisovans” represented by a fossil from Denisova cave, Siberia, shows no evidence of having contributed genes to present-day humans in mainland Eurasia (Reichet al. 2010,2011). The phylogeny is that ofFigure 4 if we set:

• Label Population

• A Denisova

• B Neandertal

• X French (or almost any population from Eurasia)

• C Yoruba

• O Chimpanzee

HereB′ is the population of Neandertals that admixed, which forms a clade with the Neandertals from Vindija that were sequenced byGreenet al. (2010). So for this example, we obtain an estimate ofα, the proportion of Neandertal gene flow into French as 0.022 ± 0.007 (seeReichet al. 2010, SI8, for more detail).

Simulations to test the accuracy off- andD-statistic-based historical inferences

We carried out coalescent simulations of five populations related according toFigure 4, usingms (Hudson 2002). Detailed information about the simulations is given in Appendix D.

Table 1 shows that using the three-population test,D-statistics, andF₄-ratio estimation, we reliably detect mixture events and obtain accurate estimates of mixture proportions, even for widely varied demographic histories and strategies for discovering polymorphisms.

The simulations also document important features of our methods. As mentioned earlier, the only case where thef₃-statistic for a population that is truly admixed fails to be negative is when the population has experienced a high degree of population-specific genetic drift after the admixture occurred. Further, theD-statistics show a substantial deviation from 0 only when an admixture event occurred in the history of the four populations contributing to the statistic. Finally, the estimates of admixture proportions usingF₄-ratio estimation are accurate for all ascertainment strategies and demographies.

Effect of ascertainment process onf- andD-statistics

So far, we have assumed that we have sequence data from all populations and ascertainment is not an issue. However, the ascertainment of polymorphisms (for example, enriching the set of analyzed SNPs for ancestry informative markers) can modulate the magnitudes ofF₃,F₄, andD. Empirically, we observe that in commercial SNP arrays developed for genome-wide association studies (like Affymetrix 6.0 and Illumina 610-Quad), ascertainment does indeed affect the observed magnitudes of these statistics, but importantly, does not cause them to be biased away from zero if this is their expected value in the absence of complex ascertainment (e.g., for complete genome sequencing data). This is key to the robustness of our tests for admixture: since our tests are largely based on evaluating whether particularf- orD-statistics are consistent with zero, and SNP ascertainment almost never causes a deviation from zero, the ascertainment process does not appear to be contributing to spuriously significant signals of admixture. We have verified this through two lines of analysis. First, we carried out simulations showing that tests of admixture (as well asF₄-ratio estimation) perfomed using these methods are robust to very different SNP ascertainment strategies (Table 1). Second, we report analyses of data from a new SNP array with known ascertainment that we designed specifically for studies of population history. Even when we use radically different ascertainment schemes, and even when we use widely used commercial SNP arrays, inferences about history are indistinguishable (Table 9).

Interpretation and limitations ofqpGraph

• 1.

  A major use ofqpGraph is to show that a hypothesized phylogeny must be incorrect. This generalizes ourD-statistic test, which is testing a simple tree on four populations.

• 2.

  After fitting parameters, study of whichf-statistics fit poorly can lead to insights as to how the model must be wrong.

• 3.

  Overfitting can be a problem, especially if we hypothesize many admixing events, but only have data for a few populations.

Simulations validate the performance ofqpGraph

We show inFigure 5 an example in which we simulated a demography with five observed populationsOut,A,B,C, andX and one admixture event. We simulated 50,000 unlinked SNPs, ascertained as heterozygous in a single diploid individual from the outgroupOut. Sample sizes were 50 in all populations and the historical population sizes were all taken to be 10,000. We show that we can accurately recover the drift lengths and admixture proportions usingqpGraph.

Choice of weight function

In many applications, we have access to two modern populationsA,B, which we can regard as surrogates for the true admixing populations, and in this context we can simply use the difference of empirical frequencies of the variant allele as our weight. For example, to study the admixture in African Americans, very good surrogates for the ancestral populations are Yoruba and North Europeans. However, a strength of rolloff is that it provides unbiased dates even without access to accurate surrogates for the ancestral populations. That is, rolloff is robust to use of highly divergent populations as surrogates. In cases when the ancestrals are no longer extant or data from the ancestrals are not available, but we have access to multiple admixed populations with differing admixture proportions (as for instance happens in India (Reichet al. 2009), we can use the “SNP loadings” generated from principal component analysis (PCA) as appropriate weights. This also gives unbiased dates for the admixture events.

Simulations to test rolloff

We ran three sets of simulations. The goals of these simulations were

• 1.

  to access the accuracy of the estimated dates, in cases for which data from accurate ancestral populations are not available,

• 2.

  to investigate the bias seen inMoorjaniet al. (2011),

• 3.

  to test the effect of genetic drift that occurred after admixture.

We describe the results of each of these investigations in turn.

• 1.

  First, we report simulation results that test the robustness of inferences of dates of admixture when data from accurate ancestral populations are not available. We simulated data for 20 individuals using phased data from HapMap European Americans (CEU) and HapMap West Africans (YRI), where the mixture date was set to 100 generations before present and the proportion of European ancestry was 20%. We ran rolloff using pairs of reference populations that were increasingly divergent from the true ancestral populations used in the simulation. The results are shown inTable 2 and are better than those of the rather similar simulations inMoorjaniet al. (2011). Here we use more SNPs (378K instead of 83K) and 20 admixed individuals rather than 10. The improved results likely reflect the fact that we are analyzing larger numbers of admixed individuals and SNPs in these simulations, which improves the accuracy of rolloff inferences by reducing sampling noise in the calculation of theZ-score. In analyzing real data, we have found that the accuracy of rolloff results improves rapidly with sample size; this feature of rolloff contrasts markedly with allele frequency correlation statistics likef-statistics where the accuracy of estimation increases only marginally as sample sizes increase above five individuals per population.

• 2.

  Second, we report simulation results investigating the bias seen inMoorjaniet al. (2011).Moorjaniet al. (2011) showed that low sample size and admixture proportion can cause a bias in the estimated dates. In our new simulations, we generated haplotypes for 100 individuals using phased data from HapMap CEU and HapMap YRI, where the mixture date was between 50 and 800 generations ago (Figure 6) and the proportion of European ancestry was 20%. We ran rolloff with two sets of reference populations: (1) the true ancestral populations (CEU and YRI) and (2) the divergent populations Gujarati (F_st(CEU, Gujarati) = 0.03 and Maasai (F_st(YRI, Maasai) = 0.03). We show the results for one run and the mean date from each group of 10 runs inFigures 6, A and B. These results show no important bias, and the date estimates, even in the more difficult case where we used Gujarati and Maasai as assumed ancestrals, are tightly clustered near the “truth” up to 500 generations (around 15,000 years). This shows that the bias is removed with larger sample sizes.

• 3.

  The simulations reported above sample haplotypes without replacement, effectively removing the impact of genetic drift after admixture. To study the effect of drift post-dating admixture, we performed simulations using the MaCS coalescent simulator (Chenet al. 2009). We simulated data for one chromosome (100 Mb) for three populations (say,A,B, andC). We set the effective population size (N_e) for all populations to 12,500, the mutation rate to 2 × 10⁻⁸/bp/generation, and the recombination rate to 1.0 × 10⁻⁸/bp/generation. Consider the phylogeny inFigure 1C.G is an admixed population that has 80%/20% ancestry fromE andF, with an admixture time (t) set to be 30, 100, or 200 generations before the present. PopulationsA,B,C are formed by drift fromE,F,G, respectively.F_st(A,B) = 0.16 (similar to that ofF_st(YRI, CEU)). We performed rolloff analysis withC as the target (n = 30) andA andB as the reference populations. We estimated the standard error using a weighted block jackknife where the block size was set to 10 cM. The estimated dates of admixture were 28 ± 4, 97 ± 10, and 212 ± 19 corresponding the true admixture dates of 30, 100, and 200 generations, respectively. This shows that the estimated dates are not measurably affected by genetic drift post-dating the admixture event.

A SNP array designed for population genetics

We conclude our presentation of our methods by describing a new experimental resource and publicly available data set that we have generated for facilitating studies of human population history and that we use in many of the applications that follow.

For studies that aim to fit models of human history to genetic data, it is highly desirable to have an exact record of how polymorphisms were chosen. Unfortunately, conventional SNP arrays developed for medical genetics have a complex ascertainment process that is nearly impossible to reconstruct and model (but seeWollsteinet al. 2010). While the methods reported in our study are robust in theory and also in simulation to a range of strategies for how polymorphisms were ascertained (Table 1), we nevertheless wished to empirically validate our findings on a data set without such uncertainties.

Here, we report on a novel SNP array that we developed that is now released as theAffymetrix Human Origins array. This includes 13 panels of SNPs, each ascertained in a rigorously documented way that is described inFile S1, allowing users to choose the one most useful for a particular analysis. The first 12 are based on a strategy used inKeinanet al. (2007), discovering SNPs as heterozygotes in a single individual of known ancestry for whom sequence data are available (fromGreenet al. 2010;Reichet al. 2010) and then confirming the site as heterozygous with a different assay. After the validation steps described inFile S1 (which serves as technical documentation for the new SNP array), we had the following number of SNPs from each panel: San, 163,313; Yoruba, 124,115; French, 111,970; Han, 78,253; Papuan (two panels), 48,531 and 12,117; Cambodian, 16,987; Bougainville, 14,988; Sardinian, 12,922; Mbuti, 12,162; Mongolian, 10,757; Karitiana, 2,634. The 13th ascertainment consisted of 151,435 SNPs where a randomly chosen San allele was derived (different from the reference Chimpanzee allele) and a randomly chosen Denisova allele (Reichet al. 2010) was ancestral (same as chimpanzee). The array was designed so that all sites from panels 1–13 had data from chimpanzee as well as from Vindija Neandertals and Denisova, but the values of the Neandertal and Denisova alleles were not used for ascertainment (except for the 13th ascertainment).

Throughout the design process, we avoided sources of bias that could cause inferences to be affected by genetic data from human samples other than the discovery individual. Our identification of candidate SNPs was carried out entirely using sequencing reads mapped to the chimpanzee genome (PanTro2), so that we were not biased by the ancestry of the human reference sequence. In addition, we designed assays blinded to prior information on the positions of polymorphisms and did not take advantage of prior work that Affymetrix had done to optimize assays for SNPs already reported in databases. After initial testing of 1,353,671 SNPs on two screening arrays, we filtered to a final set of 542,399 SNPs that passed all quality-control criteria. We also added a set of 84,044 “compatibility SNPs” that were chosen to have a high overlap with SNPs previously included on standard Affymetrix and Illumina arrays, to facilitate coanalysis with data collected on other SNP arrays. The final array contains 629,443 unique and validated SNPs, and its technical details are described inFile S1.

We successfully genotyped the array in 934 samples from the HGDP and made the data publicly available on August 12, 2011, atftp://ftp.cephb.fr/hgdp_supp10/. The present study analyzes a curated version of this data set in which we have used principal component analysis (Pattersonet al. 2006) to remove samples that are outliers relative to others from their same populations; 828 samples remained after this procedure. This curated data set is available for download from the Reich laboratory website (http://genetics.med.harvard.edu/reich/Reich_Lab/Datasets.html).

Initial application to data: South African Xhosa

The Xhosa are a South African population whose ancestors are mostly Bantu speakers from the Nguni group, although they also have some Bushman ancestors (Pattersonet al. 2010). We first ran our three-population test with San (HGDP) (Cannet al. 2002) and Yoruba (HapMap) (International Hapmap 3 Consortium 2010) as source populations and 20 samples of Xhosa as the target population, a sample set already described inPattersonet al. (2010). We obtain anf₃-statistic of −0.009 with aZ-score of −33.5, as computed with the weighted block jackknife (Businget al. 1999).

Note that the admixing Bantu-speaking population is known to have been Nguni and certainly was not Nigerian Yoruba. However, as explained earlier, this is not crucial, if the actual admixing population is related genetically (Bantu speakers have an ancient origin in West Africa). Ifα is the admixing proportion of San here, we obtain using our bounding technique with Han Chinese as an outgroup,

Although this interval is wide, it does show that the Bushmen have made a major contribution to Xhosa genomes.

Xhosa:rolloff

We then applied our rolloff technique, using San and Yoruba as the reference populations, obtaining a very clear exponential admixture LD curve (Figure 7A). We estimate a date of 25.3 ± 1.1 generations, yielding a date of about 740 ± 30 years before present (YBP) assuming 29 years per generation (we also assume this generation time in the analyses that follow) (Fenner 2005).

Archeological and linguistic evidence show that the Nguni are a population that migrated south from the Great Lakes area of East Africa. For the dating of the migration we quote:

From an archaeological perspective, the first appearance of Nguni speakers can be recognized by a break in ceramic style; the Nguni style is quite different from the Early Iron Age sequence in the area. This break is dated to about AD 1200 (Huffman 2010).

More detail on Nguni migrations and archeology can be found inHuffman (2004).

Our date is slightly more recent than the dates obtained from the archeology, but very reasonable, since gene flow from the Bushmen into the Nguni plausibly continued after initial contact.

Admixture of the Uygur

The Uygur are known to be historically admixed, but we wanted to try our methods on them. We analyzed a small sample (nine individuals from HGDP) (Cannet al. 2002). Our three-population test using French and Japanese as sources and Uygur as target gives aZ-score of −76.1, a remarkably significant value. Exploring this a little further, we get the results shown inTable 3.

Using Han instead of Japanese is historically more plausible and statistically not significantly different. Our bounding methods suggest that the West Eurasian admixtureα is in the range

We used French and Han for the source populations here. Russian as a source is significantly weaker than French. We believe that the likely reason is that our Russian samples have more gene flow from East Asia than the French samples, and this weakens the signal. We confirm this by finding thatD(Yoruba, Han; French, Russian) = 0.192,Z = 26.3. The fact that we obtain very similar statistics when we substitute a very different sub-Saharan African population (HGDP San) for Yoruba (D = 0.189,Z = 23.9) indicates that the gene flow does not involve an African population, and instead the findings reflect gene flow between relatives of the Han and Russians.

Uygur:rolloff

Applying rolloff we again get a very clear decay curve (Figure 7B). We estimate a date of 790 ± 60 YBP.

Uygur genetics has been analyzed in two articles by Xu, Jin, and colleagues (Xuet al. 2008;Xu and Jin 2008), using several sets of samples, one of which is the same set of HGDP samples we analyze here. Xu and Jin, primarily using ancestry informative markers (AIMs), estimate West Eurasian admixture proportions of ∼50%, in agreement with our analysis, but also an admixture date estimate using STRUCTURE 2.0 (Falushet al. 2003) of more than 100 generations that is substantially older than ours.

Why are the admixture dates that we obtain so much more recent than those suggested by Xu and Jin? We suspect that STRUCTURE 2.0 systematically overestimates the admixture date, when the reference populations (source populations for the admixture) are not close to the true populations, so that the assumed distribution of haplotypes is in error. It has been suggested (Mackerras 1972) that the West Eurasian component was Tocharian, an ancient Indo-European-speaking population, whose genetics are essentially unknown. Xu and Jin used 60 European American (HapMap CEU) samples to model the European component in the Uygur, and if the admixture is indeed related to the Tocharians it is plausible that they were substantially genetically drifted relative to the CEU, providing a potential explanation for the discrepancy.

Our date of ∼800 years before present is not in conformity withMackerras(1972), who places the admixture in the eighth century of the common era. Our date though is rather precisely in accordance with the rise of the Mongols under Genghis Khan (1206–1368), a turbulent time in the region that the Uygur inhabit. Could there be multiple admixture events and we are primarily dating the most recent?

Northern European gene flow into Spain

While investigating the genetic history of Spain, we discovered an interesting signal of admixture involving Sardinia and northern Europe. We made a data set by merging genotypes from samples from the population reference sample (POPRES) (Nelsonet al. 2008), HGDP (Liet al. 2008), and HapMap Phase 3 (International Hapmap 3 Consortium 2010). We ran our three-population test on triples of populations using Spain as a target (admixed population). We had 137 Spanish individuals in our sample. With Sardinian fixed as a source, we find a clear signal using almost any population from northern Europe.Table 4 gives the topf₃-statistics with correspondingZ-scores. The high score for the Russian and Adygei is likely to be partially confounded with the effect discussed in the section on flow from Asia into Europe (below).

A geographical structure is clear, with the largest magnitudef₃-statistics seen for source populations that are northern European or Slavic. TheZ-score is unsurprisingly more significant for populations with a larger sample size. (Note that positiveZ-scores are not meaningful here.) We were concerned that the Slavic scores might be confounded by a central Asian component and therefore decided to concentrate our attention on Ireland as a surrogate for the ancestral population as they have a substantial sample size (n = 62).

Spain:rolloff

We applied rolloff to Spain using Ireland and Sardinians as the reference populations. InFigure 7C we show a rolloff curve. The rolloff of signed LD out to about 2 cM is clear and gives an admixture age of 3600 ± 400 YBP (the standard error was computed using a block jackknife with a block size of 5 cM).

We have detected here a signal of gene flow from populations related to present-day northern Europeans into Spain around 2000 B.C. We discuss a likely interpretation. At this time there was a characteristic pottery termed “bell-beakers” believed to correspond to a population spread across Iberia and northern Europe. We hypothesize that we are seeing here a genetic signal of the “Bell-Beaker culture” (Harrison 1980). Initial cultural flow of the Bell-Beakers appears to have been from South to North, but the full story may be complex. Indeed one hypothesis is that after an initial expansion from Iberia there was a reverse flow back to Iberia (Czebreszuk 2003); this “reflux” model is broadly concordant with our genetic results, and if this is the correct explanation it suggests that this reverse flow may have been accompanied by substantial population movement. (SeeFigures 8,9, and10.)

It is important to point out that we are not detecting gene flow from Germanic peoples (Suevi, Vandals, Visigoths) into Spain even though it is known that they migrated into Iberia around 500 A.D. We believe such migration must have occurred, based on the historical record (and perhaps is biasing our admixture date to be too recent), but any accompanying gene flow must have occurred at a lower level than the much earlier flow we discuss.

An example of theoutgroup case

Populations closely related geographically often mix genetically, which leaves a clear signal in PCA plots. An example is that isolation-by-distance effects dominate much of the genetic patterning of Europe (Laoet al. 2008;Novembreet al. 2008). This can lead to significantf₃-statistics and is related to theoutgroup case we have already discussed. Here is an example. We find

[YRI are HapMap Yoruba Nigerians (International Hapmap 3 Consortium 2010).] Sub-Saharan populations (including HGDP San) all give aZ < −4.0 when paired with Albania, and evenf₃(Greece; Albania, Papuan) = −0.0033 (Z = −3.5). There may be a low level of sub-Saharan ancestry in our Greek samples, contributing to our signal, but the consistent pattern of highly significantf₃-statistics suggests that we are primarily seeing an outgroup case.

We attempted to date Albanian-related gene flow into Greece using rolloff [with HapMap Yoruba and Albanian as the source populations (Figure 7D)]. However, the technique evidently fails here. Formally we get a date of 62 ± 77 generations, which is not significantly different from zero. It is possible that the admixture is very old (>500 generations) or the gene flow was continuous at a low level, and our basic rolloff model does not work well here.

Admixture events detected in Human Genome Diversity Panel populations

We ran ourf₃-statistic on all possible triples of populations from the HGDP, genotyped on an Illumina 650Y array (Table 5) (Rosenberg 2006;Liet al. 2008).

Here we show for each HGDP target population (column 3) the two-source populations with the most negative (most significant)f₃-statistic. We computeZ using the block jackknife as we did earlier and just show entries withZ < −4. We bound α, the mixing coefficient involving the first source population, as

whereα_L,α_U are computed with HGDP San as outgroup using the methodology of estimating mixing proportions that we have already discussed.

In four cases indicated by an asterisk in the last column,α_L >α_R, suggesting that our three-population phylogeny is not feasible. We suspect (and in some cases the table itself proves) that here the admixing (source) populations are themselves admixed.

It is likely that there are other lines in our table where our source populations are admixed, but that this has not been detected by our rather coarse admixing bounds. In such situations our bounds may be misleading.

Many entries are easily interpretable, for instance the admixture of Uygur (Xuet al. 2008;Xu and Jin 2008) (which we have already discussed), Hazara, Mozabite (Corander and Marttinen 2006;Liet al. 2008), and Maya (Maoet al. 2007) are historically attested. The entry for Bantu-SouthAfrica is likely detecting the same phenomenon that we already discussed in connection with the Xhosa.

However, there is much of additional interest here. Note, for example, the entry for Tu, a people with a complex history and clearly with both East Asian and West Eurasian ancestry. It is important to realize that the finding here by no means implies that the target population is admixed from the two given source populations. For example, in the second line, we do not believe that Japanese, or modern Italians, have contributed genes to the Hazara. Instead one should interpret this line as meaning that an East Asian population related genetically to a population ancestral to the Japanese has admixed with a West Eurasian population. As another example, the most negativef₃-statistic for the Maya arises when we use as source populations Mozabite (north African) and Surui (an indigenous population of South America in whom we have detected no post-Colombian gene flow). The Mozabites are themselves admixed, with sub-Saharan and West Eurasian gene flow. We think that the Maya samples have three-way admixture (European, West African, and Native American) and the incorrect two-way admixture model is simply doing the best it can (Table 5).

Insensitivity to the ascertainment of polymorphisms

In theMaterials and Methods section we described a novel SNP array with known ascertainment that we developed specifically for population genetics (now available as the Affymetrix Human Origins array). The array contains SNPs ascertained in 13 different ways, 12 of which involved ascertaining a heterozygote in a single individual of known ancestry from the HGDP. We genotyped 934 unrelated individuals from the HGDP (Cannet al. 2002) and here report the value off₃-statistics on either SNPs ascertained as a heterozygote in a single HGDP San individual, or at SNPs ascertained in a single Han Chinese (Table 6). We showZ-statistics for these two ascertainments in the last two columns. The number of SNPs used is reduced relative to the 644,247 analyzed inLiet al. (2008); we had 124,440 SNPs for the first ascertainment and 59,251 for the second ascertainment, after removing SNPs at hypermutable CpG dinucleotides. Thus, we expect standard errors onf₃ to be larger and theZ-scores to be smaller, as we observe. The correlation coefficient between theZ-scores for the 2008 data (Z₂₀₀₈) and our newly ascertained data are in each case ∼0.99. We were concerned that this correlation coefficient might be inflated by the very largeZ-statistics for some populations, such as the Hazara and Uygur, but the correlation coefficients remain very large if we divide the table into two halves and analyze separately the most significant and least significant entries.

Ascertainment on aSan heterozygote or aHan heterozygote are very different phylogenetically, and theSan are unlikely to have been used in the construction of the 2008 SNP panel, so the consistency of findings for these distinct ascertainment processes provides empirical evidence, confirming our expectations from theory and findings from simulation (Table 1) that the SNP ascertainment process does not have a substantial effect on inferences of admixture from thef₃-statistics (Table 6).

Evidence for Northeast Asian-related genetic material in Europe

We single out fromTable 5 the score for French arising as an admixture of Karitiana, an indigenous population from Brazil, and Sardinians. TheZ-score of −18.4 is unambiguously statistically significant. We do not of course think that there has been substantial gene flow back into Europe from Amazonia.

The only plausible explanation we can see for our signal of admixture into the French is that an ancient northern Eurasian population contributed genetic material to both the ancestral population of the Americas and the ancestral population of northern Europe. This was quite surprising to us, and in the remainder of the article this is the effect we discuss.

We are not dealing here with theoutgroup case, where the effect is simply caused by Sardinian-related gene flow into the French. If that were the case, then we would expect to see that (French, Sardinian) are approximately a clade with respect to sub-Saharan Africa and Native Americans. There is some modest level of sub-Saharan (probably West African related) gene flow from Africa into Sardinia as is shown by analyses inMoorjaniet al. (2011), but no evidence for gene flow from the San (Bushmen), which is indeed historically most unlikely. But if we computeD(San, Karitiana, French, Sardinian) we obtain a value of −0.0178 and aZ-score of −18.1. Thus we have here gene flow “related” to South America into mainland Europe to a greater extent than into Sardinia.

Further confirmation

We merged two SNP array data sets that included data from Europeans and other relevant populations: POPRES (Nelsonet al. 2008) and HGDP (Liet al. 2008). We considered only populations with a sample size of at least 10.

We considered European populations with Sardinian and Karitiana as sources and computed the statisticf₃(X; Karitian, Sardinian), whereX is various European populations. We also added Druze, as a representative population of the Middle East (Table 7). The effect is pervasive across Europe, with nearly all populations showing a highly significant effect. Orcadians and Cyprus are island populations with known island-specific founder events that could plausibly mask admixture signals produced by the three-population test, so the absence of the signal in these populations does not provide compelling evidence that they are not admixed. Our Cypriot samples are also likely to have some proportion of Levantine ancestry (like the Druze) that does not seem to be affected by whatever historical events are driving our negativef₃-statistic. We can use any Central American or South American population to demonstrate this effect, in place of the Karitiana.

If we replace the Sardinian population by Basque as a source, the effect is systematically smaller, but still enormously statistically significant for most of the populations of Europe (Table 7). We note that in our three populations from mainland Italy [TSI (Hapmap Tuscans), Tuscan, and Italian] the effect essentially disappears when using Basque as a source, although it is quite clear and significant with Sardinian. This is not explored further here, but suggests that further investigation of the genetic relationships of Basque, Sardinian, and other populations of Europe might be fruitful.

Replication using a novel SNP array

The signal above is overwhelmingly statistically significant but we found the effect quite surprising, especially as on common-sense grounds one would expect substantial recent gene flow from the general Spanish and French populations into the Basque, and from mainland Italy into Sardinia, which would weaken the observed effect. We wanted to exclude the possibility that what we are seeing here is an effect of how SNPs were chosen for the medical genetics array used for genotyping. Could the ascertainment be producing false-positive signals of admixture? If, for example, SNPs were chosen specifically so that the population frequencies were very different in Sardinia and northern Europe, an artifactual signal would be expected to arise. This seemed implausible but we had no way to exclude it.

We therefore returned to analysis of data from the Affymetrix Human Origins SNP array with known ascertainment. We show statistics forf₃(French; Karitiana, Saridinian) for all 13 ascertainments and compare them to the statistics for the genotype data from the Illumina 650Y array developed for medical genetics (Liet al. 2008) (Table 8).

All ourZ-scores are highly significant with a very wide range of ascertainments, except for the ascertainment consisting of finding a heterozygote in a Karitiana sample, where the number of SNPs involved is small (thus reducing power). We can safely conclude that the effect is real and that the French have a complex history.

There is evidence that the effect here is substantially stronger in northern than in southern Europe. We confirm this using the statisticD(San, Karitiana; French, Italian), which has aZ-score of −6.4 on the Illumina 650Y SNP array panel and −3.5 on our population genetics panel ascertained with a San heterozygote. These results show that the Karitiana are significantly more closely related to the French than to the Italians. The Italian samples here are from Bergamo, northern Italy. A likely explanation for these findings is discussed below where we apply rolloff to date this admixture event.

As an aside we have repeatedly assumed that backmutations (or recurrent mutations) are not importantly affecting our results. As evidence that this assumption is reasonable, inTable 9 we compute two of our most importantD-statistic-based tests for treeness using a variety of increasingly distant outgroups ranging from modern human outgroups to chimpanzee, gorilla, orangutan, and macaque. Results are entirely consistent across this enormous range of genetic divergence. For example, for the crucial statisticD(Outgroup, Karitiana; Sardinian, French), which demonstrates the signal of Northeast Asian-related admixture in northern Europeans, we find thatZ-scores are consistently positive with high significance whichever outgroup is used. As a second example, when we test if the San are consistent with being an outgroup to two Eurasian populations through the statisticD(Outgroup, San; Sardinian, Han) we detect no significant deviation from zero whichever outgroup is used.

Siberian populations

We obtained Illumina SNP array data fromHancocket al. (2011) from the Naukan and Chukchi, Siberian peoples who live in extreme northeastern Siberia. After merging with the 2008 Illumina 650Y SNP array data on HGDP samples (Liet al. 2008) we obtain thef₃-statistics inTable 10.

We can assume here that we have a common admixture event to explain. Although the statistics for Chukchi are (slightly) weaker than those in the Native Americans, we obtain better bounds on the mixing coefficientα of between 5% and 18%. We caution that if the Sardinians are themselves admixed with Asian ancestry although less so than other Europeans (a scenario we think is historically plausible), then we will have underestimated the Asian-related mixture proportion in Europeans.

We wanted to test if (French, Sardinian) form a clade relative to (Karitiana, Chukchi), which would, for example, be the case if the admixing population to northern Europe had a common ancestor with an ancestor of Karitiana and Chukchi. In our data set,

while this hypothesis predictedD = 0. Thus, we can rule out this alternative hypothesis.

One possible explanation for these findings is that the ancestral Karitiana were closer genetically to the northern Eurasian population that contributed genes to northern Europeans than are the Chukchi. The original migration into the Americas occurred at least 15,000 YBP, so there is ample time for some population inflow into the Chukchi peninsula since then. However, the Chukchi and Naukan samples show no evidence of recent west Eurasian admixture, and we specifically tested for ethnic Russian admixture, finding nothing.

We carried out a rolloff analysis in which we attempted to learn about the date of the admixture events in the history of northern Europeans. We pooled samples from CEU, a population of largely northern European origin (International Hapmap 3 Consortium 2010) with HGDP French to form our target admixed population, wishing to maximize the sample size. The surrogate ancestral populations for this analysis are Karitiana and Sardinian.

The admixture date we are analyzing here is old, and to improve the performance of rolloff here and in the analysis of northern European gene flow into Spain reported above, we filtered out two regions of the genome that have substantial structural variation that is not accurately modeled by rolloff, which assumes Poisson-distributed recombination events between two alleles (Millset al. 2011). The two regions we filtered out were HLA on chromosome 6 and thep-telomeric region on chromosome 8, which we found in practice contributed to anomalous rolloff signals in some of our analyses. Our signals should be robust to removal of small genomic regions.

InFigure 7E we show the rolloff results. The signal is clear enough, although noisy. We estimate an admixture date of 4150 ± 850 YBP. Our standard errors computed using a block jackknife (block size of 5 cM) are uncomfortably large here.

However, this date must be treated with great caution. We obtained a data set from the Illumina iControl database (http://www.illumina.com/science/icontroldb.ilmn) of “Caucasians” and after curation have 1232 samples of European ancestry genotyped on an Illumina SNP array panel. We merged the data with the HGDP Illumina 650Y genotype data obtaining a data set with 561,268 SNPs. Applying rolloff to this sample with HGDP Karitiana and Sardinians as sources, we get a much more recent date of 2200 ± 762 YBP. We think that this is not a technical problem with rolloff, but rather, it is an issue of interpretation that is a challenge for all methods for estimating dates of admixture events.

Our admixture signal is stronger in northern Europe as we showed above in the context of discussing the statisticD(San, Karitiana; French, Italian). It seems plausible that the initial admixture might have been exclusively in northern Europe, but since this ancient event, there has been extensive gene flow within Europe, as shown, for example, inLaoet al. (2008) andNovembreet al. (2008). But if northern and southern Europe have differing amounts of “Asian” admixture, this intra-European flow is confounding to our analysis. The more recent gene flow between northern and southern Europe will contribute to our inferring too recent a date. Admixture into one section of a population, followed by slow mixing within the population, may be quite common in human history and will substantially complicate the dating for any genetic method.

Interpretation in light of ancient DNA

Ancient DNA studies have documented a clean break between the genetic structure of the Mesolithic hunter–gatherers of Europe and the Neolithic first farmers who followed them. Mitochondrial analyses have shown that the first farmers in central Europe, belonging to the linear pottery culture (LBK), were genetically strongly differentiated from European hunter–gatherers (Bramantiet al. 2009), with an affinity to present-day Near Eastern and Anatolian populations (Haaket al. 2010). More recently, new insight has come from analysis of ancient nuclear DNA from three hunter–gatherers and one Neolithic farmer who lived roughly contemporaneously at about 5000 YBP in what is now Sweden (Skoglundet al. 2012). The farmer’s DNA shows a signal of genetic relatedness to Sardinians that is not present in the hunter–gatherers who have much more relatedness to present-day northern Europeans. These findings suggest that the arrival of agriculture in Europe involved massive movements of genes (not just culture) from the Near East to Europe and that people descending from the Near Eastern migrants initially reached as far north as Sweden with little mixing with the hunter–gatherers they encountered. However, the fact that today, northern Europeans have a strong signal of admixture of these two groups, as proven by this study and consistent with the findings of (Skoglundet al. 2012), indicates that these two ancestral groups subsequently mixed.

Combining the ancient DNA evidence with our results, we hypothesize that agriculturalists with genetic ancestry close to modern Sardinians immigrated into all parts of Europe along with the spread of agriculture. In Sardinia, the Basque country, and perhaps other parts of southern Europe they largely replaced the indigenous Mesolithic populations, explaining why we observe no signal of admixture in Sardinians today to the limits of our resolution. In contrast, the migrants did not replace the indigenous populations in northern Europe and instead lived side-by-side with them, admixing over time (perhaps over thousands of years). Such a scenario would explain why northern European populations today are admixed and also have a rolloff admixture date that is substantially more recent than the initial arrival of agriculture in northern Europe.

An alternative history that could produce the signal of Asian-related admixture in northern Europeans is admixture from steppe herders speaking Indo-European languages, who after domesticating the horse would have had a military and technological advantage over agriculturalists (Anthony 2007). However, this hypothesis cannot explain the ancient DNA result that northern Europeans today appear admixed between populations related to Neolithic and Mesolithic Europeans (Skoglundet al. 2012), and so even if the steppe hypothesis has some truth, it can explain only part of the data.

We show an admixture graph that corresponds to our hypothesis inFigure 9.

To test the predictions of our hypothesized historical scenario, we downloaded the recently published DNA sequence of the Tyrolean “Iceman” (Kelleret al. 2012). The Iceman lived (and died) in the Tyrolean Alps close to the border of modern Austria and Italy. From isotopic analysis (Mulleret al. 2003) he was probably born within 60 miles of the site at which he was found. To analyze the Iceman data, we applied similar filtering steps as those applied in the analysis of the Neandertal genome (Greenet al. 2010). After filtering on map quality and sequence quality of a base as described in that study, we chose a random read covering each base of the Affymetrix Human Origins array. This produced nearly 590,000 sites for analysis.

OurD-statistic analysis suggests that the Iceman and the HGDP Sardinians are consistent with being a clade, providing formal support for the findings ofKelleret al. (2012) who reported that the Iceman is close genetically to modern Sardinians based on PCA. Concretely, our test for whether they are a clade is

ThisD-statistic shows no significant deviation from zero, in contrast with the highly significant evidence that the Iceman and French are not a clade:

Our failure to detect a signal of admixture using theD-statistic is not due to reduced power on account of having only one sample, since when we recompute the statistic of (10) using each of the 26 French individuals in turn in place of Iceman, theZ-scores are all significant, ranging from −3.1 to −8.5. These results imply that Iceman has less northeast Asian-related ancestry than a typical modern North European, but the data are consistent with Iceman having the same amount of northeast Asian-related ancestry as Sardinians. Further confirmation for this interpretation comes from the very similar magnitudef₃-statistics that we observe when using either Sardinians or Iceman as a source for the admixture:

TheZ-score for Iceman is of smaller magnitude than that for the Sardinian samples, because with a single individual we have much more sampling noise. However, the important quantity in this context is the magnitude of thef₃-statistic. Thus the Iceman harbors less northeast Asian-related genetic material than modern French, and the northeast Asian-related genetic material is not detectably different in Iceman and the HGDP Sardinians, to the limits of our resolution.

A caveat to these analyses is that the relatively poor quality and highly fragmented DNA sequence fragments from Iceman may occasionally align incorrectly to the reference human genome sequence (and in particular, may do so at a rate higher than that of the comparison data from present-day humans), which could in theory bias theD-statistics. However, our point here is simply that to the limits of the analyses we have been able to carry out, Iceman and modern Sardinians are consistent with forming a clade, supporting the hypothesis we sketched out above.

Although the Iceman lived near where he was found, it cannot be logically excluded that his genetic ancestry was unusual for the region. For instance, his parents might have been migrants from ancient Sardinia. However, the Iceman does not carry the signal of northeast Asian ancestry that we have detected in northern Europeans, and lived at least 2000 years after the arrival of farming in Europe. If his genome was typical of the region in which he lived, the northeast Asian-related genetic material that is currently widespread in northern Italy and southern Austria must be due to admixture events and/or migrations that occurred well after the advent of agriculture in the region, supporting the hypothesis, presented above, that Neolithic farmers of near eastern origin initially largely replaced the indigenous Mesolithic population of southern Europe and that only well afterward did they develop the signal of major admixture that they harbor today.

Summary of inferences about European history from our methods

Our methods for analyzing genetic data have led to several novel inferences about history, showing the power of the approaches. In particular, we have presented evidence suggesting that the genetic history of Europe from around 5000 B.C. includes:

• 1.

  the arrival of Neolithic farmers probably from the Middle East,

• 2.

  nearly complete replacement of the indigenous Mesolithic southern European populations by Neolithic migrants and admixture between the Neolithic farmers and the indigenous Europeans in the north,

• 3.

  substantial population movement into Spain occurring around the same time as the archeologically attested Bell-Beaker phenomenon (Harrison 1980),

• 4.

  subsequent mating between peoples of neighboring regions, resulting in isolation-by-distance (Laoet al. 2008;Novembreet al. 2008). This tended to smooth out population structure that existed 4000 years ago.

Further, the populations of Sardinia and the Basque country today have been substantially less influenced by these events.

Software

We release a software package, ADMIXTOOLS, that implements five methods: the three-population test,D-statistics,F₄-ratio estimation, admixture graph fitting, and rolloff. In addition, it computes lower and upper bounds on admixture proportions based onf₃-statistics. ADMIXTOOLS can be downloaded from the following URL:http://genetics.med.harvard.edu/reich/Reich_Lab/Software.html.

Data sets used

The following data sets are used:

• HapMap Phase 3 (International Hapmap 3 Consortium 2010), HGDP genotyped on the Illumina 650K array (Liet al. 2008), HGDP genotyped on the Affymetrix Human Origins array, POPRES (Nelsonet al. 2008), Siberian data (Hancocket al. 2011), and Xhosa data (Pattersonet al. 2010).