Neanderthals and humans interbred in Europe for 5,400 years

Around 600 kya, Europe was invaded by large-brained hominins using Acheulean stone tools (1, 2). They were probably African immigrants, because similar fossils and tools occur earlier in Africa. They have been called archaic Homo sapiens, Homo heidelbergensis, and early Neanderthals, yet they remain mysterious. They may have been ancestors of Neanderthals and modern humans (3), or ancestors of Neanderthals only (4, 5), or an evolutionary dead end. According to this last hypothesis, they were replaced later in the Middle Pleistocene by a wave of African immigrants that separated Neanderthals from modern humans and introduced the Levallois stone tool tradition to Europe (6, 7). To address this controversy, we introduce a statistical method and use it to study genetic data of Africans, Eurasians, Neanderthals, and Denisovans.

Our method extends an idea introduced by Reich et al. (8, 9). Their “ABBA-BABA” statistics infer admixture from the frequency with which derived alleles are shared by pairs of samples. As we have shown (10), these estimators have large biases when populations receive gene flow from more than one source. The magnitudes of these biases depend on the sizes and separation times of ancestral populations. Our method avoids bias by estimating these parameters simultaneously.

To accomplish this, our method uses an expanded dataset. ABBA-BABA statistics summarize allele sharing by pairs of samples. We extend this approach to include larger subsets, such as trios of samples, and to use all available subsets. This opens a rich and heretofore unused window into population history.


Fig. 1.
(A) Population tree representing an African population, XX; a Eurasian population, YY; Neanderthals, NN; and Denisovans, DD. The model involves admixture, mNmN; time parameters, TiTi; and population sizes, NiNi. (B) Population tree with embedded gene tree. A mutation on the solid red branch would generate site pattern ynyn (shown in red at the base of the tree). One on the solid blue branch would generate yndynd. Mutations on the dashed black branches would be ignored. “0” and “1” represent the ancestral and derived alleles.

Although our method can accommodate complex models, we work here with a four-population model of history (Fig. 1A), which has broad empirical support (11, 12). In this model, Neanderthals (NN) contribute genes to Eurasians (YY) but not to Africans (XX). The model allows no gene flow from Denisovans (DD), for reasons explained below. Combinations of uppercase letters, such as NDND, refer to the population ancestral to NN and DD. Lowercase letters, such as nn and dd, refer to individual haploid genomes sampled from these populations.

The gene tree describes how genes coalesce within the tree of populations. Fig. 1B illustrates one of many possible gene trees. Although closely linked nucleotide sites tend to share the same gene tree, this is not the case for sites farther apart on the chromosome, and any set of autosomal sequence data will encompass a multitude of gene trees.

The gene tree determines opportunities for allele sharing among samples. For example, a mutation on the solid red branch in Fig. 1B would be present in yy and nn but absent in xx and dd. We refer to this as the “ynyn site pattern.” Similarly, a mutation on the solid blue branch would generate site pattern yndynd. In a four-population model, there are 10 polymorphic site patterns, excluding singletons. We can tabulate their frequencies in sequence data and calculate their probabilities given particular population histories. Our program, legofit (described in Section S1), estimates parameters by fitting observed to expected frequencies. Whereas ABBA-BABA statistics use only 2 site patterns (“ABBA” and “BABA”), legofit uses all 10. This allows it to estimate additional parameters and avoid the biases discussed above.

Fig. 3.
Covariation of estimates of mNmN and 2NN2NN across bootstrap replicates. Data are as in Fig. 2.

All four analyses—YRI.CEU, YRI.CHB, LWK.CEU, and LWK.CHB—yield similar results. Estimates of Neanderthal admixture (mNmN) and Neanderthal–Denisovan separation time (TNDTND) appear in Fig. 4. The admixture estimates are 1–3%, in broad agreement with previous results. Our results do not, however, support the view that East Asians carry more Neanderthal DNA than Europeans (12, 14, 17⇓⇓⇓–21). This view may be an artifact of ascertainment bias (17) or of the biases documented by Rogers and Bohlender (10). On the other hand, the East Asian excess may be real, but hidden by the broad confidence intervals surrounding our estimates of mNmN.



Fig. 4.
Estimates of Neanderthal admixture (mNmN) and the Neanderthal–Denisovan separation time (TNDTND). The vertical line (Lower) shows TXYNDTXYND. Horizontal lines show 95% confidence intervals based on 50 moving-blocks bootstrap replicates. All point estimates and confidence intervals are based on stage 2 of the analysis.

All estimates of TNDTND, the separation time of Neanderthals and Denisovans, are close to 25,600 generations ago—only about 300 generations after the separation of archaics from moderns. Furthermore, this separation time is estimated with high confidence, judging from the narrow confidence intervals in Fig. 4, Lower. During the interval between the two separation events, the ancestral archaic population was apparently very small. Our point estimates of 2NND2NND range from about 100 to about 1,000, with narrow confidence intervals. Following the Neanderthal–Denisovan separation, our results imply a relatively large Neanderthal population, with 2N2N in the tens of thousands. Fig. S3 graphs the history of effective population size of Neanderthals, moderns, and their ancestors, as implied by the YRI.CEU analysis.