Ayta Negritos: Highest Levels of Denisovan DNA

Sequencing Analyses of Ayta Magbukon Negritos
We have sequenced and analyzed the genomes of five Ayta Magbukon individuals and present the results in Figure 4. In order to omit potential genome-sequence-processing biases, we started from the raw sequence data (FASTQs) for all individuals (including the comparative human sequence data), and processed all the sequence data in one batch, using standard mapping and filtering criteria (see Method details section for a full description).



D and f4-ratio statistics on sequenced data
Denisovan ancestry was estimated using the test D(Outgroup;Denisovan,EastAsian,X) or the f4-ratio statistic f4(Outgroup; Neanderthal,EastAsian,X) / f4(Outgroup;Neanderthal, EastAsian,Denisovan), with Altai and Vindija individuals representing the Neanderthals, and Kankanaey Cordilleran representing East Asians. Ayta Magbukon consistently demonstrate to have the highest level of Denisovan ancestry (Figure 4A; Data S2A and S2B). To assess whether Ayta Magbukon possess higher Denisovan ancestry than highland Papuans or Australians, we implemented the test D(Outgroup;Denisovan,Papuan/Australian,AytaMagbukon), with Chimp, Juhoansi, Biaka, Khomani San, or Yoruba used as outgroups. We confirm that the Ayta Magbukon possess higher levels of Denisovan ancestry relative to that of Australians or Papuans, with D values of 0.0168 - 0.0256 and all Z score of > 3 ((Figures S6A and S6B; Data S2C). To assess the percentage difference of Denisovan ancestry in Ayta Magbukon versus Australians or Papuans, we utilized the statistic D(Mbuti,Denisovan,EastAsian,AytaMagbukon) / D(Mbuti,Denisovan,East Asian,Papuan/Australian). We find that the Ayta Magbukon possess 29%–39% more Denisovan ancestry than do Australians and Papuans (Figure 4B; Data S2A).

S’ method
For identifying archaic introgressed sequences without the use of reference archaic genomes, we utilized the recently developed S-prime (S’) method15. The method is designed to detect divergent haplotypes that are in strong LD and that are of very low frequency or absent in outgroup African population (who are known not to harbor archaic ancestry). The S’ method was demonstrated to have better power relative to other archaic reference-free methods, and is appropriate for large-scale genomic data15. For our analysis, we utilized 15 SGDP Sub-Saharan African individuals as the outgroup population (which includes the Juhoansi, Mbuti, Khomani, Yoruba, and Biaka individuals), and 5 newly sequenced Ayta Magbukon and 15 SGDP Papuan individuals as the test populations.

Following retrieval of archaic sequences using S’, we leveraged the sequences to the reference Altai Denisovan and Vindija Neanderthal genomes. First, using the masks downloaded from cdna.eva.mpg.de/neandertal/, we filtered out the following sites from the archaic genomes: indels, sites within tandem repeats, sites with poor mapability, depth of coverage of < 10x, or mapping quality of < 25. For each filtered site, we compared the genomes of the test individuals to the reference Denisovan and Neanderthal genomes, and report if they match or mismatch. The match rate per archaic segment was then estimated by dividing the total number of matches over the total number of sites compared. Contour densities of the Denisovan versus Neanderthal match rates were plotted using the kde2d function of the R v3.3.1 MASS package. We find that both Ayta Magbukon and Papuans display similar patterns of affinity to the reference Denisovan and Neanderthal genomes, moderate affinity or match rate of ∼50% to the Altai Denisovan and high affinity or match rate of ∼80% to the Vindija Neanderthal (Figures 5A and 5B).

Following a previous approach15, we identify the segments as putatively Denisovan with the following criteria: presence of at least 30 introgressed alleles that can be compared with the reference archaic genomes; less than 30% match rate to the Vindija/Altai Neanderthal genome, and greater than 30% match rate to the Denisovan genome. Based on this criteria, we find that the Ayta Magbukon exhibit significantly higher magnitude of Denisovan ancestry relative to that of Papuans. The average amount of Denisovan sequence detected in Ayta Magbukon and Papuans are 100.38 Mb (95% CI: 89.66 – 111.11 Mb) and 79.35 Mb (95% CI: 76.89 – 81.79 Mb), respectively (p = 3.0x10-4); a difference of ∼27%.

To classify the introgressed haplotypes into Denisovan versus Neanderthal source, we applied the criteria used in a previous work15. We classified segments as putatively Neanderthal with the following conditions: presence of at least at least 30 introgressed alleles that can be compared with the reference archaic genomes; greater than 60% match rate to the Vindija Neanderthal genome; and less than 40% match rate to the Denisovan genome. We then identified segments as putatively Denisovan with the following conditions: presence of at least at least 30 introgressed alleles that can be compared with the reference archaic genomes; greater than 40% match rate to the Denisovan genome; and less than 30% match rate to the Vindija Neanderthal genome. Any segment that indistinguishable from Neanderthal or Denisovan are categorized as ambiguous.
The summary of archaic sequences per individual is presented in Figure S6C. Using the aforementioned criteria that distinguishes between different archaic sources, Ayta Magbukon and Papuans display similar amounts of Neanderthal sequences [(45.35 Mb (95% CI: 41.71 – 48.99 Mb) versus 46.47 Mb (95%CI: 44.96 – 47.97 Mb)] and ambiguous sequences [(7.52 Mb (95% CI: 5.47 – 9.57 Mb) versus 6.92 Mb (95%CI: 6.25 – 7.60 Mb)]. Consistently, we find Ayta Magbukon to present with significantly higher levels of Denisovan ancestry relative to Papuans (Figure 4D). The average amount of Denisovan sequences detected in Ayta Magbukon and Papuans are 51.94 Mb (95% CI: 44.62.66 – 59.25 Mb) and 41.96 Mb (95% CI: 39.54 – 44.37 Mb), respectively (p = 3.7x10-3); a difference of ∼24%.

Correlation Between Denisovan and Australasian ancestry
We plotted the correlation between Denisovan and Australasian-related ancestry among Negrito and Papuan-related groups. Denisovan ancestry was estimated using f4-ratio statistics, while Australasian-related ancestry was estimated using qpAdm. Consistent with the correlation analysis presented earlier, Negritos fall outside the slope formed by the Papuan-related populations. Furthermore, an important aspect of this study is the availability of Negritos with varying levels of admixture with East Asian-related populations. This allowed us for the first time to plot and extrapolate the level of Denisovan ancestry in an ‘unadmixed’ ancestral Negrito population (D = 4.47), which demonstrably have at least 46% more Denisovan ancestry than the highlander Papuans (Figure 3D). In addition, our analysis also demonstrated for the first time that Negritos have significantly higher levels of ‘Denisovan richness’ in their Australasian ancestry than Oceanians, which we define as the ratio of Denisovan ancestry over Australasian ancestry (Figure 3E).

Altogether, the evidence presented here, that Denisovan ancestry in Negritos is not associated with Papuan-related ancestry and that Negritos have higher levels of Denisovan richness, provide support to a demographic model where Negritos receive an additional and/or an independent Denisovan introgression event.

Evaluating Recent Admixture in Philippine Negritos versus Papuans
The lower levels of Denisovan ancestry in Papuans could be explained by dilution via recent admixture with non-Papuan populations with no detectable or low levels of Denisovan ancestry. To determine the presence of gene flow from populations with little or no Denisovan ancestry into Papuans or Philippine Negritos, we implemented the test D(Outgroup;X,Papuan, Philippine Negrito); with Balangao, Han, or Dai as a surrogate for East Asian; Even and Yakut as a surrogate for Siberian or North Asian; Karitiana or Mixe as a surrogate for Native American; Irula as a surrogate for South Asian; French or Sardinian as a surrogate for European; and Luhya or Mbuti as a surrogate for African (Figures S7A–S7L). Based on the various combination of D tests listed above, we do not find any population that shares more alleles with Papuans relative to Negritos (Figures S7A–S7L). In fact, it is the Negritos that consistently displayed evidence for recent admixture, given the high levels of allele sharing between Negritos and East Asians, relative to Papuans.

Admixture Graph Analysis
We constructed admixture graphs using the qpGraph tool of the AdmixTools v5.0 software package to determine which population-history model best fits the data that accounts for the number of Denisovan introgression events into Philippine Negritos and Papuans. Using the Phil_1KGP_SGDP_Ancient_Transv_317K dataset, we incorporarted in the model the chimpanzee as an outgroup, French to represent West Eurasians, Balangao Cordilleran and Liangdao 2 to represent East Asians, and Australian, Papuan, and Ayta Magbukon to represent Australasians (Figures 6A–6C). The qpGraph tool was ran with the following parameters: outpop (NULL), blgsize (0.05), diag (0.0001), hires (YES), and lsqmode (YES). A model is rejected when the worst fitting f statistical test has a significant Z score of > | 3 |. A model with a single Denisovan introgression event shared between Philippine Negritos and Australopapuans is rejected (worst-fitting f score Z = −3.6). Models that allow independent Denisovan introgression events into Negritos and Papuans are not rejected (worst-fitting f score of Z = −1.18 for both models), providing support that Philippine Negritos likely admixed with Denisovans independently from the Denisovan introgression event that occurred in Australopapuans.

Simulation of Demographic Scenarios
The coalescent simulator msprime29 was utilized to evaluate whether a null model (Figure 7F) of single Denisovan introgression event into the common ancestor of Ayta Negritos and Papuans or alternative models (Figures 7K, 7P, 7U, and 7Z) of multiple Denisovan admixture events are consistent with our empirical data: higher levels of Denisovan ancestry in Ayta Negritos than Papuans (Figures 7A and 7B), similar Denisovan tract lengths between Ayta Negritos and Papuans (Figure 7C), higher ratio of Denisovan over Australasian ancestry in Ayta Negritos relative to Papuans (Figure 7D), and Ayta Negritos is an outlier to the Papuan-Denisovan slope formed after correlating Denisovan and Australasian ancestry (Figure 7E). For each of the models shown in Figures 7F, 7K, 7P, 7U, and 7Z, 500 replicate simulations were performed, with 10 samples of 100M sequence taken from each of 8 populations. The parameters fixed across all models are shown in Table S1, and parameters that differed between models are shown in Table S2. All scripts used to generate data are available at github.com/jammc313/Philippine_Denisovan_Ancestry.

The null model (Figure 7F) includes an admixture event coming from an unsampled Denisovan population that diverged from the ‘Altai’ Denisovan 300 kya, into the shared ancestral population of Papuan and Ayta (4%, 48 kya). An alternative model (Alt1) follows an identical setup (Figure 7K), plus an additional admixture event (2%) coming from a distinct and unsampled Denisovan population, (also diverged from the Altai Denisovan 300 kya), into the (ancestors of) an Ayta population (40 kya), resulting in a total of 4% + 2% = 6% Denisovan admixture in the (ancestral) Ayta population. A second alternative model (Alt2) includes Denisovan admixture events private to both Papuans (4%, 45kya) and Ayta (6%, 35 kya) (Figure 7P). Again, the sources of Denisovan admixture are distinct Denisovan populations both diverged from the ‘Altai’ Denisovan 300 kya. A third alternative model (alt3) is similar to alt2, but with the Denisovan introgression event into Papuans occurring more recently at 25 kya (Figure 7U). The last alternative model (alt4) is similar to alt2, but with the independent Denisovan introgression events into Ayta and Papuans occurring around the same at 25 kya (Figure 7Z). The parameters of all models were kept as general as possible, with population sizes, admixture dates, and admixture proportions all chosen in accordance with the general consensus from previous studies24,27,28.

For all the simulations, we plotted the levels of Denisovan ancestry based on f4-ratio (Africa;Neanderthal,EastAsian,X)/(Africa;Neanderthal,EastAsian,Denisovan), mean tract lengths of Denisovan ancestry, and the ratio of Denisovan ancestry over Australasian ancestry. All pairwise comparisons were tested for significance using the Mann-Whitney U test with Holm-Bonferroni correction. The correlation between Denisovan and Australasian ancestry were also plotted following additional simulations with incremental 10% East Asian admixture into Ayta and Papuans.

The simulated models were leveraged against the empirical data using the following four parameters: level of Denisovan ancestry, mean Denisovan tract length, ratio of Denisovan ancestry over Australasian ancestry, and correlation between Denisovan ancestry and Australasian ancestry. Except for mean Denisovan tract length, the null model produced patterns of Denisovan ancestry that are not consistent with the observed data (Figures S7F–S7J). On the other hand, all the alternative models consistently replicate three pieces of empirical evidence: i) higher level of Denisovan ancestry in Ayta Negritos than Papuans, ii) higher ratio of Denisovan over Australasian ancestry in Ayta Negritos than Papuans, and iii) Ayta Negritos falls outside the slope formed by Papuan-related groups when plotting the correlation between Denisovan ancestry and Australasian ancestry (Figures S7K–S7Ad). Among the alternative models, only Alt4 model produced patterns of mean Denisovan tract length that is approximating similar levels in Ayta Negritos and Papuans, supporting a demographic scenario where the timing of the separate Denisovan introgression event into Ayta Negritos and Papuans happened around the same time, (Figures 7C and 7Ab).