Archaic human ancestry in East Asia

We performed principal component analysis (PCA) (44) by defining the first two principal components (PCs) using the Denisova, the Neandertal, and the chimpanzee and projected extant humans on the resulting axes of variation (24). This setup resulted in PC1 describing general genetic similarity to archaic humans (represented by both the Neandertal and Denisova genomes) and PC2 contrasting genetic similarity between Neandertal and Denisova. Under the assumption of a common shared history between Neandertal and Denisova (24) as well as no admixture between archaic populations and the ancestors of extant human populations since the diversification of modern humans (Fig. 1A), extant individuals are expected to be homogeneously distributed between archaic human and chimpanzee variations (24) (SI Materials and Methods).

We recovered the previously reported (24) pattern—that extant human variation is largely organized in three clusters corresponding to Africans, Oceanians, and other non-Africans, respectively (Fig. 1B and Fig. S1). In the worldwide collection of extant populations, we used Procrustes superimposition (45) to compare PC1 and PC2 with geographic coordinates (longitude and latitude) for each population and found that sampling location was mirrored, to some extent, by archaic ancestry (individuals: Procrustes correlation = 0.127, P << 10⁻⁶; population means: Procrustes correlation = 0.309, P < 0.01). This pattern was expected given a model with two episodes of archaic gene flow involving non-Africans and Oceanians (Fig. 2B) (23, 24). However, PC1 and PC2 were also correlated with geography in a region comprising Eurasia and the Americas (individuals: Procrustes correlation = 0.104, P < 10⁻⁴; population means: Procrustes correlation = 0.335, P = 0.017), which seemed incompatible with previously suggested admixture scenarios postulating that archaic ancestry is homogeneous in people of European, Asian, and Native American descent (23, 24).

By comparing the geographic locations of individuals from different regions with the archaic ancestry signal (PC1 values), East Asian (two-sided t test: P < 10⁻⁶) and Native American (P < 10⁻⁸) populations were found to be more similar to archaic hominins compared with European and Central/South Asian populations (Fig. 1B and Table S1). In addition, the archaic ancestry signal was positively correlated with distance from East Africa (individuals: Spearman's rank correlation coefficient r_s = 0.189, P < 10⁻⁹; population means: r_s = 0.600, P < 10⁻⁴) (Fig. 1C). Because it is well-known that PCs of modern human genetic variation capture geography to some extent (35, 46), we also computed PCs using modern human genetic variation in Eurasia (SI Results) and found that the top two PCs of differentiation between Europe and East Asia were correlated with the archaic ancestry signal (PC1: r_s = 0.138, P < 10⁻⁵; PC2: r_s = 0.090, P = 0.002). This result suggests that, if separated along the major axes of differentiation between Europe and East Asia, individuals on the eastern end of the spectrum tend to be more similar to archaic human genomes. In contrast, after correction for multiple tests, we found no correlations between the archaic ancestry signal and intraregional PCs within Europe, and we did not find patterns of variation in archaic ancestry within Africa and America that could not be explained by more recent admixture with non-African populations (SI Results and Table S2).

To disentangle the signal of Denisova admixture from the signal of Neanderthal admixture, we investigated the distribution of PC2, which separated Denisovans and Neandertals, and found that proximity to Neandertal increases with distance from Africa (individuals: r_s = 0.078, P = 0.010). However, East Asians were an exception to this trend in that they were significantly closer to Denisovans compared with Europeans (P = 0.003), all West/Central Eurasians (P = 0.006), and Americans (P < 10⁻⁴). In a spatial interpolation of log-transformed PC2 values, the affinity to the Denisova genome seemed to be strongest in Southern China and Southeast Asia (Fig. 1D), which is noteworthy considering the supposed East Eurasian distribution of Denisovans (24). To investigate this geographical pattern using an alternative approach, we identified SNPs in chimpanzees and Neandertals that shared the ancestral allele, and Denisova had the derived allele; we computed the frequency of the derived Denisova allele in global modern human populations. Using the same method of spatial interpolation, we found that East Asian populations, particularly Southeast Asian populations, had, on average, a greater frequency of the derived Denisova allele compared with other populations (except for Oceanians) (Fig. 1E). For example, although the greatest Denisova allele frequency was found in Papuans (53.5%), Yizu from Southern China had a greater frequency of the Denisova allele than Melanesians from Bougainville (53.0% vs. 52.9%) (Table S3).

Although the above described spatial patterns of variation in archaic ancestry among extant humans suggest that the distribution of archaic ancestry is more complex than previously suggested, it is also well-known that population differentiation, genetic diversity, and frequencies of derived alleles correlate with distance from Africa (1, 3, 34, 35), which is in line with anatomically modern human expansion out from Africa (3, 21). To investigate the impact of demography on signals of archaic admixture, we simulated a model of human expansion out of Africa, which is similar to the approach used in ref. 21. Briefly, the model comprises serial founder events with successive bottlenecks initiated ∼51,000 years ago (kya) and leaves 100 descendant populations (colonies). We also model two archaic hominin populations that may (or may not) contribute genetic material to the extant populations (Materials and Methods and Fig. 2 A–C). Based on simulated data from the model, we projected individuals from the colony populations onto PCs defined by a hypothetical chimpanzee and one sample from each archaic population.

Assuming no admixture with archaic populations, we obtain no distinct clustering of different colony populations (Fig. 2 A, D, and G). In contrast, if we assume two separate admixture events involving each archaic population, akin to the conclusions of the works in refs. 23 and 24, we observe similar clustering patterns as in the empirical data but without regional intracluster patterns (Fig. 2 B and E). However, if we filter out rare alleles to mimic some aspects of ascertainment bias affecting SNP-chip data, there is a strong correlation between PC1/PC2 and colony number among the hypothetical Eurasian colonies—the colonies that have only been involved in one episode of admixture (Fig. 2H). Thus, this model predicts increasing affinity to the hypothetical Neandertal with increasing distance from the founder population, which is caused by the combination of ascertainment bias and genetic drift. This observation is in line with the correlations of the signal of archaic ancestry and distance from Africa that we reported above but in stark contrast to the observation that East Asians are significantly closer to Denisova relative to Neandertal. Interestingly, when we add a third archaic admixture event representing a Denisova-related contribution to the ancestral population of East Asians, we obtain a pattern that is qualitatively more similar to the empirical data (Fig. 2 C, F, and I).

To test the hypothesis of Denisovan ancestry in East Asians using an alternative approach, we performed 4-population tests (23, 24, 40, 47, 48) on diploid genotype data from different regions. Note that using population allele frequencies is confounded by demographic effects (48) but may also provide more power and a more diverse set of samples compared with using low-coverage shotgun sequences from a smaller set of individuals (23, 24). To investigate fine-scale geographical patterns suggested by the PC analyses (Fig. 1 D and E), we divided the East Asian populations in our dataset into two groups, South and North of Beijing (with the North subgroup including Beijing and Japan). We then tested the hypothesis of no gene flow in the population topology (Denisova, chimpanzee), (Northeast Asians, Southeast Asians) and found that allele frequencies in the Southeastern Asian group were significantly more similar to the Denisova data (D = 0.55 ± 0.23% SE, Z = 2.40) (Table S4). In addition, we found that allele frequencies in the Southeastern Asian group were significantly more similar to the Denisova genome compared with the Neandertal genome (D = 0.66 ± 0.30%, Z = 2.22). This differentiation is not predicted by the effects of ascertainment bias, which was observed in our simulations, because ascertainment bias is expected to artificially increase similarity to Neandertal under the model suggested by previous studies (23, 24). Indeed, the only other pairwise regional comparisons for which we see a strong pattern of similarity to Denisova compared with both chimpanzee and Neandertal include Oceanians (Fig. 3 and Table S4).

The reason why an affinity to Denisova is not as clearly seen in 4-population tests between Southeast Asians and other populations (except Northeast Asians) (Fig. 3 and Table S4) could be the joint effect of ascertainment bias and genetic drift, which is indicated by our simulations. For example, a test of the topology [archaic population A, archaic population B, (colonies 50–70), (colonies 97–100)] fails to detect a fraction of 5% ancestry contributed by archaic population B to colonies 97–100 in simulation model ii (Fig. 2B), with filtering for minor allele frequency > 5%, and in fact, it suggests a skew toward colonies 97–100 being more similar to archaic population A (D = 1.0 ± 3.7%, Z = 0.27). In general, a signal of archaic ancestry is detected only if tests are performed between populations in which the magnitude of genetic drift has been similar. For example, if colonies 50–70 are substituted with colonies 90–96 in the test above, the test statistic significantly deviates from zero, and the admixture is detected (D = −7.0 ± 2.7%, Z = −2.6), although colonies 50–70 have exactly the same true fraction of archaic ancestry as colonies 90–96. However, the admixture is detected in both cases (Z < −8) if the data are unascertained (no minor allele frequency filtering).

To circumvent the effects of ascertainment bias, we identified polymorphisms between complete genomes sequenced to high coverage from seven individuals of Korean, Chinese, European, and African ancestry (49–52) and retained positions that overlapped with data from Neandertal and Denisova. We then tested the hypothesis of no gene flow in the population topology (Denisova, Neandertal), [East Asian, (European or African)] using haploid sets of SNPs from one individual from each population (23, 24, 48). Similar to previous analysis using this approach (23, 24), we did not observe any significant deviations from the null hypothesis in tests between Europeans and East Asians (Table S5). However, we note that this approach offers less power compared to using multiple individuals from each population (SEs were 3.4–4.4 times larger than in the population-based test between Southeast and Northeast Asia) (23, 24), and neither of the two complete East Asian genomes were from Southeast Asia, the region where we observe the strongest signal of Denisova ancestry.

We obtained haploid autosomal variants from the Neandertal (23) and Denisova (24) genomes and removed bases with quality <40 and from within 5 and 1 bp from the 5′ end of the sequence reads from Neandertal and Denisova, respectively (24). We filtered reads with mapping quality <90 (Neandertal) and <37 (Denisova) (24) and randomly chose a single read from positions covered by multiple reads. We obtained phased HapMap 3 genotypes from 630 individuals (41, 42) and phased genotypes for 938 globally distributed individuals from the Human Genome Diversity Project (35). We intersected this dataset with chimpanzee (43), Neandertal, and Denisova genotypes, resulting in 228,984 overlapping SNPs. To avoid the effect of postmortem nucleotide misincorporations in analyses that included ancient genomes, we followed the information in refs. 23 and 24 by removing all transition SNPs (C/T and G/A), resulting in 40,656 SNPs.

From the seven complete genome sequences, we identified SNPs using genotypes called by the 1,000 Genomes Project (49). We excluded hypermutable CpG sites and transition substitutions and retained SNPs for which there was both Neandertal and Denisova genome data, using the same genotyping criteria for the ancient genomes as described above.

For PC analyses, we used a dataset created by randomly sampling one allele from each individual at each position to allow comparison with the single-pass ancient data. PCA was performed with EIGENSOFT 4.0 (44) on a dataset where one SNP from each pair of SNPs with r² > 0.2 had been excluded. In the dataset that included only transversion SNPs, the number of SNPs used was 38,848, but the full dataset comprised 183,166 SNPs. We first performed PCA using chimpanzee, Denisova, and Neandertal, and thereafter, projected all other samples on the resulting components. In the transversion analysis, (PC1, PC2) loadings were chimpanzee (−0.81, 0.072), Denisova (0.34, −0.74), and Neandertal: (0.47, 0.67).

For Procrustes analyses, two-sided t tests, and Pearson correlation analyses on the PCA results, we excluded Middle Easterners and Balochi because of possible gene flow from Africa. Procrustes analysis was performed using the vegan R package (65), and significance was assessed with 1,000,000 permutations. All tests were performed on raw PC scores unless explicitly stated otherwise.

We used the computer program ms (66) to simulate a serial founder model as in ref. 21, with the main modification being inclusion of two archaic human populations A and B (hypothetical Neandertal and Denisova, respectively). We assumed that the two archaic populations diverged from each other 200 kya and had a divergence time of the hypothetical Neandertal/Denisova ancestral population and modern humans of 350 kya (18, 23, 24). Outgroup (hypothetical chimpanzee) population divergence was assumed to be 4 million years (y). We assumed a generation time of 25 y. We sampled 10 haploid samples from each of 100 extant populations. The age was set to ∼50 kya for both the archaic samples by moving their lineages to an isolated population at that specific time point. This movement prevents the lineages from coalescing with other lineages in the simulation, and any private mutations do not contribute to any tests for admixture under an infinite sites mutation model. We simulated three scenarios: (i) no admixture with archaic populations, (ii) admixture fraction from archaic population A of 2.5% in colony 25 and admixture fraction of 5% from archaic population B in colony 97, and (iii) admixture fractions as in scenario ii with the addition of a 1% admixture fraction from archaic population B in colony 75. We based our analyses on two different datasets, each consisting of 38,848 independent SNPs. In the first dataset, only singleton SNPs were excluded (because they are not informative in our analyses). In the second dataset, SNPs with minor allele frequency > 5% were excluded (to investigate the effect of ascertainment bias).

4-population tests on allele frequencies (23, 24, 40, 47, 48) were implemented as in ref. 24 and performed using diploid genotype data from Africa, Europe, Oceania, America, the Middle East, Southeast Asia, and Northeast Asia. We computed SEs using a block jackknife, dropping 1 of 114 contiguous blocks with 600 SNPs in each block. Tests on complete genome sequence data (23, 24, 48) used a randomly sampled haploid copy from each individual, and SEs were computed using a block jackknife over contiguous blocks with 200 informative SNPs (ABBA or ABAB) in each block. Tests on the simulated data were performed by computing a block jackknife SE over 129 contiguous windows of 600 SNPs each. We used computed SEs to obtain Z scores for the test statistic D and interpreted |Z| > 2.0 as a statistically significant deviation from zero.