1 of 7EVOLUTIONARY BIOLOGYNeanderthal-Denisovan ancestors interbred with a distantly related homininAlan R. Rogers*, Nathan S. Harris, Alan A. AchenbachPrevious research has shown that modern Eurasians interbred with their Neanderthal and Denisovan predecessors. We show here that hundreds of thousands of years earlier, the ancestors of Neanderthals and Denisovans interbred with their own Eurasian predecessorsâ€”members of a â€œsuperarchaicâ€ population that separated from other humans about 2 million years ago. The superarchaic population was large, with an effective size between 20 and 50 thousand individuals. We confirm previous findings that (i) Denisovans also interbred with superarchaics, (ii) Neanderthals and Denisovans separated early in the middle Pleistocene, (iii) their ancestors endured a bottleneck of population size, and (iv) the Neanderthal population was large at first but then declined in size. We provide qualified support for the view that (v) Neanderthals interbred with the ancestors of modern humans. INTRODUCTIONDuring the past decade, we have learned about interbreeding among hominin populations after 50 thousand years (ka) ago, when modern humans expanded into Eurasia (1â€“3). Here, we focus farther back in time, on events that occurred more than a half million years ago. In this earlier time period, the ancestors of modern humans separated from those of Neanderthals and Denisovans. Somewhat later, Neanderthals and Denisovans separated from each other. The paleon tology and archaeology of this period record important changes, as large-brained hominins appear in Europe and Asia and Acheulean tools appear in Europe (4 , 5 ). It is not clear, however, how these large-brained hominins relate to other populations of archaic or modern humans (6 â€“ 9 ). We studied this period using genetic data from modern Africans and Europeans and from two archaic populations, Neanderthals and Denisovans. Figure1 illustrates our notation. Uppercase letters refer to populations, and combinations such as XY refer to the population ancestral to X and Y . X represents an African population (the Yorubans), Y is a European population, N is Neanderthals, and D is Denisovans. S is an unsampled â€œsuperarchaicâ€ population that is distantly related to other humans. Lowercase letters at the bottom of Fig.1 label â€œnucleotide site patterns.â€ A nucleotide site exhibits site pattern xyn if random nucleotides from populations X, Y, and N carry the derived allele, but those sampled from other populations are ancestral. Site pattern probabilities can be calculated from models of population history, and their frequencies can be estimated from data. Our Legofit (10) software estimates parameters by fitting models to these relative frequencies. Nucleotide site patterns contain only a portion of the information available in genome sequence data. This portion, however, is of particular relevance to the study of deep population history. Site pattern frequencies are unaffected by recent population history because they ignore the within-population component of variation (10). This reduces the number of parameters we must estimate and allows us to focus on the distant past. The current data include two high-coverage Neanderthal genomes: one from the Altai Mountains of Siberia and the other from Vindija Cave in Croatia (11). Rather than assigning the two Neanderthal fossils to separate populations, our model assumes that they inhabited the same population at different times. This implies that our estimates of Neanderthal population size will refer to the Neanderthal metapopulation rather than to any individual subpopulation. The Altai and Vindija Neanderthals appear in site pattern labels as â€œaâ€ and â€œvâ€. Thus, av is the site pattern in which the derived allele appears only in nucleotides sampled from the two Neanderthal genomes. Figure2 shows the site pattern frequencies studied here. In contrast to our previous analysis (12), the current analysis includes singleton site patterns, x, y, v, a, and d, as advocated by Mafessoni and Prfer (13). A simpler tabulation, which excludes the Vindija genome, is included as fig. S2. Greek letters in Fig.1 label episodes of admixture. We label models by concatenating Greek letters to indicate the episodes of admixture they include. For example, model â€œabâ€ includes only episodes a and b. Our model does not include gene flow from Denisovans into moderns because there is little evidence of such gene flow into Europeans (14, 15). Two years ago, we studied a model that included only one episode of admixture: a , which refers to gene flow from Neanderthals into Europeans (12). The left panel of Fig.3 shows the residuals from this model, using the new data. Several are far from zero, suggesting that something is missing from the model (16). Recent literature suggests some of what might be missing. There is evidence for admixture into Denisovans from a superarchaic popu lation, which was distantly related to other humans (2 , 11, 17â€“ 19), and also for admixture from early moderns into Neanderthals (19). These episodes of admixture appear as b and g in Fig.1. Adding b and/or g to the model improved the fit, yet none of the resulting models were satisfactory. For example, model abg implied (implausibly) that superarchaics separated from other hominins 7 million years (Ma) ago. To understand what might still be missing, consider what we know about the early middle Pleistocene, around 600 ka ago. At this time, large-brained hominins appear in Europe, along with Acheulean stone tools (4, 5). They were probably African immigrants, because similar fossils and tools occur earlier in Africa. According to one hypothesis, these early Europeans were Neanderthal ancestors (6, 7). Somewhat earlierâ€”perhaps 750 ka ago [( 8 ), table S12.2]â€”the â€œneandersovanâ€ ancestors of Neanderthals and Denisovans separated from the lineage leading to modern humans. Neandersovans may have separated from an African population and then expanded into Eurasia. Department of Anthropology, University of Utah, Salt Lake City, UT, USA. *Corresponding author. Email: email@example.comCopyright 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).
2 of 7If so, then they would not have been expanding into an empty continent, for Eurasia had been inhabited since 1.85Ma ago (20). Neandersovan immigrants may have met the indigenous superarchaic population of Eurasia. This suggests a fourth episode of admixture, from superarchaics into neandersovans, which appears as d in Fig.1.RESULTSWe considered eight models, all of which include a, and including all combinations of b, g, and/or d. In choosing among complex models, it is important to avoid overfitting. Conventional methods such as Akaikeâ€™s information criterion (21) are not available because we do not have access to the full likelihood function. Instead, we use the bootstrap estimate of predictive error (bepe) (10, 22, 23). The best model is the one with the lowest value of bepe. When no model is clearly superior, it is better to average across several than to choose just one (24). For this purpose, we used bootstrap model averaging (booma) (10, 24). The booma weight of the ith model is the fraction of datasets (including the real data and 50 bootstrap replicates) in which that model â€œwins,â€ i.e., has the lowest value of bepe. The bepe values and booma weights of all models are in Table. The best model is abgd, which includes all four episodes of admixture. It has smaller residuals (Fig. 3, right), the lowest bepe value, and the largest booma weight. One other model, abd, has a positive booma weight, but all others have zero weight. To understand what this means, recall that bootstrap replicates approximate repeated sampling from the process that generated the data. The models with zero weight lose in all replicates, implying that their disadvantage is large compared with variation in repeated sampling. On this basis, we can reject these models. Neither of the two remaining models can be rejected. These results provide strong support for two episodes of admixture (b and d) and qualified support for a third (g). Not only does this support previously reported episodes of gene flow but it also reveals a much older episode, in which neandersovans interbred with superarchaics. Model-averaged parameter estimates, which use the weights in Table 1, are graphed in Fig. 4 and listed in table S1. Episode d , which proposes gene flow from superarchaics into neandersovans, is a novel hypothesis. Before accepting it, we should ask whether the evidence in its favor could be artifactual, reflecting a bias in site pattern frequencies caused by sequencing error or somatic mutations. Sequencing error adds a positive bias to the frequency of each singleton site pattern proportional to the per-nucleotide error rate in the corresponding population (see the Supplementary Materials). Somatic mutations have a similar effect. These biases might explain evidence for episode d, if it were true that larger values of md (the fraction of superarchaic admixture in neandersovans) imply larger frequencies of singleton site patterns. However, Table2 shows that this is not the case. There is no consistent tendency for singleton frequencies to increase with md. Indeed, three of them decrease. Consequently, evidence that md>0 cannot be the result of a positive bias in the frequencies of singleton site patterns. The evidence for d admixture cannot be an artifact of sequencing error or somatic mutations. The superarchaic separation time, TXYNDS, has a point estimate of 2.3 Ma ago. This estimate may be biased upward because our molecular clock assumes a fairly low mutation rate of 0.38 10 per nucleotide site per year. Other authors prefer slightly higher rates ( 25). Although this rate is apparently insensitive to generation time among the great apes, it is sensitive to the age of male puberty. If the average age of puberty during the past 2 Ma were halfway between those of modern humans and chimpanzees, the yearly mutation rate would be close to 0.45 10 [(26), Fig. 2B], and our estimate of TXYNDS would drop to 1.9 Ma, just at the origin of the genus Homo. Under this clock, the 95% confidence interval is 1.8 to 2.2 Ma. If superarchaics separated from an African population, then this separation must have preceded the arrival of superarchaics in Eurasia. Nonetheless, our 1.8 to 2.2Ma interval includes the 1.85Ma date of yv ad xv ad xy ad xyvd xyv a va d ya d yvd yv a xad xvd xv a xyd xy a xyv ad vd va yd ya yv xd xa xv xy d a v y x 0.00 0.05 0.10 0.15 Site p atter n frequency Fig. 2. Observed site pattern frequencies. Horizontal axis shows the relative frequency of each site pattern in random samples consisting of a single haploid genome from each of X, Y, V, A, and D, representing Africa, Europe, Vindija Neanderthal, Altai Neanderthal, Denisovan, and superarchaic. Horizontal lines (which look like dots) are 95% confidence intervals estimated by a moving blocks bootstrap (35). Data: Simons Genome Diversity Project (SGDP) (14) and Max Planck Institute for Evolutionary Anthropology (11). Fig. 1. A population network including four episodes of gene flow, with an em bedded gene genealogy. Upper case letters (X , Y , N , D , and S ) represent populations (Africa, Europe, Neanderthal, Denisovan, and superarchaic). Greek letters label episodes of admixture. d and xyn illustrate two nucleotide site patterns, in which 0 and 1 represent the ancestral and derived alleles. A mutation on the red branch would generate site pattern d. One on the blue branch would generate xyn. For simplicity, this figure refers to Neanderthals with a single letter. Elsewhere, we use two letters to distinguish between the Altai and Vindija Neanderthals.
3 of 7the earliest Eurasian archaeological remains at Dmanisi (20). Thus, superarchaics may descend from the earliest human dispersal into Eurasia, as represented by the Dmanisi fossils. On the other hand, some authors prefer a higher mutation rate of 0.5 10 per year (2 ). Under this clock, the lower end of our confidence interval would be 1.6Ma ago. Thus, our results are also consistent with the view that superarchaics entered Eurasia after the earliest remains at Dmanisi. Parameter NS is the effective size of the superarchaic population. This parameter can be estimated because there are two sources of superarchaic DNA in our sample (b and d ), and this implies that coalescence time within the superarchaic population affects site pattern frequencies. Although this parameter has a broad confidence interval, even the low end implies a fairly large population of about 20,000. This does not require large numbers of superarchaic humans, because effective size can be inflated by geographic population structure (27). Our large estimate may mean that neandersovans and Denisovans received gene flow from two different superarchaic populations. Parameter TND is the separation time of Neanderthals and Denisovans. Our point estimate, 737 ka ago, is remarkably old. Furthermore, the neandersovan population that preceded this split was remarkably small: NND 500. This supports our previous results, which indicated an early separation of Neanderthals and Denisovans and a bottleneck among their ancestors (12). Because our analysis includes two Neanderthal genomes, we can estimate the effective size of the Neanderthal population in two separate epochs. The early epoch extends from TN0 =455 ka to TND=737 ka, and within this epoch, the effective size was large: NN0 16,000. It was smaller during the later epoch: NN1 3400. These results support previous findings that the Neanderthal population was large at first but then declined in size (2,11).DISCUSSIONThis project began with a puzzle. We had argued in 2017 that Neanderthals and Denisovans separated early, that their neandersovan ancestors endured a bottleneck of population size, and that the postseparation Neanderthal population was large (12). That analysis omitted singleton site patterns. Mafessoni and Prfer (13) pointed out that introducing singletons led to different results. In response, Rogers etal. (16) agreed, but also observed that the withsingleton analysis implied that the Denisovan fossil was only 4000 years oldâ€”a result that is plainly wrong. Furthermore, a residual analysis showed that neither of the models under discussion in 2017 fit the data very well (16). Something was apparently missing from both modelsâ€”but what? The present paper provides an answer to that question. yv ad xv ad xy ad xyvd xyv a va d ya d yvd yv a xad xvd xv a xyd xy a xyv ad vd va yd ya yv xd xa xv xy d a v y x .002 0.0000 .002 Observ ed minus fittedSite p atternModel yv ad xv ad xy ad xyvd xyva va d ya d yvd yv a xad xvd xv a xyd xy a xyv ad vd va yd ya yv xd xa xv xy d a v y x .002 0.000 0.002 Obser v ed minus fitted Site p atter n Model Fig. 3. Residuals from models a and abgd. Key: red asterisks, real data; blue circles, 50 bootstrap replicates. Table 1. bepe values and booma weights. Model bepe Weight a 1:16 100 ad 0:87 100 ag 0:62 100 agd 0:44 100 ab 0:18 100 abg 0:17 100 abd 0:15 100.16 abgd 0:13 100.84
4 of 7Our results shed light on the early portion of the middle Pleistocene, about 600 ka ago, when large-brained hominins appear in the fossil record of Europe along with Acheulean stone tools. There is disagreement about how these early Europeans should be interpreted. Some see them as the common ancestors of modern humans and Neanderthals (28), others as an evolutionary dead end, later replaced by immigrants from Africa (29, 30), and others as early representatives of the Neanderthal lineage (6, 7). Our estimates are most consistent with the last of these views. They imply that by 600 ka ago, Neanderthals were already a distinct lineage, separate not only from the modern lineage but also from Denisovans. These results resolve a discrepancy involving human fossils from Sima de los Huesos (SH). Those fossils had been dated to at least 350 ka ago and perhaps 400 to 500 ka ago (31). Genetic evidence showed that they were from a population ancestral to Neanderthals and therefore more recent than the separation of Neanderthals and Denisovans (9). However, genetic evidence also indicated that this split occurred about 381 ka ago [(2), table S12.2]. This was hard to reconcile with the estimated age of the SH fossils. To make matters worse, improved dating methods later showed that the SH fossils are even older, about 600 ka, and much older than the molecular date of the Neanderthal-Denisovan split (32). Our estimates resolve this conflict because they push the date of the split back well beyond the age of the SH fossils. Our estimate of the Neanderthal-Denisovan separation time con flicts with 381 ka ago estimate discussed above (2 , 13). This discrepancy results, in part, from differing calibrations of the molecular clock. Under our clock, the 381-ka date becomes 502 ka (12), but this is still far from our own 737-ka estimate. The remaining discrepancy may reflect differences in our models of history. Misspecified models often generate biased parameter estimates. Our new results on Neanderthal population size differ from those we published in 2017 (12). At that time, we argued that the Neanderthal population was substantially larger than others had estimated. Our new estimates are more in line with those published by others (2, 11). The difference does not result from our new and more elaborate model because we get similar results from model a , which (as in our 2017 model) allows only one episode of gene flow (table S2). Instead, it was including the Vindija Neanderthal genome that made the difference. Without this genome, we still get a large estimate (NN1 11,000), even using model abgd (table S3). This implies that the Neanderthals who contributed DNA to modern Europeans were more similar to the Vindija Neanderthal than to the Altai Neanderthal, as others have also shown (11). Our results revise the date at which superarchaics separated from other humans. One previous estimate put this date between 0.9 and 1.4Ma [(2), p.47], which implied that superarchaics arrived well after the initial human dispersal into Eurasia around 1.9Ma. This required a complex series of population movements between Africa and Eurasia [(33), pp.66 to 71]. Our new estimates do not refute this reconstruction, but they do allow a simpler one, which involves only three expansions of humans from Africa into Eurasia: an expansion of early Homo at about 1.9Ma ago, an expansion of neandersovans at about 700 ka ago, and an expansion of modern humans at about 50 ka ago. Our results indicate that neandersovans interbred with superarchaics early in the middle Pleistocene, shortly after expanding into Eurasia. This is the earliest known admixture between hominin populations. Furthermore, the two populations involved were more distantly related than any pair of human populations previously known to interbreed. According to our estimates, neandersovans and superarchaics had been separate for about 1.2Ma. Later, when superarchaics exchanged genes with Denisovans, the two populations had been separate even longer. By comparison, the Neanderthals and Denisovans who interbred with modern humans had been separate less than 0.7Ma. Table 2. Effect on singleton site pattern frequencies of gene flow (md) from superarchaics into neandersovans. Column 2 shows expected frequencies of singleton site patterns in a model in which md = 0, and all other parameters are as fitted under model abgd. In column 3, all parameters including md are as fitted under this model. Column 4 is obtained by subtracting column 2 from column 3. Expected site pattern frequencies were estimated using legosim with 107 iterations. Site pattern Frequency md = 0 md = 0.034 Difference x 0.15583 0.15174 .00409 y 0.15176 0.14778 .00398 v 0.04974 0.04942 .00032 a 0.03795 0.03798 0.00003 d 0.16051 0.16444 0.00393 m m m m 0.00 0.02 0.04Admixture fractionTDTVTATN0TNDTXYTXYNDS0.1 0.5 12 3Million years agoNN1NN0NNDNXYNXYNDNS0 10,000 20,000 30,000 40,000 50,000P opulation size Fig. 4. Model-averaged parameter estimates with 95% confidence intervals estimated by moving blocks bootstrap (35). Key: ma, fraction of Y introgressed from N; mb, fraction of D introgressed from S; mg, fraction of N introgressed from XY; md, fraction of ND introgressed from S ; TXYNDS, superarchaic separation time; TXY, separation time of X and Y; TND, separation time of N and D; TN0, end of early epoch of Neanderthal history; TA, age of Altai Neanderthal fossil; TV, age of Vindija Neanderthal fossil; TD, age of Denisovan fossil; NS, size of superarchaic population; NXYND, size of populations XYND and XYNDS; NXY, size of population XY; NND, size of population ND; NN0, size of early Neanderthal population; NN1, size of late Neanderthal population. Parameters that exist in only one model are not averaged.
5 of 7It seems likely that superarchaics descend from the initial human settlement of Eurasia. As discussed above, the large effective size of the superarchaic population hints that it comprised at least two deeply divided subpopulations, of which one mixed with neandersovans and another with Denisovans. We suggest that around 700 ka ago, neandersovans expanded from Africa into Eurasia, endured a bottleneck of population size, interbred with indigenous Eurasians, largely replaced them, and separated into eastern and western subpopulationsâ€”Denisovans and Neanderthals. These same events unfolded once again around 50 ka ago as modern humans expanded out of Africa and into Eurasia, largely replacing the Neanderthals and Denisovans.MATERIALS AND METHODSStudy design Our sample of modern genomes includes Europeans but not other Eurasians. This allowed us to avoid modeling gene flow from Denisovans because there is no evidence of such gene flow into Europeans. The precision of our estimates depends largely on the number of nucleotides studied. For this reason, we used entire high-coverage genomes. The number of genomes sampled per population has little effect on our analyses, because of our focus on the between-population component of genetic variation, i.e., on site pattern frequencies. Nonetheless, our sample of modern genomes for the Yoruban, French, and English includes all those available from the Simons Genome Diversity Project (SGDP) (14), as detailed in the Supplementary Materials. We also included all available high-coverage archaic genomes (11). These data provide extremely accurate estimates of site pattern frequencies, as indicated by the tiny confidence intervals in Fig.2. The large confi dence intervals for some parameters in Fig.4 reflect identifiability problems (discussed below) and would not be alleviated by an increase in sample size. Quality control Our quality control (QC) pipeline for the SGDP genomes excludes genotypes at which an FL value equals 0 or N. We also excluded sex chromosomes, normalized all variants at a given nucleotide site using the human reference genome, excluded sites within seven bases of the nearest insertion-deletion, and included sites only if they were monomorphic or were biallelic single-nucleotide polymorphisms. Further details are provided in the Supplementary Materials. All ancient genomes were also filtered against .bed files, which identify bases that pass the Max Planck QC filters. These .bed files are available at http://ftp.eva.mpg.de/neandertal/Vindija/FilterBed. Molecular clock calibration We assumed a mutation rate of 1.110 per site per generation (34) and a generation time of 29 yearsâ€”a yearly rate of 0.3810. To calibrate the molecular clock, we assumed that the modern and neandersovan lineages separated TXYND=25,920 generations before the present ( 12). This is based on an average of several estimates published by Prfer etal. [(2), table S12.2]. The average of their estimates is 570.25 ka, assuming a mutation rate 0.510/base pair/year. Under our clock, their separation time becomes 751.69 ka or 25,920 generations. Statistical analysis Because of our focus on deep history, we based statistical analyses on site pattern frequencies, using the Legofit statistical package ( 10). This method ignores the within-population component of genetic variation and is therefore unaffected by recent changes in population size. For example, the sizes of populations X, Y, and D (Fig.1) have no effect, so we need not complicate our model with parameters describing the size histories of these populations. This allows us to focus on the distant past. Nonetheless, our models are quite complex. For example, model abgd has 17 free parameters. To choose among models of this complexity, we need methods of residual analysis, model selection, and model averaging. Legofit provides these methods, but alternative methods generally do not. These methods are described in detail elsewhere (10), so we summarize them only briefly here. We chose among models by minimizing the bepe (22, 23). This approach was needed because we could not use methods, such as Akaikeâ€™s information criterion (21), that depend on likelihood. Bepe is analogous to cross validation but uses bootstrap replicates instead of partitions of the data. The model is fit to each bootstrap replicate and then tested against the real data, after applying a correction for bootstrap bias. Bepe estimates the mean squared difference between observed and predicted site pattern frequencies, when the model is fit to one dataset and tested against another. We also used booma (24), which assigns weights to individual models, based on their bepe values. Parameters are estimated as the weighted average of estimates from individual models. The booma weight of the i th model is the fraction of replicates (including the real data and 50 bootstrap replicates) in which that model wins, i.e., has the lowest value of bepe. Because bootstrap replicates approximate repeated sampling from the process that generated the data, a model will receive zero weight if its disadvantage (as measured by bepe) is large compared with variation in repeated sampling. Figure S3 illustrates a problem of statistical identifiability. Several parameters are tightly correlated with others, indicating that our problem has fewer dimensions than parameters. This does not lead to incorrect estimates, but it broadens the confidence intervals of the parameters involved. Legofit addresses this problem using principal components analysis to remove dimensions that account for less than a fraction 0.001 of the total variance. This narrows confidence intervals and increases the accuracy of parameter estimates. Uncertainties are estimated by moving blocks bootstrap (35), using a block size of 500 single-nucleotide polymorphisms. Our statistical pipeline is detailed in the Supplementary Materials. SUPPLEMENTARY MATERIALSSupplementary material for this article is available at http://advances.sciencemag.org/cgi/ content/full/6/8/eaay5483/DC1 Supplementary Materials and Methods Fig. S1. Heterozygosity as a function of FL value for genome SS6004468 of the SGDP (14). Fig. S2. Observed site pattern frequencies excluding the Vindija genome. Fig. S3. Associations between estimates of several pairs of parameters after second stage in analysis of model abgd. Table S1. Model-averaged parameter estimates. Table S2. Estimates under model a. Table S3. Estimates under model abgd with a data set that excludes the Vindija Neanderthal genome. References (36â€“39) View/request a protocol for this paper from Bio-protocol.REFERENCES AND NOTES 1. R. E. Green, A. S. Malaspinas, J. Krause, A. W. Briggs, P. L. F. Johnson, C. Uhler, M. Meyer, J. M. Good, T. Maricic, U. Stenzel, K. Prfer, M. Siebauer, H. A. Burbano, M. Ronan, J. M. Rothberg, M. Egholm, P. Rudan, D. Brajkovi, . Kuan, I. Gui, M. Wikstrm,
6 of 7L. Laakkonen, J. Kelso, M. Slatkin, S. Pbo, A complete Neandertal mitochondrial genome sequence determined by high-throughput sequencing. Cell 134, 416 (2008). 2. K. Prfer, F. Racimo, N. Patterson, F. Jay, S. Sankararaman, S. Sawyer, A. Heinze, G. Renaud, P. H. Sudmant, C. de Filippo, H. Li, S. Mallick, M. Dannemann, Q. Fu, M. Kircher, M. Kuhlwilm, M. Lachmann, M. Meyer, M. Ongyerth, M. Siebauer, C. Theunert, A. Tandon, P. Moorjani, J. Pickrell, J. C. Mullikin, S. H. Vohr, R. E. Green, I. Hellmann, P. L. F. Johnson, H. Blanche, H. Cann, J. O. Kitzman, J. Shendure, E. E. Eichler, E. S. Lein, T. E. Bakken, L. V. Golovanova, V. B. Doronichev, M. V. Shunkov, A. P. Derevianko, B. Viola, M. Slatkin, D. Reich, J. Kelso, S. Pbo, The complete genome sequence of a Neanderthal from the Altai Mountains. Nature 505, 43 (2014). 3. M. Meyer, M. Kircher, M.-T. Gansauge, H. Li, F. Racimo, S. Mallick, J. G. Schraiber, F. Jay, K. Prfer, C. de Filippo, P. H. Sudmant, C. Alkan, Q. Fu, R. Do, N. Rohland, A. Tandon, M. Siebauer, R. E. Green, K. Bryc, A. W. Briggs, U. Stenzel, J. Dabney, J. Shendure, J. Kitzman, M. F. Hammer, M. V. Shunkov, A. P. Derevianko, N. Patterson, A. M. Andrs, E. E. Eichler, M. Slatkin, D. Reich, J. Kelso, S. Pbo, A high-coverage genome sequence from an archaic Denisovan individual. Science 338, 222 (2012). 4. R. G. Klein, Anatomy, behavior, and modern human origins. J. World Prehist. 9, 167 (1995). 5. R. G. Klein, The Human Career: Human Biological and Cultural Origins (University of Chicago Press, ed. 3, 2009). 6. J. J. Hublin, The origin of Neandertals. Proc. Natl. Acad. Sci. U.S.A. 106, 16022 (2009). 7. J. J. Hublin, in Neandertals and Modern Humans in Western Asia, T. Akazawa, K. Aoki, O. Bar-Yosef, Eds. (Kluwer, 1998), pp. 295. 8. H. Li, S. Mallick, D. Reich, Population size changes and split times, Supplementary Information 12 of Prfer et al. (2) (2014). 9. M. Meyer, J.-L. Arsuaga, C. de Filippo, S. Nagel, A. Aximu-Petri, B. Nickel, I. Martnez, A. Gracia, J. M. Bermdez de Castro, E. Carbonell, B. Viola, J. Kelso, K. Prfer, S. Pbo, Nuclear DNA sequences from the Middle Pleistocene Sima de los Huesos hominins. Nature 531, 504 (2016). 10. A. R. Rogers, Legofit: Estimating population history from genetic data. bioRxiv 613067 [Preprint]. 18 April 2019. https://doi.org/10.1101/613067. 11. K. Prfer, C. de Filippo, S. Grote, F. Mafessoni, P. Korlevi, M. Hajdinjak, B. Vernot, L. Skov, P. Hsieh, S. Peyrgne, D. Reher, C. Hopfe, S. Nagel, T. Maricic, Q. Fu, C. Theunert, R. Rogers, P. Skoglund, M. Chintalapati, M. Dannemann, B. J. Nelson, F. M. Key, P. Rudan, . Kuan, I. Gui, L. V. Golovanova, V. B. Doronichev, N. Patterson, D. Reich, E. Eichler, M. Slatkin, M. H. Schierup, A. Andrs, J. Kelso, M. Meyer, S. Pbo, A high-coverage Neandertal genome from Vindija Cave in Croatia. Science 358, 655 (2017). 12. A. R. Rogers, R. J. Bohlender, C. D. Huff, Early history of Neanderthals and Denisovans. Proc. Natl. Acad. Sci. U.S.A. 114, 9859 (2017). 13. F. Mafessoni, K. Prfer, Better support for a small effective population size of Neandertals and a long shared history of Neandertals and Denisovans. Proc. Natl. Acad. Sci. U.S.A. 114, E10256â€“E10257 (2017). 14. S. Mallick, H. Li, M. Lipson, I. Mathieson, M. Gymrek, F. Racimo, M. Zhao, N. Chennagiri, S. Nordenfelt, A. Tandon, P. Skoglund, I. Lazaridis, S. Sankararaman, Q. Fu, N. Rohland, G. Renaud, Y. Erlich, T. Willems, C. Gallo, J. P. Spence, Y. S. Song, G. Poletti, F. Balloux, G. van Driem, P. de Knijff, I. G. Romero, A. R. Jha, D. M. Behar, C. M. Bravi, C. Capelli, T. Hervig, A. Moreno-Estrada, O. L. Posukh, E. Balanovska, O. Balanovsky, S. Karachanak-Yankova, H. Sahakyan, D. Toncheva, L. Yepiskoposyan, C. Tyler-Smith, Y. Xue, M. S. Abdullah, A. RuizLinares, C. M. Beall, A. Di Rienzo, C. Jeong, E. B. Starikovskaya, E. Metspalu, J. Parik, R. Villems, B. M. Henn, U. Hodoglugil, R. Mahley, A. Sajantila, G. Stamatoyannopoulos, J. T. Wee, R. Khusainova, E. Khusnutdinova, S. Litvinov, G. Ayodo, D. Comas, M. F. Hammer, T. Kivisild, W. Klitz, C. A. Winkler, D. Labuda, M. Bamshad, L. B. Jorde, S. A. Tishkoff, W. S. Watkins, M. Metspalu, S. Dryomov, R. Sukernik, L. Singh, K. Thangaraj, S. Pbo, J. Kelso, N. Patterson, D. Reich, The simons genome diversity project: 300 genomes from 142 diverse populations. Nature 538, 201 (2016). 15. G. S. Jacobs, G. Hudjashov, L. Saag, P. Kusuma, C. C. Darusallam, D. J. Lawson, M. Mondal, L. Pagani, F.-X. Ricaut, M. Stoneking, M. Metspalu, H. Sudoyo, J. S. Lansing, M. P. Cox, Multiple deeply divergent Denisovan ancestries in Papuans. Cell 177, 1010.e32 (2019). 16. A. R. Rogers, R. J. Bohlender, C. D. Huff, Reply to Mafessoni and Prfer: Inferences with and without singleton site patterns. Proc. Natl. Acad. Sci. U.S.A. 114, E10258â€“E10260 (2017). 17. P. J. Waddell, J. Ramos, X. Tan, Homo denisova, correspondence spectral analysis, finite sites reticulate hierarchical coalescent models and the Ron Jeremy hypothesis. arXiv:1112.6424 [q-bio.PE] (29 December 2011). 18. P. J. Waddell, Happy New Year Homo erectus? More evidence for interbreeding with archaics predating the modern human/Neanderthal split. arXiv:1312.7749 [q-bio.PE] (30 December 2013). 19. M. Kuhlwilm, I. Gronau, M. J. Hubisz, C. de Filippo, J. Prado-Martinez, M. Kircher, Q. Fu, H. A. Burbano, C. Lalueza-Fox, M. de la Rasilla, A. Rosas, P. Rudan, D. Brajkovic, . Kucan, I. Guic, T. Marques-Bonet, A. M. Andrs, B. Viola, S. Pbo, M. Meyer, A. Siepel, S. Castellano, Ancient gene flow from early modern humans into Eastern Neanderthals. Nature 530, 429 (2016). 20. R. Ferring, O. Oms, J. Agust, F. Berna, M. Nioradze, T. Shelia, M. Tappen, A. Vekua, D. Zhvania, D. Lordkipanidze, Earliest human occupations at Dmanisi (Georgian Caucasus) dated to 1.85.78 Ma. Proc. Natl. Acad. Sci. U.S.A. 108, 10432 (2011). 21. H. Akaike, A new look at the statistical model identification. IEEE Trans. Automat. Contr. 19, 716 (1974). 22. B. Efron, Estimating the error rate of a prediction rule: Improvement on cross-validation. J. Am. Stat. Assoc. 78, 316 (1983). 23. B. Efron, R. J. Tibshirani, An Introduction to the Bootstrap (Chapman and Hall, 1993). 24. S. T. Buckland, K. P. Burnham, N. H. Augustin, Model selection: An integral part of inference. Biometrics 53, 603 (1997). 25. A. Kong, M. L. Frigge, G. Masson, S. Besenbacher, P. Sulem, G. Magnusson, S. A. Gudjonsson, A. Sigurdsson, A. Jonasdottir, A. Jonasdottir, W. S. Wong, G. Sigurdsson, G. B. Walters, S. Steinberg, H. Helgason, G. Thorleifsson, D. F. Gudbjartsson, A. Helgason, O. T. Magnusson, U. Thorsteinsdottir, K. Stefansson, Rate of de novo mutations and the importance of fatherâ€™s age to disease risk. Nature 488, 471 (2012). 26. G. Amster, G. Sella, Life history effects on the molecular clock of autosomes and sex chromosomes. Proc. Natl. Acad. Sci. 113, 1588 (2016). 27. M. Nei, N. Takahata, Effective population size, genetic diversity, and coalescence time in subdivided populations. J. Mol. Evol. 37, 240 (1993). 28. G. P. Rightmire, Human evolution in the Middle Pleistocene: The role of Homo heidelbergensis. Evol. Anthropol. 6, 218 (1998). 29. R. Foley, M. M. Lahr, Mode 3 technologies and the evolution of modern humans. Camb. Archaeol. J. 7, 3 (1997). 30. P. Endicott, S. Y. W. Ho, C. Stringer, Using genetic evidence to evaluate four palaeoanthropological hypotheses for the timing of Neanderthal and modern human origins. J. Hum. Evol. 59, 87 (2010). 31. J. L. Bischoff, D. D. Shamp, A. Aramburu, J. L. Arsuaga, E. Carbonell, J. B. De Castro, The Sima de los Huesos hominids date to beyond U/Th equilibrium (> 350kyr) and perhaps to 400kyr: New radiometric dates. J. Archaeol. Sci. 30, 275 (2003). 32. J. L. Bischoff, R. W. Williams, R. J. Rosenbauer, A. Aramburu, J. L. Arsuaga, N. Garca, G. Cuenca-Bescs, High-resolution U-series dates from the Sima de los Huesos hominids yields 600 + kyrs: Implications for the evolution of the early Neanderthal lineage. J. Archaeol. Sci. 34, 763 (2007). 33. D. Reich, Who We Are and How We Got Here: Ancient DNA and the New Science of the Human Past (Pantheon Books, 2018). 34. K. R. Veeramah, M. F. Hammer, The impact of whole-genome sequencing on the reconstruction of human population history. Nat. Rev. Genet. 15, 149 (2014). 35. R. Y. Liu, K. Singh, Moving blocks jackknife and boostrap capture weak dependence, in Exploring the â€œLimitsâ€ of the Bootstrap, R. LePage, L. Billard, Eds. (Wiley, 1992), pp. 225. 36. Chimpanzee Sequencing and Analysis Consortium, Initial sequence of the chimpanzee genome and comparison with the human genome. Nature 437, 69 (2005). 37. A. Scally, J. Y. Dutheil, L. W. Hillier, G. E. Jordan, I. Goodhead, J. Herrero, A. Hobolth, T. Lappalainen, T. Mailund, T. Marques-Bonet, S. McCarthy, S. H. Montgomery, P. C. Schwalie, Y. A. Tang, M. C. Ward, Y. Xue, B. Yngvadottir, C. Alkan, L. N. Andersen, Q. Ayub, E. V. Ball, K. Beal, B. J. Bradley, Y. Chen, C. M. Clee, S. Fitzgerald, T. A. Graves, Y. Gu, P. Heath, A. Heger, E. Karakoc, A. Kolb-Kokocinski, G. K. Laird, G. Lunter, S. Meader, M. Mort, J. C. Mullikin, K. Munch, T. D. O'Connor, A. D. Phillips, J. Prado-Martinez, A. S. Rogers, S. Sajjadian, D. Schmidt, K. Shaw, J. T. Simpson, P. D. Stenson, D. J. Turner, L. Vigilant, A. J. Vilella, W. Whitener, B. Zhu, D. N. Cooper, P. de Jong, E. T. Dermitzakis, E. E. Eichler, P. Flicek, N. Goldman, N. I. Mundy, Z. Ning, D. T. Odom, C. P. Ponting, M. A. Quail, O. A. Ryder, S. M. Searle, W. C. Warren, R. K. Wilson, M. H. Schierup, J. Rogers, C. Tyler-Smith, R. Durbin, Insights into hominid evolution from the gorilla genome sequence. Nature 483, 169 (2012). 38. D. Gordon, J. Huddleston, M. J. Chaisson, C. M. Hill, Z. N. Kronenberg, K. M. Munson, M. Malig, A. Raja, I. Fiddes, L. W. Hillier, C. Dunn, C. Baker, J. Armstrong, M. Diekhans, B. Paten, J. Shendure, R. K. Wilson, D. Haussler, C.-S. Chin, E. E. Eichler, Long-read sequence assembly of the gorilla genome. Science 352, aae0344 (2016). 39. K. Price, R. M. Storn, J. A. Lampinen, Differential Evolution: A Practical Approach to Global Optimization (Springer Science and Business Media, 2006). Acknowledgments: We thank R. Bohlender, E. Cashdan, F. Mafessoni, N. Rogers, J. Seger, and T. Webster for comments. This project was declared exempt (IRB_00093972) from review by the Institutional Review Board of the University of Utah on 13 July 2016. Funding: This work was supported by NSF BCS 1638840 (A.R.R.), NSF GRF 1747505 (A.A.A.), and the Center for High Performance Computing at the University of Utah
7 of 7(A.R.R.). Author contributions: A.R.R. designed the study, did the statistical analyses, and wrote the paper. N.S.H. and A.A.A. developed and used the QC pipeline. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper, the Supplementary Materials, or at osf.io/vrwna. The Legofit software is available at https://github.com/alanrogers/legofit. Additional data related to this paper may be requested from the authors. Submitted 27 June 2019 Accepted 20 December 2019 Published 21 February 2020 10.1126/sciadv.aay5483 Citation: A. R. Rogers, N. S. Harris, A. A. Achenbach, Neanderthal-Denisovan ancestors interbred with a distantly related hominin. Sci. Adv. 6, eaay5483 (2020).