On the limits of fitting complex models of population history to f-statistics

1Admixture graph models provide a powerful intellectual framework for describing the relationships among populations that allows not only branching of populations from a common ancestor but also mixture events. An admixture graph (abbreviated below as AG), as fit in the widely used software packages ADMIXTOOLS and TreeMix, is a directed acyclic bifurcating graph with two types of edges: those representing genetic drift, and those representing gene flow. Each admixture event is represented as a confluence of two gene flow edges. Nodes of such a graph represent unsampled intermediate populations, and terminal nodes (leaves) represent sampled present-day or ancient groups (see a mathematical definition in ). An attractive feature of AGs is that they can summarize important features of population history without requiring specification of all parameters such as population sizes, split times, mixture times, and distinguishing between sudden splits or drawn-out separations. All these parameters describe important features of demographic history and are fit by many methods for fitting demographic models. However, the fact that it is possible to factor this difficult problem by first inferring important aspects of the topology (AGs fitted to allele frequency correlation statistics), and then fitting additional demographic parameters to data such as site frequency spectra, simplifies demographic inference. AGs thus serve both as conceptual frameworks that allow us to think about the relationships of populations deep in time, and as mathematical models we can fit to genetic data.

2AGs are fitted to f-statistics. For convenience, below we use a concise definition of f-statistics by: ‘The most general definition is that of the f4-statistic f4(A, B; C, D), which measures the average correlation in allele frequency differences between (1) populations A and B and (2) populations C and D that is, (pA – pB)∗ (pC – pD), for allele frequencies p, typically averaged over many biallelic single-nucleotide polymorphisms. This f4-statistic is the same as the D-statistic up to a normalization factor.’ The other f-statistics (f2 and f3) can be defined as special cases of f4-statistics: f2(A, B) = f4(A, B; A, B) and f3(A; B, C) = f4(A, B; A, C). f4-Statistics can be written as linear combinations of f3- or f2-statistics, and f3-statistics can be written as linear combinations of f4- and f2-statistics. f2-, f3-, and f4-statistics have straightforward interpretations in terms of drift edges along the tree, see Figure 2 in and Appendix 1—figure 1b. A challenge for fitting AG models is that they are often not uniquely constrained by the data, with many providing equally good fits to the f2-, f3-, and f4-statistics used to constrain them within the limits of statistical resolution. Previously published methods for finding fitting AGs (mainly qpGraph, and TreeMix,; ) were not well equipped to handle the large range of equally well-fitting models for three reasons: (1) They did not reliably provide information on whether there is a uniquely fitting parsimonious model or alternatively whether there are many models that fit equally well to the limits of statistical resolution, (2) they did not provide formal goodness-of-fit tests, and related to this, (3) they did not provide tests for whether the difference between the fits of any two models is statistically significant. As a consequence, and as we demonstrate in what follows, many published AG models have been interpreted as providing more confidence than is merited about the extent to which genetic data allows us to disentangle ancestral relationships.

3To appreciate these problems, we first need to consider two main approaches that were utilized to study demographic history with AGs.

4The first approach is to identify AGs automatically, either without human intervention or with guidance. It is possible in theory to exhaustively test all possible graphs for a given set of populations and pre-specified number of admixture events, as implemented, for example, in the admixturegraph R package. An exhaustive approach can provide a complete view of the range of models that are consistent with the data for a specified level of parsimony (total number of admixture events allowed in the graph), which is not biased by the algorithm used to explore the space of possible AGs. However, this approach is limited to small graphs (typically up to six groups, two admixture events) due to the rapid increase in the number of possible AGs as the number of populations and admixture events grows. As we show in our discussion of case studies, the simple models explored with an exhaustive approach can lead to misleading conclusions about population history because not including additional populations can blind users to additional mixture events that occurred (and whose existence is revealed by examining data from additional populations). Furthermore, models with additional admixture events that are qualitatively different to the best-fitting parsimonious graph and that capture the true history, will sometimes be completely missed when constraining the number of gene flows. Alternatively, the programs TreeMix, MixMapper, miqoGraph, and AdmixtureBayes ( preprint) all address the problem of how to rapidly explore the vast space of AGs relating a set of populations by applying algorithmic ideas or heuristics; all of these methods speed up model search by orders of magnitude.

5The second approach to fitting AGs is to manually build them up by grafting additional populations onto simpler smaller graphs that fit the data. This approach involves stepwise addition of populations in an order that is chosen based on the best judgment of the user, and for each newly added population involves adding admixture events or tweaks in the graph until a fit is obtained; the user then moves on to adding the next population (see;;;;;;;;;;;;;;;;;;;; for examples). The program qpGraph in the ADMIXTOOLS package has been the most common computational method used for testing fits of individual AGs. Most AGs in the literature have been constructed manually in this way, often acknowledging the existence of alternative models by presenting plausible models side-by-side, and this approach has been the basis for many claims about population history. A strength of this approach is that it takes advantage of human judgment and outside knowledge about what graphs best fit the history of the human or animal populations being analyzed. This external information is powerful as it can incorporate non-genetic evidence such as geographic plausibility and temporal ordering of populations or linguistic similarity, or other genetic data such as estimates of population split times, or shared Y chromosomes, or rejection of proposed scenarios based on joint analysis of much larger numbers of populations than can reasonably be analyzed within a single AG. Thus, while manual approaches explore many orders of magnitude fewer topologies than automatic approaches often do, they still may provide inferences about population history that are more useful than those provided by automatic approaches. These methods’ strength is also their weakness: by relying on intuition, following a manual approach has the potential to validate the biases users have as to what types of histories are most plausible (these may be the only types of histories that will be carefully explored). This can blind users to surprises: to profoundly different topologies that may correspond more closely to the true history, and we discuss examples of this in the Results section.

6In this study, we introduce a new method, findGraphs, that belongs to the first class of algorithms (those for automated AG topology inference). Algorithmic innovations and speedups in findGraphs enable us to explore a much larger proportion of plausible AG space than many other methods reported to date. The findGraphs method combines the advantages of automated and manual topology exploration by allowing users to encode various sources of information as constraints on the space of AGs, which is then explored automatically. However, the main innovations in findGraphs are not computational, but instead conceptual. Instead of finding one or a few AGs fitting the data well, we use findGraphs for exploring AG spaces and assessing if any reliable information on population history can be extracted from a given AG space (defined by a population set and parsimony constraints) in the first place.

7Regardless of the approach used to search through the space of possibly fitting AGs, a challenge in the effort to find a uniquely well-fitting AG (or group of topologically similar AGs) is that it has been difficult to quantify the absolute goodness of fit of a model to date. We have not been entirely successful with this and are not aware of other work that has been successful. It is also difficult to assess the relative fits of multiple models, especially if they differ in complexity. Performance gains relative to the original implementation of qpGraph allow us to address this problem by obtaining bootstrap confidence intervals and p-values for estimated parameters of single models, as well as for the difference in fit quality of two models (see Appendix 1, Sections 1.B.3 and 2.E). In combination with the approach to automating the search of well-fitting AGs, this leads to a situation where we are able to find and test a large number of models, many of which fit equally well despite often having very different topological features. Published approaches to comparing the fits of AG models based on Akaike information criterion (AIC) or Bayesian information criterion (BIC), see; have the problem that it is often not clear what the effective number of degrees of freedom is in the two models being compared since in the case of AGs it depends not only on the number of graph edges, but also on graph topology.

8The methods for automated graph topology inference and model comparison relying on bootstrap resampling are implemented in ADMIXTOOLS 2, a comprehensive platform for learning about population history from f-statistics. It is built to provide a stand-alone workspace for research in this area and is implemented as an R package. For all computations, ADMIXTOOLS 2 exhibits large speedups relative to previously published platforms for f-statistic analysis (e.g., popstats and ADMIXTOOLS version 6.0 which we call ‘Classic ADMIXTOOLS’ in what follows to distinguish it from updated ADMIXTOOLS version 7.0.2 which implements some of the speedup ideas also implemented in ADMIXTOOLS 2). This is achieved by deploying a series of algorithmic improvements, most notably storage of precomputed f-statistics in random access memory, which avoids having to rely on reading in extremely large genotype matrices to perform most computations. In addition to the new algorithmic ideas allowing efficient searching through the space of AGs and comparing the fits of two AGs, ADMIXTOOLS 2 also provides a solution to the question of which parameters of an AG are identifiable in the limit of infinite data. Methodological details are presented in Appendix 1, and below we focus on documenting problems of AG inference on simulated data and revisiting AGs from the literature to understand the extent to which methodological challenges with AG fitting biased previous studies.

9First, we explored the performance of the findGraphs method for automated topology inference on simulated AGs of random topology, focusing on the following questions: (1) among findGraphs results, how common are AGs fitting nominally or significantly better than the true one but different topologically; (2) what is the degree of topological diversity among these models fitting the data better than the true one? For this purpose, we simulated AGs of four complexity classes using msprime v.1.1.1: eight or nine non-outgroup populations, and four or five admixture events. Only simulations where pairwise FST for groups were in the range characteristic for anatomically modern and archaic humans were selected for further analysis, resulting in 20 random topologies per complexity class, each including a distant outgroup that facilitates automated exploration of the topology space.

10We ran findGraphs on each simulated dataset starting from random graphs and pre-specifying the true number of admixture events (n), or n − 1, or n + 1 events. For each of these graph complexity levels, we performed 100 independent findGraphs runs and recorded 5 AGs from each run having the best log-likelihood (LL) scores. Topologically redundant AGs were discarded, and for the remaining AGs we calculated worst f4-statistic residuals (WR) and tested if the newly found models fit significantly better than the true model, using the bootstrap model comparison method developed in this study (see Appendix 1, Sections 1.B.3 and 2.E). In Figure 1a–c, we show the following statistics for each simulated AG, summarized across simulated complexity classes and parsimony levels allowed at the stage of topology exploration: fraction of topologies found with findGraphs that fit better than the true AG (according to LL score), or that fit significantly better than the true AG, or those with plausible absolute fits (WR < 3 SE). It is clear that for the great majority of simulated datasets, even a shallow exploration of the topology space with findGraph (100 independent runs) uncovers AGs that fit nominally better than the true topology (Figure 1a and d) and are topologically diverse (see Figure 1e for examples). When allowing for n admixture events, at least one AG fitting significantly better than the true one was found for 60% of simulated datasets. When n + 1 admixture events were allowed, this grew to 100% (all 80 datasets). It should be noted that some admixture events are indistinguishable with f-statistics; for instance, successive gene flows between two lineages, with no other edges branching off between the gene flows. If such gene flows were included in the random topologies we simulated, AGs with n events were overly complex for representing the true history. Thus, if we are dealing with random histories, choosing an optimal complexity class for topology search is not straightforward.

11These results on simulated data raise concerns about the extent to which fitting AG topologies provide reliable information about population history. Even for histories including eight or nine groups, an outgroup, and four or five pulse-like admixture events, perfect diploid data, and groups as differentiated as Neanderthals and anatomically modern humans—a complexity class that is simpler than many models fitted to real genetic data in published papers—models fitting the data as well as or better than the true one are common, and their topological diversity is in most cases so high that it precludes consensus inference of topology by analysis of multiple topologies. As we demonstrate on another set of simulated data in Appendix 1 (Section 1.B.1), the probability of finding a ‘wrong’ model that fits better than the true one grows with increasing graph complexity, and that effect is reproduced with both findGraphs and TreeMix (Appendix 1—figure 2). We expect this problem to be even more acute when researchers are dealing with realistic complex histories. However, geneticists often rely on external constraints on AG topologies (such as temporal plausibility of a topology, results from qpAdm modelling, f4-statistics, PCA, ADMIXTURE, geographical, archaeological, and linguistic considerations) that were not used for filtering the results of our topology searches on simulated data. Thus, it is possible in principle that published AG models are more robust than our results on simulated data, and we explore this issue in depth in the next section and in Appendix 2.

12We studied AGs from eight publications with the goal of comparing published models to models identified by our algorithm for automatically inferring optimal (best-fitting) AGs (Table 1). In all but one study, qpGraph or its automated reimplementation (admixturegraph, ) was used for fitting topologies to genetic data, while relied on the automated OrientAGraph method. The main question we were interested in is whether we can find alternative models which (1) fit as well as, or better than the published graph, (2) differ in important ways from the published graph, and (3) cannot immediately be rejected based on other evidence such as temporal plausibility. The studies were selected according to the criterion that an AG model inferred in the study is used as primary evidence for at least one statement about population history in the main text of the study. In other words, the AG method was used in the original studies to support new conclusions about population history, and not simply to show that there is a model that exists that does not contradict results of other genetic analyses, an approach that is a valid use of AGs and has been taken in some studies (e.g.,;; ). There are many published studies that could have been included in our re-evaluation exercise as they meet our key criterion (e.g.,;;;;,;;; ). However, critical re-evaluation of each published graph is an intensive process, and the sample of studies we revisited is diverse enough to identify some general patterns.

13Here we present a high-level summary of these analyses. Discussion of individual case studies follows below, and for details see the exposition in Appendix 2.

14For 19 out of 22 published graphs we examined, we were able to find at least one, but usually many, graphs of the same complexity (number of groups and admixture events), with an LL score that was nominally better than that of the published graph (see results for 11 selected graphs in Table 1 and full results for all 22 graphs in Supplementary file 1). The 22 graphs were drawn from the 8 publications as there were multiple final graphs presented in some of the publications, or we examined selected intermediates in the model construction process, or we introduced an outgroup not used in the original study, or we tested additional graph complexity classes dropping ‘unnecessary’ admixture events.

15These alternative graphs often fit not significantly better than the published one after taking into account variability across single-nucleotide polymorphisms (SNPs) via bootstrapping. In the following cases, at least one model that fits significantly better than the published one according to our bootstrap model comparison method was found: the Bergström et al. and Lazaridis et al. 7-population graphs; the graph with 3 admixture events; the Hajdinjak et al. graphs with or without adding a chimpanzee outgroup; the. intermediate graphs with 7 groups and 4 admixture events and with 10 groups and 8 admixture events; the Wang et al. 12-population graph; and the Sikora et al. graphs for West Eurasians and for East Eurasians with 10 or 6 admixture events (Supplementary file 1). In nearly all cases (except for the Lazaridis et al. six-population graph, Shinde et al. graph with eight populations and three admixture events, and the. graph with four admixture events), we also identified many additional graphs that fit the data not significantly worse than the published ones. In every example, some of these graphs have topologies that are qualitatively different in important ways from those of the published graphs. Features such as which populations are admixed or unadmixed, direction of gene flow, or the order of split events, if not constrained a priori, are generally not the same between alternative fitting models for the same populations. This result agrees with the expectation from our exploration of simulated AGs (Figure 1). While some of these graphs can be rejected since their topologies appear highly unlikely because of non-genetic or unrelated genetic evidence, for all of the publications except one, there are alternative equally-well-or-better-fitting graphs we identified and examined manually that differ in qualitatively important ways with regard to the implications about history, are temporally plausible (for instance, very ancient populations do not receive gene flows from sources closely related to much less ancient groups), and not obviously wrong based on other lines of evidence. These findings and the results on simulated AGs suggest that complex AG models, even with a very good fit to the data, often differ in important ways from true population histories.

16The previous statements are valid if the original parsimony constraints are applied, that is, if the graph complexity (the number of admixture events) is not altered. Below in selected case studies (Shinde et al.,.) we also explore the effect of relaxing the parsimony constraint. Table 1 and Figure 2 summarize these results for one or a few graphs from each publication, while Supplementary file 1 contains the full results for all studied graphs and setups. Table 2 summarizes our assessment of inferences in the original publications that were supported by the published graphs.

17To identify alternative models, we ran many iterations of findGraphs for each set of input populations, constraints, and the number of admixture events being fit to the data, and we selected the best-fitting graph in each iteration, that is, a graph with the lowest LL score. Each iteration was initiated from a random graph. The algorithm is non-deterministic so that in each iteration it takes a different trajectory through graph space, possibly terminating in a different final best graph. The number of admixture events in the initial random graphs and in the output graphs was always kept equal to that of the published graph. For each example, we counted how many distinct topologies were found with significantly or non-significantly better or worse LL scores than that of the published graph (Table 1, Supplementary file 1). To obtain a formally correct comparison of model fit, the published graph and each alternative model were fitted to resampled replicates of the dataset and the resulting LL score distributions were compared (see Appendix 1, Sections 1.B.3 and 2.E). As shown in Figure 2, for four of the eight publications we re-analyzed, the LL score of the published graph run on the full data is better than almost all the bootstrap replicates on the same data (it falls below the fifth percentile), which is a sign of overfitting, and underscores the importance of applying bootstrap to assess the robustness of fitted models and conclusions drawn from them.

18The fraction of graphs with scores better than the score of the published graph should not be overinterpreted, as it is influenced by the findGraphs algorithm, which does not guarantee ergodic sampling from the space of well-fitting AGs. In particular, it is possible that despite findGraph’s strategies for efficiently identifying classes of well-fitting AGs (see Appendix 1, Sections 1.B.1 and 2.C), it has a bias toward missing particular classes of graph topologies. However, even one alternative graph which is not significantly worse-fitting than the published graph suggests that we are not able to identify a single best-fitting model. Many of these alternatives, despite providing a good fit to the data, appear unlikely, for example, because they suggest that Paleolithic-era humans are mixed between different lineages closely related to present-day humans. We were mainly interested in alternative models which are also plausible, and so we constrained the space of allowed topologies in findGraphs to those we considered plausible a priori, in cases where this was necessary for reducing the search space size. Constraints were either integrated into the topology search itself, or were applied to outcomes of unconstrained searches, as detailed below.

19Below we summarize our key findings and the methodological implications from our re-analysis of the eight published datasets. For more detailed discussions see Appendix 2.

20The AG for ancient and present-day dogs in Figure 1e of includes an outgroup, six other groups and three admixture events (Figure 3a, Figure 3—source data 1). A best-fitting newly found graph fits the data nominally better than the published one (two-tailed empirical p-value = 0.332), and it bears a closer resemblance to the human population history (Figure 3—source data 2). In this new seven-population model (Figure 3a), both American and Siberian dog lineages represent a mixture between groups related to the Asian and East European dog lineages, and robust genetic results suggest that in the time horizon investigated in the original publication (after ca. 10,900 years ago) nearly all Siberian and all American human populations were admixed between groups most closely related to Europeans and East Asians. According to this model, East Mediterranean dogs are modeled as a mixture of a basal branch (splitting deeper than the divergence of the Asian and European dogs) and West European dogs, again in agreement with current models of genetic history of West Asian human populations who are modeled as a mixture of ‘basal Eurasians’ and West European hunter–gatherers. Although greater congruence with human history increases the plausibility of findGraph’s newly identified model relative to the published model, to make unbiased comparisons between the history of the two species, model selection should be done strictly independently for each species, and so the genetic data alone does not favor one model more than another.

21To explain why the original paper on the population history of dogs missed the model that findGraphs identified, we observe that the AG search was exhaustive under the parsimony constraint (no more than two admixture events for six populations, with the seventh group added at a later stage without an exhaustive topology search), and thus missed the potentially true topology including three admixture events for these six populations. This case study also illustrates that even in a relatively low complexity context (seven groups and three admixture events) applying manual approaches for finding optimal models is risky. When any new group such as an Early Neolithic dog from Germany is added to the model, it may introduce crucial new constraints into the system, and re-exploring the whole graph space in an automated way is necessary to avoid missing the true model. In contrast, mapping a newly added group on a simple skeleton graph (even when that skeleton is a uniquely best-fitting model like in ) may yield a topology that is at odds with the true history.

22The graph in Figure 3 (in ) suggested that present-day Europeans are derived from at least three populations that are very much differentiated genetically: West European hunter–gatherers (WHG), early European farmers (EEF), and Siberian hunter–gatherers from the same lineage as that of the Mal’ta boy who lived about 24,000 years ago (MA1). For seven-population graphs with four admixture events, we found 40 out of 306 distinct graphs with a score better than that of the published graph (10 of those graphs are shown in Figure 3—source data 3). The best-fitting newly found model and two other models fit the data significantly better than the published model (Supplementary file 1), but their topology is qualitatively very similar to that of the published graph (Figure 3—source data 3). In the best-fitting newly found model, French and Karitiana share some drift to the exclusion of MA1, while in the published model the source of MA1-related ancestry in French is closer to MA1 than to Karitiana.

23It is important to point out that not all of the 40 alternative graphs that fit nominally or significantly better than the published one are consistent with the conclusion that modern European populations are admixed between three different ancestral populations (Figure 3—source data 3). According to the fifth alternative graph in Figure 3—source data 3 that fits nominally better than the published model (p-value = 0.464), the present-day European population was formed by admixture of an MA1-related lineage and a European Neolithic-related lineage, with no WHG contribution. Of course, other lines of evidence make it clear that LBK Stuttgart is a mixture of Anatolian farmer- and WHG Loschbour-related ancestry (e.g.,; ), thus providing external information in favor of the model. The use of such ancillary information in concert with graph exploration is important in order to obtain more confident inferences about population history taking advantage of AGs. We also note that a large group of newly found models (247 graphs) fits not significantly worse than the published one (Supplementary file 1), and those are topologically diverse. Thus, strictly speaking, the AG method on the given dataset cannot be used to prove that the published model is the only one fitting the data.

24The skeleton AG in the original study was constructed manually, and subsequently all possible branching orders (105) within the five-population Iranian farmer/herder-related clade were tested. The published model (Figure 3 in that study) included nine groups and three admixture events, but one group (Belt Cave Mesolithic) had a very high missing data rate. Following the approach of the published paper, we repeated findGraphs analysis both with and without the Belt Cave individual. Thus, we initially explored the following topology classes: 9 groups with 3 admixture events on ca. 19,000 polymorphic sites and 8 groups with 3 admixture events on ca. 470,000 sites (Figure 3—source data 4, Supplementary file 1). The finding that the predominant ancestry component of the Indus Periphery group was the most basal branch in the Iranian farmer clade was a prominent claim of the original study. This finding if correct is important, since it implies that the Iranian-related ancestry in the Indus Valley Civilization genetic grouping (which is the same group as Indus Periphery or IP) split from the Iranian-related ancestry in the first Iranian plateau farmers before the date of the Hajji Firuz farmers, who at ~8000 years ago are among the earliest people living on the Iranian plateau known to have grown West Asian crops. The ancient DNA record combined with radiocarbon dating evidence suggests that beginning around the time of the Hajji Firuz farmers, both West Asian domesticated plants such as wheat and barley, and Anatolian farmer-related admixture, began spreading eastward across the Iranian plateau. If the Iranian-related ancestry in IP was spread eastward into the Indus Valley across the Iranian plateau as part of the same agriculturally associated expansion—perhaps brought by people speaking Indo-European languages as well as introducing West Asian crops—then we would expect to see at least some of the Iranian-related ancestry in IP being a clade with that in Hajji Firuz relative to Ganj Dareh.

25The following groups were admixed by default in the graph models compared in the original study: Hajji Firuz Neolithic (labeled ‘Chalcolithic’ in that study but the dates are Neolithic) and Tepe Hissar Chalcolithic were considered as mixtures of an Anatolian farmer-related lineage and an Iranian farmer-related lineage. Indus Periphery was assumed to be a mixture of an Andamanese-related lineage representing ancient South Indians (ASI) and an Iranian farmer-related lineage. The original study differed from ours since these constraints were introduced manually, but we wanted our topology search to be automatic and to explore a wider range of parameter space.

26To provide power to detect negative f3-statistics useful for constraining the model search, we introduced several modifications to the original group composition (described in Appendix 2) and a new algorithm that makes it possible to compute negative f3-statistics on pseudo-haploid data, but at a cost of removing sites with only one chromosome genotyped in any non-singleton population (see Appendix 1, Section 1.A). We repeated topology search with this set of f-statistics providing additional constraints, performing 4,000 runs of the findGraphs algorithm. The Mota ancient African individual was set as an outgroup and three admixture events were allowed in the eight-population graph. Among 4,000 resulting graphs (one from each findGraphs run), 144 were distinct topologically, and the published model was recovered in 13 runs of 4,000 (Supplementary file 1). Only four distinct topologies fitting nominally better than the published one were found, and those had LL scores almost identical to that of the published AG. These four alternative models (Figure 3—source data 5b) shared all topologically important features of the published model (Figure 3—source data 5). Five other topologies differed in important ways from the published one and emerged as fitting the data worse, but not significantly worse, than the published one (Figure 3—source data 5c).

27But in fact, the AG analysis reported above may not be an adequate exploration of the problem. Although absolute fits of the best models found are good (WR = 2.5 SE), the parsimony constraint allowing only three admixture events precluded correct modeling of basal Eurasian ancestry shared by all West Asian groups or of the Indus Periphery group itself, for which a more complex 3-component admixture model was proposed elsewhere. Concerned that this oversimplification could be causing our search to miss important classes of models, we explored qpAdm models for the Indus Periphery group, following the protocol by (see the dataset composition in Supplementary file 4). Our qpAdm results (Supplementary files 5 and 6, Appendix 2) show that the parsimony assumption that was made when constructing the AG analysis in is contradicted by f-statistic evidence since the simplest fitting qpAdm model for the IP group includes four ancestry sources, not two (Indus Periphery = Ganj Dareh Neolithic + Onge (ASI) + WSHG + Anatolia Neolithic), and indeed Narasimhan et al. themselves showed this when they presented a qpAdm model that was more complex (Ganj Dareh Neolithic + Onge (ASI) + WSHG) than the one used for constraining the AG model comparison (Ganj Dareh Neolithic + Onge (ASI)).

28To explore how the parsimony constraint influences results, we allowed four admixture events in the eight-population graph (Supplementary file 1). Among 4,000 resulting graphs (one from each findGraphs run), 443 were distinct topologically, and 270 had WRs between 2 and 3 SE, that is, fitted the data well. In Figure 3—source data 6b, we show four graphs with four admixture events that model the Indus Periphery group as a mixture of three or four sources, with a significant fraction of its ancestry derived from the Hajji Firuz Neolithic or Tepe Hissar Chalcolithic lineages including both Iranian and Anatolian ancestries. The fits of these models are just slightly different (e.g., LL = 11.7 vs 9.3, both WRs = 2.4 SE) from that of the best-fitting model (Figure 3—source data 6a), and are similar to that of the simpler published graph. Besides these four illustrative graphs, dozens of topologies with very different models for the Indus Periphery group fit the data approximately equally well, suggesting that there is no useful signal in this type of AG analysis when the parsimony constraint is relaxed (this finding is similar to that in our re-analysis of the dog AG in, where relaxation of the parsimony constraint identified equally well-fitting AGs that were very different with regard to their inferences about population history). These results show that at least with regard to the AG analysis, a key historical conclusion of the study (that the predominant genetic component in the Indus Periphery lineage diverged from the Iranian clade prior to the date of the Ganj Dareh Neolithic group at ca. 10 kya and thus prior to the arrival of West Asian crops and Anatolian genetics in Iran) depends on the parsimony assumption, but the preference for three admixture events instead of four is hard to justify based on archaeological or other arguments.

29Why did the AG analysis find support for the IP Iranian-related lineage being the first to split, while our findGraphs analysis did not? study sought to carry out a systematic exploration of the AG space in the same spirit as findGraphs—one of only a few papers in the literature where there has been an attempt to do so—and thus this qualitative difference in findings is notable. We hypothesize that the inconsistency reflects the fact that the deeply diverging WSHG-related ancestry present in the IP genetic grouping at a level of ca. 10% was not taken into account explicitly neither in the AG analysis nor in the admixture-corrected f4-symmetry tests also reported in. The difference in qualitative conclusions may also reflect the fact that the Shinde et al. study was distinguishing between fitting models relying on an LL difference threshold of 4 units (based on the AIC). As discussed in Appendix 1, AIC is not applicable to AGs where the number of independent model parameters is topology-dependent even if the numbers of groups and admixture events are fixed, and models compared with AIC should have the same number of parameters. Thus, we believe that the analysis by Shinde et al. was over-optimistic (as compared to the bootstrap model comparison method we use) about being able to reject models that were in fact plausible using its AG fitting setup.

30In contrast to the other studies revisited in our work, the AG published by was inferred automatically using OrientAGraph. Models with three (Figure 3b in that study) and zero to five (Ext. Data Figure 5a–d in that study) admixture events were shown. The dataset included 10 populations (9 horse populations and donkey as an outgroup) and was based on 7.4 million polymorphic transversion sites with no missing data at the group level. Unlike all the other AGs we re-evaluate in this study whose fits to the data were evaluated in the published studies using qpGraph, the topologies published in were not evaluated for statistical goodness-of-fit, and in fact fit the f-statistic data so poorly that even simple statistics show they cannot be correct (Figure 3b, Figure 3—source data 7a, c, e,, Supplementary file 1). In this case, the approach of using findGraphs to identify alternative topologies with the same number of admixture events that fit the data better is meaningless, as both the published models and the alternative models do not have enough degrees of freedom to accommodate the complexity present in the real data (Figure 3—source data 7). In particular, we found that WR of the published model with three admixture events is 23.9 SE (Figure 3—source data 7a).

31For this reason, we moved to topology searches in more complex model spaces incorporating six to nine admixture events. Temporally plausible models with even a modest fit (WR between 3 and 4 SE) were encountered only among models with eight and nine admixture events (Figure 3—source data 7j–r). discussed five inferences relying fully or partially on their published AGs reported in that study (Table 2). The simplest temporally plausible and best-fitting (WR = 3.4 SE) model we found (eight admixture events, see Figure 3b and the first model in Figure 3—source data 7j) supports inferences 2 and 4, and is incompatible with inferences 1, 3, and 5 (Table 2). We consider this model to be plausible also from the geographical perspective (see Appendix 2 for an interpretation of this topology). We are not arguing here that this AG represents the true history; in fact, it is highly unlikely to be the truth, given how large the space of all possible admixture events is and how much admixture evidently occurred relating all these groups (which makes finding the true model extremely unlikely, see the results on simulated data in Figure 1 and Appendix 1—figure 2). However, our set of 16 temporally plausible and fitting (WR < 4 SE) models with eight or nine admixture events (Figure 3—source data 7j–r) is consistent with some features of the published graph being stable: the features (2) that DOM2 and C-PONT are sister groups, and (4) that there was a gene flow from a deep-branching ghost group to NEO-ANA (Table 2).

32Equally important is our finding that there are plausible models that are inconsistent with other inferences in. (Table 2). For example, 13 of these 16 models are inconsistent with the suggestion that there was no gene flow connecting the CWC group and the cluster maximized in the Western steppe (DOM2, C-PONT, and TURG) (Figure 3—source data 7j–r). In the eight-admixture-event best-fitting AG (Figure 3b, the first model in Figure 3—source data 7j), CWC actually derives appreciable ancestry from the early domestic horse lineage (DOM2) associated with the Sintashta culture to the exclusion of the more distant Yamnaya-associated TURG and C_PONT horses. This scenario presents a parallel to the one observed in humans, with individuals associated with the CWC receiving admixture from Steppe pastoralists albeit in different proportions: ~75% for humans, versus ~20% in horses. These models specifying a substantial Steppe horse contribution to CWC horses are inconsistent with the inference in. that ‘Our results reject the commonly held association between horseback riding and the massive expansion of Yamnaya steppe pastoralists into Europe around 3000 BC’. We are not aware of other lines of evidence in the paper (apart from the fitted AG) that support the claim of no Yamnaya horse impact on CWC horses.

33Another example of a feature of the published graph that turned out to be unstable is the model for the Tarpan horse. Only 8 of 16 temporally plausible and fitting models (Figure 3—source data 7j–r) support the conclusion by Librado et al. that the Tarpan is a mixture of a DOM2-related and a CWC-related lineage. The other eight models suggest that Tarpan is a mixture of a deep lineage and a DOM2-related lineage (Figure 3b, the first model in Figure 3—source data 7j), echoing a hypothesis that Tarpan may be a hybrid with Przewalski-related horses not represented in the AG. Again, we are not arguing here that our alternative model is right—indeed we are nearly certain it is wrong in important aspects—but we are merely pointing out that the complexity of the AG space means that qualitatively quite different conclusions are compatible with the statistics fitted in the published paper.

34The AG inferred by Hajdinjak et al. was constructed manually and incorporated 11 groups and 8 admixture events (Figure 2d in the original study). Most (71.4%) models found with findGraphs fit nominally better, and 15.7% fit significantly better than the published model (Table 1, Supplementary file 1, Figure 2), which has a poor absolute fit on this set of sites and groups (WR = 4.8 SE, Figure 3c, Figure 3—source data 8). The statistics described above and the fact that LL scores on all sites lie outside of the LL distribution on resampled datasets (Figure 2) suggest that models in this complexity class are overfitted, but the published topology emerged as fitting relatively poorly. Overfitting arises naturally during manual graph construction as performed in many studies (not only in, but also in;;;;;;;;;;;;;; ). The graph grew one group at a time, and each newly added group was mapped on to the pre-existing skeleton graph as unadmixed or as a two-way mixture. Another requirement was that all intermediate graphs have good absolute fits (WR below 3 or 4 SE). When the model-building process is constrained in a particular path and fits of all intermediates are required to be good, unnecessary admixture events are often added along the way, and the resulting graph belongs to a complexity class in which models are overfitted.

35’s published graph had three notable features that were interpreted by the authors and used to support some conclusions of the study (Table 2), with the following feature considered the most important: there are gene flows from the lineage found in the ~45,000- to 43,000-year-old Bacho Kiro Initial Upper Paleolithic (IUP) individuals to the Ust’-Ishim, Tianyuan, and GoyetQ116-1 lineages. We identified 1,421 topologies fitting nominally or significantly better than the published model and moved on to inspect 50 best-fitting topologies for temporal plausibility (all of them fitting significantly better than the published model). All non-African individuals included in the model are Upper Paleolithic and their dates are not drastically different in relative terms, from ca. 45 to 30 kya (1,000 years before present). Nevertheless, we considered most gene flows from later- to earlier-attested lineages as temporally implausible (for instance, GoyetQ116-1 (~35 kya) → Ust’-Ishim (~44 kya), Sunghir III (34.5 kya) → Tianyuan (40 kya), etc.) since they imply great antiquity of the later-attested lineages and of all lineages derived from them at least partially.

36Of the 50 topologies inspected, 32 were considered temporally plausible. Of those topologies, none supported the finding of gene flows from the Bacho Kiro IUP lineage specifically into all three of the Ust’-Ishim, Tianyuan, and GoyetQ116-1 lineages. A total of 17 topologies supported features 2 and 3 but were inconsistent with feature 1; and 14 topologies supported feature 3 only (Table 2). Best-fitting representatives of each of these topology classes are shown along with the published model in Figure 3—source data 8. Considering topological diversity among models that are temporally plausible, conform to current knowledge about relationships between modern and archaic humans, and fit significantly better than the published model, we conclude that feature 3 is probably robust but other details of the fitted AG in Hajdinjak et al.—for example, gene flows to the Ust’-Ishim, Tianyuan, and Goyet Q116-1 lineages from sources sharing drift exclusively with the Upper Paleolithic Bacho Kiro lineage—should not be interpreted as providing meaningful inferences about population history of Upper Paleolithic modern humans.

37A central finding of Hajdinjak et al. is that the Bacho Kiro IUP group shares more alleles with present-day East Asians than with Upper Paleolithic Holocene Europeans despite coming from Europe. Specifically, the study documents significantly positive statistics of the form D(an Asian group, Kostenki14; Bacho Kiro IUP, Mbuti). Hajdinjak et al.’s interpretation of this observation is that ‘there was at least some continuity between the earliest modern humans in Europe [Bacho Kiro IUP] and later people in Eurasia [East Asians]’. However, a significant D- or f4-statistic can have multiple explanations. The statistic f4 (Tianyuan, Kostenki14; Bacho Kiro IUP, Mbuti) is fitted equally well by the published 12-population AG (Z-score for the difference between the observed and fitted statistics = 0.64) and by, for example, the AG in Figure 3c (Z-score = 0.94). Under the latter model that fits the data significantly better than the published model (p-value = 0.02), the Bacho Kiro IUP and Tianyuan branches are not connected by a gene flow and do not receive gene flows from a third common source, but the common ancestor of Ust’-Ishim and all European Paleolithic lineages receives an 8% gene flow from a divergent modern human lineage splitting deeper than Bacho Kiro IUP and Tianyuan (Figure 3c, Figure 3—source data 8c). This scenario or some version of it seems archaeologically and geographically plausible and is not disproven by any other line of genetic or non-genetic evidence of which we are aware. It could correspond to a scenario where a primary modern human expansion out of West Asia contributed serially to the major lineages leading to Bacho Kiro, then later East Asians, then Ust’-Ishim, and finally the primary ancestry in later European hunter–gatherers. This has a very different interpretation from the scenario of distinctive shared ancestry between the earliest modern humans in Europe such as Bacho Kiro IUP and later people in East Asia—to the exclusion of later European hunter–gatherers—that is suggested by the Hajdinjak et al. published graph.

38We are not claiming that this specific alternative model is correct—indeed, it is almost certainly not the correct one given the topological complexity of the set of all AGs consistent with the data—but the existence of it and many other models that fit the data makes it clear that we do not yet have a unique historical explanation for the excess sharing of alleles that has been documented between some Upper Paleolithic European groups (Bacho Kiro IUP,, and GoyetQ116-1, and ) and all East Asians.

39The AG in the original study was constructed manually (Extended Data Figure 4 in that study) and is very complex (12 groups and 12 admixture events): it exists in a space of ~1044 topologies of this complexity. We note that one admixture event was added by to account for potential modern DNA contamination in ancient Shum Laka individuals, and removing it caused a negligible difference in the fit of the published model (Supplementary file 1). Thus, to decrease the complexity of the graph search space, we considered graphs with 12 groups and 11 admixture events. Among 2,000 newly found models, 11.9% fit nominally (but not significantly) better than the published model (Table 1, Supplementary file 1, Figure 2), and absolute fits of 36.7% of novel models are good (WR <3 SE). These metrics, along with the fact that LL scores on all sites lie outside of the LL distribution on resampled datasets (Figure 2), suggest that models in this complexity class, including the published model, are overfitted. Of the AGs we re-evaluate in this study, the graph from shares with the graphs from;;, evidence of being overfitted (Figure 2).

40Below we discuss four prominent features of the AG published in the original study (that were interpreted by the authors and used to support some conclusions of the study) and the extent to which these features consistently replicate across the large number of fitting 12-population graphs with 11 admixture events (Table 2). High topological diversity is observed among temporally plausible newly found AGs (see an example in Figure 4a and further topologies in Figure 4—source data 1). Considering extreme cases, two AGs completely lacked support for three features of the published graph (Figure 4a, Figure 4—source data 1c), and one graph supported all four features of the published graph fully (Figure 4—source data 1q, the second model). There are some graphs where defining two distinct ancestral lineages maximized in West Africans and in Mbuti and Biaka (features 1 and 2, Table 2) is essentially impossible since all or nearly all Africans are modeled as a mixture of at least two deep lineages (see alternative graph no. 4 shown in Figure 4—source data 1d, the second model). In some graphs there is no single lineage specific to rainforest hunter–gatherers (Biaka, Mbuti, and Shum Laka) since the primary ancestries in these groups form independent deep branches in the African graph (see Figure 4a and graph no. 16 shown in Figure 4—source data 1j, the second model). The ghost modern and super-archaic gene flows to Africans also had no universal support in the set of alternative graphs we examined (see, for example, Figure 4a and Figure 4—source data 1c).

41Considering the high degree of topological diversity among models that are temporally plausible, conform to known findings about relationships between modern and archaic humans, and fit nominally better than the published model, we conclude that none of the four AG features from the original study are consistently supported by our re-analysis (Table 2). This situation may be attributed to (1) overfitting and/or to (2) the lack of information in the dataset (in the combination of groups and SNP sites) and/or to (3) inherent limitations of f-statistics, when distinct topologies predict identical f-statistics. Our results highlight the mystery around the highly distinctive genetic ancestry of the Shum Laka individuals themselves, who represent the newly reported data in the study, and represent a highly important set of genetic datapoints that was not available prior to the study. The ancestral relationships of these four individuals to rainforest hunter–gatherers, and to the primary lineage in present-day West Africans, remains an open question, one whose resolution promises meaningful new insights into modern human population history.

42The AG inferred by was constructed manually, and the final graph (Extended Data Figure 6 in the original study) included 12 groups and 8 admixture events. We applied several constraints on the graph space exploration process all of which were shared with the Wang et al. graphs (Supplementary file 1). An important feature of the published graphs in that was remarked upon in the study is admixture from a source related to Andamanese hunter–gatherers that is almost universal in East Asians (Table 2). For example, the abstract states ‘Hunter-gatherers from Japan, the Amur River Basin, and people of Neolithic and Iron Age Taiwan and the Tibetan Plateau are linked by a deeply splitting lineage that probably reflects a coastal migration during the Late Pleistocene epoch.’ We performed 2,000 findGraphs iterations and obtained 1,778 distinct topologies satisfying all the constraints, nearly all of them (1,724) fitting nominally better than the published model, and 12.6% fitting significantly better (Table 1, Supplementary file 1). The models were ranked by LL scores, and 56 highest-ranking topologies, all of them fitting significantly better than the published one, were assessed for temporal plausibility, and 20 topologies were considered temporally plausible (all of them are shown in Figure 4—source data 2). According to these topologies, 0–2 East Asian groups had a fraction of their ancestry derived from a source specifically related to Onge, and 19 topologies included gene flows from the European (Loschbour)-related branch to all 8 East Asian groups (Figure 4—source data 2). The inferred topological relationships among East Asians are variable in this group of AGs, and we decided to apply further constraints that guided model ranking and elimination by Wang et al., based on considerations from archaeological evidence, Y chromosome haplogroup divergence patterns, and population split time estimation (see Appendix 2 for details). Applying these three additional constraints, we identified two models (among the 56 subjected to manual inspection) that satisfied all of them. The highestranking of those models is shown in Figure 4b and Figure 4—source data 2c (the second model), and it includes a 13% (deeply) European-related gene flow to the common ancestor of all East Asians, and gene flows from the Onge-related branch to just two East Asian groups: Nepal Chokhopani and China WLR LN. This model fits the data significantly better than the published model (p-value = 0.028). We do not claim that this is the correct model (indeed we are almost certain that it is not given the high topological diversity of fitting models), but it is not obviously wrong and differs in qualitatively important ways from the published one.

43The AG provides an illuminating example that helps us to understand the value added by AG construction. The AG construction process in Wang et al. followed a philosophy of not relying entirely on the allele frequency correlation data (not treating the genetic data as independent to explore how much new insight could come from genetic data alone). Instead, the study integrated other lines of genetic evidence as well as linguistic and archaeological insights explicitly into the AG construction process, with the goal of identifying models consistent with multiple lines of evidence. The fact that after this procedure a fitting graph was obtained is not of great interest, as it is essentially always possible to obtain a fit to allele frequency correlation data when enough admixture events are added. The important question is whether any of the emergent features of the graph that were not applied as constraints in the construction process—for example the evidence of ubiquitous Andamanese-related gene flow throughout East Asia suggesting a coastal route expansion that admixed with an interior route expansion proxied by Tianyuan—were stably inferred. Our analysis does not come to this finding consistently among well-fitting and plausible AGs. We conclude that this important feature of the published graph is not supported by f-statistic analysis alone (Table 2), and indeed we are not aware of a single feature of the AG that is stably inferred beyond the constraints applied to build it.

44Two AGs inferred by were constructed manually based on an SNP set derived from whole-genome shotgun data and incorporated 12 or 13 groups and 10 admixture events (Extended Data Figure 3f in the original study). One graph was focused on West Eurasians, and the other one on East Eurasians, and both included a Neanderthal, a Denisovan, and an African group (Dinka). Although the chimpanzee outgroup was not included in the original graphs, we added it as it drastically constrains the topology search space. In contrast to most other published graphs discussed above, gene flows in the graphs inferred by Sikora et al. do not have equal standing: four low-level gene flows (0–1%) connect the Neanderthal lineage to Upper Paleolithic lineages. We repeated each topology search under two alternative settings: either keeping the number of admixture events at 10 to match the published graphs, or at 6 to match simplified versions of the published graphs lacking these low-level Neanderthal gene flows. We performed that modification to simplify the search space and to alleviate the overfitting problem which becomes severe if 10 gene flows across the graph are allowed (Supplementary file 1).

45In the case of the "Western" graphs with 6 admixture events, 1,000 topology search iterations were performed, 894 distinct topologies were found, 4 models fit significantly better, and 151 models fit nominally better than the published one (Table 1, Supplementary file 1). We inspected those 155 topologies and identified 29 topologies (Figure 4—source data 3) that are temporally plausible. Sikora et al. came to the following striking conclusion relying on the "Western" AG (Table 2): the Mal’ta (MA1_ANE) lineage received a gene flow from the Caucasus hunter–gatherer (CHG) lineage. However, in our findGraphs exploration this direction of gene flow (CHG → Mal’ta) was supported by two of the 29 topologies, and the opposite gene flow direction (from the Mal’ta and East European hunter–gatherer lineages to CHG) was supported by the remaining 27 plausible topologies (Figure 4—source data 3). The highest-ranking plausible topology (Figure 4c) has a fit that is not significantly different from that of the simplified published model with six admixture events (p-value = 0.392). We note that the gene flow direction contradicting the graph by Sikora et al. was supported by published qpAdm analyses, and qpAdm is not affected by the same model degeneracy issues that are the focus of this study. Considering the topological diversity among models that are temporally plausible, conform to robust findings about relationships between modern and archaic humans, and fit nominally better than the published model, we conclude that the direction of the Mal’ta-CHG gene flow cannot be resolved by AG analysis (Table 2).

46Some important conclusions based on the "Eastern" graph also do not replicate across all plausible AGs (Table 2). In the case of the "Eastern" graphs with 6 admixture events, 4,446 topology search iterations were performed, and 2,785 distinct topologies were found. Only 3 topologies fit significantly and 13 nominally better than the published one, and 9.8% of topologies fit not significantly worse than the published one (Table 1, Supplementary file 1). Of the AGs belonging to these groups, we inspected 116 best-fitting ones and identified 97 AGs that are temporally plausible. The Sikora et al. "Eastern" AG had three distinctive features that were used to support some conclusions of the study (Table 2). Only feature 2 was universally supported by all the 97 plausible alternative models fitting significantly better, nominally better, or not significantly worse than the simplified published model, while feature 3 was supported by 83 of 97 plausible models, and feature 1 was supported by 28 of 97 plausible models (Table 2). We plotted 14 plausible graphs as examples of topologies supporting all three features, two features, or one feature of the published graph (Figure 4—source data 4). We note that all the "Eastern" graphs discussed here, both the published and alternative ones, have relatively poor absolute fits with WR above 4 or 5 SE. Increasing the number of gene flows to 10 allowed us to reach much better absolute fits (with WR as low as 2.42 SE), but that resulted in high topological diversity (on a par with some other case studies discussed above).

47AGs represent a conceptually powerful framework for thinking about demographic history, but, as we demonstrate in this study (see also Appendix 2), the practice of manually constructing a small number of complex models without exploring AG space in an automated way can lead to overconfidence in the validity of these models. An ideal outcome of an AG model exploration exercise would be the identification of a model or a group of topologically very similar models which fit the data well and significantly better than all alternative models with the same number of admixture events; however, this is almost never achieved for graphs with more than eight populations and three admixture events in our experience (Appendix 2), and even this approach can lead to potentially unstable results as relaxing the assumption of parsimony (that fewer admixture events is more likely) can lead to qualitatively quite different equally well-fitting topologies as in our re-analysis of the Bergström et al. and Shinde et al. datasets. Most of the examples of AGs in eight recently published studies we revisited do not fit this ideal pattern, as we were able to identify many topologically different alternative models that could not easily be rejected based on temporal plausibility or other constraints (Figure 3—source data 4, Figure 3—source data 5, Figure 3—source data 6, Figure 3—source data 7, Figure 3—source data 8, Figure 4—source data 1, Figure 4—source data 2, Figure 4—source data 3, Figure 4—source data 4). In particular, for all studies except (under a strict parsimony assumption however), we identified AGs that were not significantly worse fitting than the published ones, and with topological features that were different in qualitatively important ways. There were also some more encouraging findings of the exercise we performed to re-evaluate published models. For example, at least one of the key inferences about population history relying on AG modeling were stable for all analyzed models for the., Hajdinjak et al., Shinde et al. (under the parsimony assumption), and Sikora et al. (simplified "Eastern" graph) studies (Appendix 2). The existence of some stable features in these graphs helps to point the way toward a protocol that we believe should be applied in all future studies that use AG fitting exercises to support claims about population history.

48We propose the following tentative protocol to identify features of fitting AGs that are stable enough to be used to make inferences about population history.

49For a given combination of populations, carry out an initial scan using findGraphs to identify reasonable parameter values for the number of allowed admixture events (graph complexity class). For example, run findGraphs allowing between zero and eight admixture events (100 algorithm iterations per graph complexity class), saving one or a few best-fitting AGs after each iteration. The smallest number of admixture events that yields models where the (negative) LL score or the worst f-statistic residual is lower than some threshold can then be explored more deeply by running more iterations of findGraphs.

50Run findGraphs on the chosen complexity class, where some of the resulting graphs should be inspected manually to determine whether they could in principle be historically plausible models. Implausible models (e.g., models where a very ancient population appears to be admixed between two modern populations) can be filtered out by imposing topological constraints. If no or only a few graphs remain, findGraphs can be run again under these constraints. This can be repeated until one or more graphs with an acceptable LL score or worst residual has been identified. At this stage, apply the bootstrap method to determine whether the best-fitting graph is significantly better than the next best-fitting graph. If it is not, identify a set of graphs which are not clearly worse than the best-fitting graph by performing bootstrap model comparison for many model pairs.

51Researchers should compare the resulting graphs to each other with the goal of identifying common features. Although ADMIXTOOLS 2 includes automated tools for cataloguing common topological features (Appendix 1, Sections 1.B.5 and 2.G), we found a manual approach to be valuable as the fitted parameters (especially admixture proportions) are as important for this task as graph topology.

52Once a set of fitting graphs and stable topological features shared between them is identified, researchers should carry out a findGraphs exploration of the space of graphs with one additional admixture event. If inferences are stable even when fitting graphs with one more level of complexity than the graphs with the minimal number of admixture events needed to fit the data, this increases confidence in the inferences. Furthermore, addition of a new population may introduce crucial information to an existing set of populations, which can change the space of fitting topologies in a profound way, as in our re-analysis of the data from (Figure 3a, Figure 3—source data 1). Thus, it is advisable that the topology optimization procedure is repeated on several alternative population sets, in addition to considering models that allow an additional admixture event beyond the minimum required for parsimony, to explore if inferences about topology change qualitatively.

53AGs fitted with f-statistics do not distinguish between time and population size as the two factors affecting genetic drift. Moreover, many different complex genetic histories for a set of populations can result in the exact same expected f-statistics. This provides an opportunity to further constrain a model fitting procedure. Methods that take advantage of information from the site frequency spectra (momi2, fastsimcoal,, ) or derived site patterns, a special case of site frequency spectra (Legofit, ), can supply alternative information not captured by f-statistics (further information can come from methods that fit haplotype divergence patterns such as MSMC, and SMC++,, or inferences based on fitted gene trees such as RELATE,, and ARGweaver,; ). These tools are too computationally intensive to explore a large number of models, but the advantages of the different approaches can be combined by first identifying a set of candidate models using findGraphs, and then testing these candidate models with other methods. This approach is also expected to help address overfitting since different data types almost always include different variable site sets.

54We believe that researchers should only begin to make strong claims about population history with AGs once a protocol such as we propose is applied.

55We see the guidelines above as analogous in spirit to the protocols that were introduced in medical genetics at a time of the reproducibility crisis in the field of candidate gene association studies. Many studies looking for risk factors for common, complex diseases resulted in publications with marginally significant p-values without correcting for multiple hypothesis testing that was implicitly performed due to many candidate genes being tested and only those with significant findings being published. Unsurprisingly, most of these claims failed to replicate in follow-up studies in independent sets of samples. The human medical genetic community addressed this challenge by coming together to support a rigorous set of commonly accepted standards for declaring genome-wide statistical significance, such as the requirement that p-values be corrected for the effective number of independent common variants in the genome and requiring correction for the known confounders of population structure and undocumented relatedness among individuals.

56Sampling AG space is a useful method for modeling population histories, but finding robust and accurate models can be challenging. As we demonstrated by revisiting a handful of published AGs and re-analyzing the datasets used to fit them, f-statistics are usually insufficient for identifying uniquely fitting AG models, making it necessary to incorporate other sources of evidence. This provides a challenge to previous approaches for automated model building. We investigated several published AG models and, in nearly all cases, found many alternative models, some of which are historically and geographically plausible but contradict conclusions that were derived from the published models. To conduct these analyses, we developed a method for automated AG topology optimization that can incorporate external sources of information as topological constraints. This method is developed in the ADMIXTOOLS 2 framework, which aside from AG modeling, implements many other methods for population history inference based on f-statistics.

57It is important to recognize that the key concern we have highlighted in this study—the fact that there can often be thousands of different topologies that are equally good fits to the allele frequency correlation patterns relating a set of populations—does not invalidate the use of allele frequency correlation testing in many other contexts in which it has been applied to make inferences about population history. For example, negative f3-statistics (‘admixture’ f3-statistics) continue to provide unambiguous evidence for a history of mixture in tested populations, and f4-and D-symmetry statistics remain powerful ways to evaluate whether a tested pair of populations is consistent with descending from a common ancestral population since separation from the ancestors of two groups used for comparison. The qpWave methodology remains a fully valid generalization of f4-statistics, making it possible to test whether a set of populations is consistent with descending from a specified number of ancestral populations (which separated at earlier times from a comparison set of populations). In addition, and the qpAdm extension of qpWave—which allows for estimating proportions of mixtures for the tested population under the assumption that we have data from the source populations for the mixture—remains a valid approach, unaffected by the concerns identified here. Instead of relying on a specific model of deep population relationships, qpAdm relies on an empirically measured covariance matrix of f4-statistics for the analyzed populations, which is highly constraining with respect to estimation of mixture proportions but can be consistent with a wide range of deep history models. All these methods are implemented in ADMIXTOOLS 2.

58Finally, approaches that use AGs to adjust for the covariance structure relating a set of populations without insisting that the particular AG model that is proposed is true with can be useful, for example for the purpose of analyzing shared genetic drift patterns of a group of populations that derive from similar mixtures. One example was a study that attempted to test for different source populations for Neolithic migrations into the Balkans after controlling for different proportions of hunter–gatherer admixture. Another example was a study that attempted to study shared ancestry between different East African forager populations after controlling for different proportions of deeply divergent source populations. However, with respect to the inferences about deep history produced by AGs themselves, our results highlight the importance of caution in proposing specific models of population history that relate a set of groups.