We further assessed the hypothesis that the lack of a clearly bimodal East ISEA mismatch distribution could be related to our power to resolve the true mismatch distribution. To address this, we re-sampled Papuan blocks to mimic the distribution of block lengths in East ISEA based on rejection sampling (with replacement). We sampled 1000 sets of 2000 blocks, and generated mismatches against the Altai Denisovan. Plotting these mismatch distributions against the observed Papuan and East ISEA mismatch distributions ( Figure S2 B) demonstrates that the signal of a dual mismatch peak would be expected to be weak in our East ISEA sample based on the smaller number of long blocks available.

This phenomenon also explains why we are able to detect the multiple peaks in East Asia and Siberia, despite having considerably shorter blocks in these populations. For example, the 2000 longest blocks in these populations have an average length of just 43 Kb and 35 Kb respectively, yet we are still able to detect two Denisovan mismatch peaks due to the probable difference in population divergence driving these of 7500 generations (5000 versus 12500). In contrast, the longest 2000 blocks have an average length of 263 Kb in Papuans, and such long blocks appear necessary to capture two populations diverging from the Altai Denisovan in a narrower time frame, approximately 2750 generations apart (9750 versus 12500). Among East ISEA, the longest 2000 blocks have an average length of 152 Kb; Figure S2 C implies that this lack of resolution may be hiding the two peaks that we detect in Papua.

To explore these phenomena, we performed a simple simulation of expected mismatch for two populations diverging from a source population at two fixed times, tand t. The distribution of mismatches corresponds to a Poisson with meanwhere= 1.45 × 10and is the mutation rate per base pair per generation, l is the block length and c is a coalescent time. The values of c are repeatedly drawn from the distribution of coalescence times – either an exponential distribution with location parameter tand mean t+ 2N, or with location parameter tand mean t+ 2N. In this way, we are able to incorporate random coalescence and mutation. Using the values of t= 9750, t= 12500 and N= 100 from our simulation fitting (see below) for the two Denisovan-like populations contributing to Papua, and t= 5000, t= 12500 and N= 100 for the two Denisovan-like populations contributing to East Asia, we confirm that longer blocks are required to detect two signals of divergence ( Figure S2 C). This is particularly true when the divergence times are less widely separated, as is the case for Papua. The value used here for the divergence time of t= 5000 generations for the Asian-specific Denisovan introgression signal is an approximation based on the location of that mismatch peak and not based on explicit modeling.

We hypothesize that the challenge in resolving mismatch peaks using short blocks is related to two factors. First, mutations are discrete. In our analysis, it might be typical to observe 0 or 1 mismatching base pairs in a small 1 Kb block, corresponding, respectively, to extremely low and extremely high mismatch compared to the average mismatch of larger, ‘higher resolution’ blocks. Equivalently, it is not possible to observe a theoretically expected mismatch of, for example, 0.0005 in a 1 Kb block as 1/1000 = 0.001 and 0/1000 = 0.0; a better approximation of the expected value is possible with larger blocks. Second, the ratio of stochasticity in the mutation process along the branches of a coalescent tree versus the difference in coalescent times driving two mismatch peaks is also an important variable.

To assess the impact of using a different minimum block length cut-off, we determined the mismatch using thresholds from 0 to 260 Kb (Main Text Figure 3 A). To ensure mismatch accuracy is not being impacted by the alignability mask, we additionally required that blocks contained total called sequence over the minimum cut-off. We did not calculate the mismatch when less than 50 entirely non-overlapping blocks are involved in the analysis. The twin peaks obtain resolution at approximately 130 Kb (based on 3556 blocks). As expected for small block sizes, the overall resolution is poor below 50 Kb. We also explored the role of the minimum block length cut-off for our combined Siberia and East Asia sample (Main Text Figure 3 B) and East ISEA continental groups ( Figure S2 ). Again, the East Asian and Siberian specific peak only resolves when a minimum block length cut off ≥ 50 Kb is applied. The two peaks observed in Papua do not clearly resolve for East ISEA, due to the reduced sample size and small Papuan ancestry leading to fewer observed blocks in this region.

Focusing on the full set of 2000 largest Papuan blocks, we instead observed 226 entirely disjoint regions using the bedtools ‘merge’ function. Of these, 151 were > 0.5 Mb from any other region. Although the top ten most introgressed regions contributed 515 blocks, this leaves 216 disjoint regions with an average frequency of 4.8% in Papuans. The two highest frequency blocks, at 94/144 and 57/144 in the 72 (diploid) Papuan samples, lie outside both of the Papuan-specific mismatch peaks, the former having an intermediate mismatch and the latter an unusually large mismatch versus the Altai Denisovan.

Second, there is a clear dual mismatch peak in the Papuan population that is not apparent in other groups, which instead show a single divergent mismatch peak with some fluctuation in its exact placement. Given these fluctuations, we sought to confirm that the twin Papuan peaks were not a result of a very small number of common blocks being overrepresented due to high frequency in the sample, which could emphasize stochasticity in the coalescent and mutation processes. For example, in an extreme case, if a set of 2000 blocks consists of just 14 blocks near fixation then only a few independent coalescent histories might be represented in the mismatch distribution.

Two primary patterns emerge. First, we replicate the results ofin identifying a mainland Asian-specific peak with relatively low divergence to the Altai Denisovan. This peak is not limited to East Asian populations (in whom it was originally detected), but also extends into Siberia and the Americas. Interestingly, our American populations show some reduction in the extent of this peak, potentially suggesting that introgression was ongoing more recently than the divergence of American and East Asian populations.

Plotting these mismatches by continent yields the distribution pictured in Main Text Figure 3 C. The variation in the number of blocks > 50 Kb between continental groups reflects varying sample size, as well as the total amount of Denisovan introgression and the time since Denisovan introgression in different groups.

An informative way of assessing the relationship of our high confidence Denisovan blocks to the Altai Denisovan genome, which was used to extract them, is by mismatch analysis. As longer blocks are better able to resolve a mismatch distribution (explored in Figure S2 and below), we extracted the 2000 longest blocks from each continental population. For each block with > 50 Kb of total unmasked sequence data (counts in Table S3 F), we calculated the number of differences compared to the Denisovan reference (with Denisovan and Neanderthal heterozygous positions masked; see STAR Methods S3). We then calculated the effective block length by subtracting the portion of each block covered by the alignability mask from the total block length, and converted these values into a mismatch (difference/bp). As the number of differences per bp will be impacted by our masking of low quality and heterozygous archaic sites, quality control decisions (e.g., call rate > 99% filter) and the SNP calling protocol, it is not possible to directly translate mismatch distance into times based on a standard human genome mutation rate. We therefore chose to scale the observed mismatches by dividing the mismatches of each block by the average genome-wide mismatch rate between the Altai Denisovan and West Eurasian samples, m. We chose West Eurasians as our baseline mismatch rate because that population has the lowest amount of Denisovan introgression. Note that the average genome-wide mismatch rate between Altai Denisovans and West Eurasians, assuming 0% Denisovan ancestry in West Eurasians, reflects both the divergence time of humans and the Neanderthal/Denisovan clade, and the ancestral population size of the common ancestor of humans, Neanderthals and Denisovans. We calculated the scale factor by summing the total number of mismatches between each of the 75 samples × 2 = 150 West Eurasian haploid genomes and the Altai Denisovan, and dividing this by 150 times the total length of unmasked sequence in the dataset.

We next considered the possibility that our masking of heterozygous sites in the archaic genomes to simplify phasing was causing the double peak. Briefly, in a tree of a single introgressing haplotype X and two haplotypes of the Altai Denisovan, Dand D, there are two possible coalescent topologies – either the first coalescent event is between two Altai Denisovan haplotypes, (X,(D,D)), or between an Altai Denisovan haplotype and the introgressing haplotype, (D,(D,X)) or (D,(D,X)). As more homozygote mismatches are expected in the former case, our masking of sites that are heterozygous in one of the archaic genomes could generate a complex mismatch signal. To confirm that the masking of heterozygotes was not causing the multiple peaks, we re-phased the data retaining loci with archaic heterozygotes, and performed the CP and mismatch analysis on this data. This time, we trimmed the CP Denisovan set by simply removing blocks that were more than 99.9% overlapped by CP Neanderthal blocks. As before, the twin peaks are clearly visible ( Figure S3 A).

Before proceeding to analyze possible causes of the dual Denisovan mismatch signal in Papuans, we sought to confirm that it is observed in block sets retrieved using a variety of methods. We first confirmed that the signal is apparent both in the raw CP Denisovan and HMM Denisovan output ( Figures S3 A and S3B) to verify that our iterative trimming approach is not causing the signal. The HMM tends to detect longer archaic blocks than CP, which incorporates sequence with greater divergence from the Altai Denisovan, consistent with its lower specificity (higher West Eurasian Denisovan signal in Table S3 A). The signal remains ( Figure S3 B) when performing a similar trimming procedure to that described above on our HMM block set (see STAR Methods S9a), this time removing blocks that were i) less than 50% overlapped by CP Denisovan blocks, ii) less than 0.1% overlapped by Swindows, or iii) over 99.9% covered by Neanderthal HMM output. These criteria were chosen to obtain a high level of enrichment of Papuan Denisovan signal over West Eurasian Denisovan signal.

As a final consistency check, we sought to make use of a second ancient Neanderthal genome, the Vindija Neanderthal (). This sequence is known to be more closely related to the Neanderthal that introgressed into modern humans than the Altai Neanderthal, and thus may be better suited for extracting introgressed Neanderthal blocks (a 10%–20% increase is reported by). We used CP to extract introgressing blocks from Papuan and East ISEA individuals, this time with either the Vindija Neanderthal genome or both Neanderthal genomes as our Neanderthal group. The mismatch of these blocks against the Altai Neanderthal again shows a unimodal distribution ( Figure S3 F). We additionally repeated our trimming procedure of the CP Denisovan blocks to create high-confidence block sets, now with CP Vindija blocks or blocks identified by CP as affiliated with either Vindija or the Altai Neanderthal removed. Again, the bimodal mismatch distribution is observed ( Figure S3 C). If the more divergent peak were Neanderthal spillover signal, then we would expect it to be reduced by removing more Neanderthal introgressed sequence; such behavior is not apparent.

We confirmed the single mismatch peak using the HMM Neanderthal block set, removing blocks that were i) less than 50% overlapped by CP Neanderthal blocks, ii) less than 0.1% overlapped by Swindows, or iii) over 99.9% covered by Denisovan HMM output. This confirmed the unimodal mismatch distribution of Neanderthal introgressing blocks against the Altai Neanderthal ( Figure S3 E).

Interestingly, despite slight fluctuations, a single dominant Neanderthal peak is observed for all populations. This strongly suggests a) that any demographic cause of the dual mismatch peak relates to events occurring in the Denisovan population, rather than among Neanderthals or the Denisovan/Neanderthal common ancestor, and b) that the dual peak is not caused by a bioinformatic error or property of our methods (see further discussion below) that might give rise to a bimodal mismatch distribution against any archaic hominin.

To better understand the potential demographic implications of the dual mismatch peaks observed in introgressed Denisovan blocks among Papuans, we generated similar mismatch distributions based on our 2000 longest high-confidence Neanderthal-specific introgressed blocks for each continental group. These were generated using the same trimming protocol described above ( STAR Methods S9a), but starting with CP Neanderthal blocks and requiring overlap with both the Neanderthal HMM output and S, and only allowing limited overlap with CP Denisovan blocks (see Table S3 E for continental distributions of these blocks). As before, we only used blocks with > 50 Kb of unmasked sequence data. A large number of blocks remained for all continental groups (1978–1999 blocks), and the average block length ranged from 168 Kb (America, with 27 samples) to 287 Kb (Papua, with 72 samples). As with the Denisovan introgressed chunks, the average length of the 2000 longest blocks reflects both sample size and introgression history; the higher levels of Neanderthal introgression in the continental groups translates to longer chunks being successfully extracted. The mismatch for each continental group is shown in Figure S3 D. For ease of comparison, the curves are again scaled to the genome-wide average mismatch between West Eurasians and the Altai Denisovan.

A simple way to achieve a bimodality in a mismatch distribution is through a demographic model involving multiple sources of archaic introgression.recently inferred that two Denisovan-like ancestries introgressed into East Asians on the basis of a bimodal mismatch distribution. Before proceeding to analyze the mismatch distribution of Denisovan introgressed blocks in Papuans in this context, we sought to exclude other factors, including bioinformatic bias, and to confirm that the driving signal is consistent with introgression from two source populations on the Denisovan clade and not potential confounders.

In total, only 6 blocks were assigned to both D1 and D2 in 6 different genomic regions. We suspected that these cases reflect occasional misclassification due to the stochasticity of the mutation process; most blocks in a region are consistent with the D1 or D2 classification, but some spill over to the other peak. We expect a negative correlation between block frequencies in this case. Calculating the frequency as in described in STAR Methods S12, we compared the correlation between D1 and D2 frequencies in our small sample of 6 overlapping blocks with 1000 randomly selected sets of 6 non-overlapping block pairs. We observed a weak negative correlation trend in the set of overlapping blocks (gradient = −1.15, lower than 92% of resampled sets), consistent with the probability of occasional block misclassification.

We used the Gaussian mixture model to assign blocks to either the less-divergent Denisovan component (D1) or the more divergent Denisovan component (D2), classifying as ambiguous any block with less than 80% support for one model over the other (i.e., 0.2 < p(D1) < 0.8). We were able to classify 1538 of 1683 blocks in this manner, 718 to D1 and 508 to D2; we did not classify blocks outside the 0.1 < m D < 0.23 bounds. These block sets were used in most subsequent analyses. When studying the demographic implications of the bimodal mismatch distribution, we studied either the full mismatch distribution (when retrieving split dates of Denisovan populations), or the entire set of high-confidence Denisovan blocks, classifying blocks with m D < 0.1 to D1 and blocks m D > 0.23 to D2. This allows us to maximize the information in those analyses. Spillover from D0 or Neanderthal – into D1 or D2 respectively – will be minimal as those introgression sources are extremely limited in our high confidence Denisovan block set.

We additionally statistically confirmed that a bimodal mismatch distribution was supported for unique chunks > 180kb in Papuans. Here, we recorded the mismatch of sets of exactly overlapping chunks as their average mismatch, to limit the impact of high frequency chunks on the distribution. The bimodal distribution remains strongly supported (AIC = −2158.53, versus unimodal −2076.80). This confirms that the bimodal mismatch distribution is unlikely to be an artifact due to selection (i.e., a small number of introgressed regions with high frequency) or sampling effects.

To statistically confirm the bimodal distribution, we fitted a Gaussian and Gaussian mixture to the mismatch distribution of long (> 180 Kb; Main Text Figure 3 A) Denisovan introgressed blocks identified in Papuans for 0.1 < M< 0.23 using the Python package sklearn v. 0.19.1 (function sklearn.mixture.GaussianMixture;). We used blocks > 180 Kb because i) large blocks are needed to resolve a complex mismatch distribution ( Figure S2 C); ii) the mismatch distribution is relatively stable with a minimum block length of 180 Kb or more ( Figure 3 A); and iii) sufficient Papuan blocks (n = 1683) remain for downstream analysis using this minimum block length. The bimodal distribution was strongly supported (AIC = –5808.85, versus unimodal –5582.72). Note that the negative AIC values occur due to the high probability density in the range of mismatch values observed, with the probability density function of the Gaussian mixture model concentrated over a mismatch range less than 1; and that the difference between AIC scores rather than their values are relevant here. The Gaussians are distributed according to N(0.146, 0.018) and N(0.199, 0.015), and are weighted [58%, 42%] respectively, where N(μ,s) indicates a Normal distribution with mean μ and standard deviation s (Main Text Figure 4 A).

We additionally confirmed that the tendency for Denisovan blocks to be called on both chromosome copies of an individual – the Denisovan haplotype homozygosity within that individual – was similar between D1 and D2 blocks (10.9% and 10.2%, corrected χ 2 = 0.066, p = 0.80).

We second counted the number of blocks showing a pattern indicative of potential incorporation of human sequence – lower mismatch against the human sample and higher mismatch against the Altai Denisovan in the middle third. 3.2% of blocks were consistent with this pattern in D1, compared to 5.0% of blocks in D2 (non-significant Chi square statistic 1.87, p = 0.17). The patterns of mismatches over individual blocks suggest that incorporation of human sequence due to phasing errors is rare, and similarly rare in D1 and D2.

Ideally, we would directly profile switch point errors by comparing our Denisovan introgressed blocks to those retrieved using experimentally phased haplotypes. However, experimentally phased data are only available for two Australian samples in our dataset, and given the low Australian sample size, we only called D1 and D2 blocks in Papuans. We therefore searched for local variation in the mismatch within D1 and D2 blocks. For CP to call a block as Denisovan that has human sequence data incorporated into it, we expect that the human sequence data would need to be closely surrounded by real introgressing Denisovan sequence. If this were not the case, then CP is expected to simply bisect the Denisovan block. A Denisovan block containing human sequence is therefore expected to have a sharp increase in the mismatch against the Altai Denisovan and a sharp dip in the mismatch against a human genome toward the center of the block. We assessed this signal by dividing D1 and D2 blocks into thirds, and calculating the mismatch of each third against the Altai Denisovan and a West Eurasian sample. For this analysis, we used a minimally masked version of the dataset where the main phased data was combined with unphased data in which archaic heterozygous/low quality sites were not masked, and the call rate filter was not applied; heterozygote/homozygote mismatches in unphased regions were half-weighted to reflect the average 50% probability of applying to the block, and heterozygote-heterozygote mismatches were not counted. The West Eurasian was chosen to be sample LP6005441-DNA_G10, a Russian with approximately average high confidence Denisovan-specific introgression as compared to other West Eurasian samples.

Such a pronounced phasing error is unexpected, and would be easily observed. As a check, we asked whether the average recombination rate differs between D1 and D2 genomic regions. A higher recombination rate in one set of blocks would be expected to lead to faster recombination between human and introgressing haplotypes and could complicate phasing. Additionally, recombination rate is expected to generally impact local genetic variation as it determines the linkage of neutral regions to any nearby selected variation. We compressed the sets of blocks in the two components to their respective overlapping subsets using the bedtools function ‘merge’, leaving 86 D1 regions and 68 D2 regions. For each region we calculated the average recombination rate using the same combined HapMap phase II b37 genetic map as used for phasing (); the distribution of recombination rate between the two chunk sets was not significantly different (standardized 2-sample Anderson-Darling test statistic –0.60, asymptotic p = 0.71).

We then considered whether phasing errors might be higher in one block set than the other. For example, if a considerable amount of human sequence were erroneously incorporated into D2 blocks but not D1 blocks, this could drive up the D2 mismatch versus the Altai Denisovan. A simple calculation based on themodel of hominin evolution, with a split between a single introgressing Denisovan clade and the Altai Denisovan t= 11998 generations ago (N= 13249) and the Denisovan/Human split t= 20225 generations ago (N= 32671), suggests a considerable human input would be required for D2 blocks to be approximately 30% more divergent from the Altai Denisovan than D1 blocks. The approximate human component required is h where:such that h ≈0.25. This can be considered a conservative lower bound, in that many introgressed Denisovan haplotypes will survive into the larger Nregime and hence have greater mismatch versus the Altai Denisovan.

We then sought to determine whether the mismatch seen in the two block sets might reflect challenges related to alignability. The overall proportion of coverage by the alignability mask was extremely similar in the genomic regions covered by the two block sets (median: D1 = 6.55%, D2 = 6.52%), indicating that alignability challenges are unlikely to be leading to biases in sequence diversity.

These checks, along with the lack of a bimodal mismatch distribution when comparing introgressing blocks retrieved using the Altai Neanderthal genome against the Altai Neanderthal ( Figures S3 C–S3E) and the fact that our filtering process to retrieve high confidence Denisovan blocks requires that blocks overlap Soutput (which does not use a reference archaic genome), suggest that genetic properties of the Altai Denisovan genome do not substantially contribute to the mismatch difference between D1 and D2 blocks.

Second, we confirmed that the number of Denisovan low-quality alleles does not differ greatly between the D1 and D2 blocks sets (D1: 0.122%, D2: 0.124%; of the 23.1 Mb and 19.2 Mb of D1 and D2 sequence identified, respectively). While low-quality SNP calls were masked in the dataset used to identify the two block sets (see STAR Methods S3), a strong bias in the amount of low-quality calls in the genomic regions covered by one block set could imply relevant biases in region quality. We did not observe this.

While the high coverage of the Altai Denisovan genome argues for high confidence in SNP calls, a bimodal distribution of mismatch distances could be generated by certain quality biases. To rule this out, we performed two checks. First, a well-documented form of ancient DNA damage is cytosine deamination leading to C-to-T substitution (). If there were a strong bias in GC content in the genomic regions of the two block sets, with higher GC content in the D2 regions, then deamination of the ancient Altai Denisovan genome could increase the D2 mismatch between the modern (introgressing blocks) and ancient sequences. We therefore assessed GC content in the two Denisovan block sets using the UCSC Table Browser (GRCh37/hg19, Mapping and Sequencing, GC percent, gc5Base). The average GC content is very similar in both sets (39.75% for D1 and 38.75% for D2) and the distributions clearly overlap.

The set of Papuan samples that we study come from three different sources – 30 newly generated sequences, 25 samples fromand 17 samples from. We used exactly the same pipeline for mapping and base calling on all three datasets, and while we did not observe any obvious differences between data from the three sources during our quality control steps, we thought it possible that the bimodal distribution could conceivably be caused by this batch effect. However, blocks assigned to D1 and D2 were common over all three sample sources, and the number of blocks assigned to D1 and D2 in sequences from the three sources showed no significant deviation from random expectations (χ= 0.0997; p = 0.95).

Based on the lack of a clear genic/non-genic division between D1 and D2 blocks, their similar B values, and no pronounced frequency skew toward rarer D1 blocks, we do not interpret the mismatch difference between D1 and D2 as being driven by selection.

As D1 bocks have lower mismatch estimates compared to D2, they could have been under strong purifying selection since the time of introgression and might show a more pronounced skew toward rare blocks. We therefore performed a Mann-Whitney U test to assess whether the frequency distribution of blocks in D1 and D2 are different in our Papuan samples. There was a significant statistical difference in the frequency distributions (U = 3304, two-tailed p = 0.031), but summary statistics describing the distributions were similar (mean D1: 0.048, D2: 0.049; median D1: 0.028, D2: 0.021; standard deviation D1: 0.062, D2: 0.069). Importantly, the proportion of rare blocks with frequency < 5% was in fact lower in D1 (70%) than D2 (76%), which is consistent with neutrality.

As an additional test for differing levels of purifying selection, we asked whether D1 and D2 genomic regions differed in their B values. B values are measures of background selection over the genome based on observed diversity in an alignment of human, chimp, gorilla, orangutan and macaque (). After converting original B values to GRCh37/hg19 genome coordinates, we calculated the average B value over each D1 and D2 region. The distributions of average B values in D1 and D2 regions were not significantly different (standardized 2-sample Anderson-Darling test statistic −0.87, asymptotic p = 0.91). The total average and standard deviation for all D1 and D2 regions were 0.48 ± 0.24 and 0.50 ± 0.23, respectively, and hence statistically overlapping.

Purifying selection is expected to be focused on genic regions. We therefore assessed whether D1 blocks have more overlap with genes than D2 blocks. As when assessing recombination rate differences, we compressed the sets of blocks in the two components to their respective overlapping subsets using the bedtools function ‘merge’, leaving 86 D1 regions and 68 D2 regions. 79 and 60 regions overlapped genes (Ensemble 91 GRCh37) respectively (non-significant χ 2 , p = 0.63).

We second considered the possibility that the two peaks represented differential selection. Under this explanation, we might expect D1 blocks to be under strong purifying selection, reducing variation and mismatch, with D2 blocks evolving under neutrality.

g. The topology of D1/D2 blocks

Wall, 2017 Wall J.D. Inferring human demographic histories of non-African populations from patterns of allele sharing. We therefore attempted to assign D1 and D2 blocks to specific coalescent topologies by counting mutation sharing patterns and assessing consistency with the fifteen possible topologies implied by the four leaves of the tree (the block, Altai Denisovan, Altai Neanderthal, and human). For this analysis, we followed the phasing assessment (see STAR Methods S10e above) in using the Russian LP6005441-DNA_G10 to represent humans and a minimally masked dataset. We counted the number of mutations in the 16 possible sharing categories. For example, with notation 0 indicating the ancestral allele and 1 indicating the derived allele and in order [Human, Neanderthal, Denisovan, Block], a derived mutation that is unique to the block would contribute to mutation motif 0001, a derived mutation shared between the block and the Altai Denisovan would contribute to motif 0011, and a fixed derived mutation would contribute to mutation motif 1111. When counting, we downweighed the contribution of variants in unphased regions such that all possible mutation motifs represented were counted, in proportion to their probability given an equal chance of each variant being on each haplotype. The approach we take – studying the frequency of ancestral and derived variants in a set of samples of interest, and assessing the consistency of these with different phylogenetic tree topologies – is allied to that taken in recent work investigating the different demographic models of modern human history ().

We note that relating topologies to demographic histories is complex; for example, given a sufficiently large ancestral population size, each coalescent topology would have an approximately equal probability even if all blocks are Denisovan-introgressed. Nevertheless, studying the difference in topology proportions in D1 and D2 to clarify the history of these block sets can be useful.

I. Mismatch order. For a block to be consistent with a given topology in terms of Mismatch order, the order of the mismatches between its leaves must be consistent with that topology. For example, given topology (H,(N,(D,X))) and with m ij indicating the mismatch between leaves i and j, the mismatch order is required to be m D X < m D N ≤ m X N < m D H ≤ m X H ≤ m N H . For topology ((H,N),(D,X)), there are two mismatch order conditions, m D X < m D N ≤ m D H ≤ m X N ≤ m X H and m H N < m H D ≤ m H X ≤ m N D ≤ m N X . The specific cause of inconsistency in the m D N < m D X . For a block to be consistent with a given topology in terms of Mismatch order, the order of the mismatches between its leaves must be consistent with that topology. For example, given topology (H,(N,(D,X))) and with mindicating the mismatch between leaves i and j, the mismatch order is required to be. For topology ((H,N),(D,X)), there are two mismatch order conditions,and. The specific cause of inconsistency in the Table S4 inset is

II. Branch length order. The Mismatch order requirement above does not consider ancestral or derived status. This makes it robust to multiple hit mutations and incorrect ancestral state inference, but may reduce accuracy if such phenomena are rare. For a block to be consistent with a given topology in terms of Branch length order, the order of the average branch lengths, measured in shared derived alleles, between its leaves must be consistent with that topology. For example, given topology (H,(N,(D,X))), and keeping mutation motif notation order [H,N,D,X], the branch length order constraints are 〈 0001,0010 〉 < 0100 and 〈 0001 + 0011,0010 + 0011,0100 〉 < 1000 , where e.g., 〈 0001,0010 〉 denotes the average of mutation motif counts 0001 and 0010. For topology ((H,N),(D,X)), there are two mismatch order conditions, 〈 1000,0100 〉 < 〈 0010 + 1100,0001 + 1100 〉 and 〈 0001,0010 〉 < 〈 0100 + 1100,1000 + 1100 〉 . The specific cause of inconsistency in the 〈 0001 + 0011,0010 + 0011,0100 〉 > 〈 1000 〉 . The Mismatch order requirement above does not consider ancestral or derived status. This makes it robust to multiple hit mutations and incorrect ancestral state inference, but may reduce accuracy if such phenomena are rare. For a block to be consistent with a given topology in terms of Branch length order, the order of the average branch lengths, measured in shared derived alleles, between its leaves must be consistent with that topology. For example, given topology (H,(N,(D,X))), and keeping mutation motif notation order [H,N,D,X], the branch length order constraints areand, where e.g.,denotes the average of mutation motif counts 0001 and 0010. For topology ((H,N),(D,X)), there are two mismatch order conditions,and. The specific cause of inconsistency in the Table S4 inset is

III. Consistency threshold. For a block to be consistent with a given topology in terms of Consistency threshold, the ratio of the number of topology-inconsistent mutation motifs (ignoring multiple hits and ancestry errors) to the number of topology-consistent mutation motifs should be under some threshold, T c . For example, given topology (H,(N,(D,X))), and keeping mutation motif notation order [H,N,D,X], ( ( 0101 + 1001 + 0110 + 1010 + 1100 + 1011 + 1101 + 1110 ) / ( 1000 + 0100 + 0010 + 0001 + 0011 + 0111 ) ) < T c . For topology ((H,N),(D,X)), ( ( 0101 + 1001 + 0110 + 1010 + 1110 + 1101 + 1011 + 0111 ) / ( 1000 + 0100 + 0010 + 0001 + 1100 + 0011 ) ) < T c . The specific cause of inconsistency in the ( ( 0110 ) / ( 1000 + 0100 + 0010 + 0001 + 0011 + 0111 ) ) > T c . For a block to be consistent with a given topology in terms of Consistency threshold, the ratio of the number of topology-inconsistent mutation motifs (ignoring multiple hits and ancestry errors) to the number of topology-consistent mutation motifs should be under some threshold, T. For example, given topology (H,(N,(D,X))), and keeping mutation motif notation order [H,N,D,X],. For topology ((H,N),(D,X)),. The specific cause of inconsistency in the Table S4 inset is

IV. Subtree balance threshold. For a block to be consistent with a given topology in terms of the Subtree balance threshold, the absolute log-ratio of two tree branches that are predicted to be equal under that topology should be under some threshold, T s . For example, given topology (H,(N,(D,X))), and keeping mutation motif notation order [H,N,D,X], we have six conditions | l n ( 0001 / 0010 ) | < T s , | l n ( ( 0001 + 0011 ) / 0100 ) | < T s , | l n ( ( 0010 + 0011 ) / 0100 ) | < T s , | l n ( ( 0001 + 0011 + 0111 ) / 1000 ) | < T s , | l n ( ( 0010 + 0011 + 0111 ) / 1000 ) | < T s , | l n ( ( 0100 + 0111 ) / 1000 ) | < T s . For topology ((H,N),(D,X)) there are six conditions | l n ( 1000 / 0100 ) | < T s , | l n ( 0001 / 0010 ) | < T s , | l n ( ( 1000 + 1100 ) / ( 0001 + 0011 ) ) | < T s , | l n ( ( 0100 + 1100 ) / ( 0001 + 0011 ) ) | < T s , | l n ( ( 1000 + 1100 ) / ( 0010 + 0011 ) ) | < T s , | l n ( ( 0100 + 1100 ) / ( 0010 + 0011 ) ) | < T s . The specific cause of inconsistency in the | l n ( ( 0001 + 0011 + 0111 ) / 1000 ) | > T s . For a block to be consistent with a given topology in terms of the Subtree balance threshold, the absolute log-ratio of two tree branches that are predicted to be equal under that topology should be under some threshold, T. For example, given topology (H,(N,(D,X))), and keeping mutation motif notation order [H,N,D,X], we have six conditions. For topology ((H,N),(D,X)) there are six conditions. The specific cause of inconsistency in the Table S4 inset is We categorized the topologies of the blocks based on four heuristic criteria (see inset in Table S4 ):

c = 0.10 and T s = ln(1.5) ( Assessing support for different topologies using each of these constraints independently with T= 0.10 and T= ln(1.5) ( Table S4 ) reveals support for only five topologies – (H,(N,(D,X))), (H,(D,(N,X))), (H,(X,(D,N))), (N,(H,(D,X))) and ((H,N),(D,X)). These topologies reveal a strong signal of human as an outgroup to Neanderthal, Denisovan and the block sequence, and/or the Denisovan and block sequence forming a clade. Support is substantially greater for the (H,(N,(D,X))) topology, which follows the proposed population tree and suggests that incomplete lineage sorting is not driving the signal of the long blocks we have identified. Importantly, topology (H,(D,(N,X))) is similarly infrequent in both D1 and D2 blocks. This topology is consistent with the introgressing block being placed as most similar to the Altai Neanderthal and is expected to be more common in the D2 block set if the excess divergence were driven by spillover signal from Neanderthal introgression. We further note that we do not observe a D2 peak in all non-African populations, and especially that we do not see this peak in West Eurasians, despite each having considerable Neanderthal introgression – again arguing against this explanation of the D2 mismatch.

Topologies with either the Altai Denisovan or D1/D2 block acting as an outgroup to other genomes – which are consistent with H. erectus introgression into blocks or into the Denisovan – are extremely uncommon. While it is helpful to confirm the absence of such topologies, they are not expected given our block identification methodology as they do not allow for excess similarity between blocks and the Denisovan genome. In this context, the ((H,N),(D,X)) topology is interesting. This topology is consistent with, but certainly not unambiguous evidence for, shared H. erectus introgression into both an introgressing block and the Altai Denisovan genome. While the topology is common, importantly we do not see clear evidence that the frequency of this topology is substantially different in D1 and D2. This may imply either that the H. erectus introgression identified in the Altai Denisovan genome in other studies (see above) occurred earlier than the divergence of D2 from Altai Denisovan or later than the divergence of D1 such that it is Altai Denisovan-specific.

Other frequency differences in topologies are observed and can be informative regarding the demographic drivers of the D1 and D2 block sets. In particular, three topologies show prevalence differences – an increase in the prevalence of (H,(N,(D,X))) in D1 versus D2, and decreases in D1 versus D2 in both (H,(X,(D,N))) and (N,(H,(D,X))) ( Table S4 , column < D2 > / < D1 > ). The increased frequency of (H,(N,(D,X))) in D1 is consistent with D1 originating from a less divergent Denisovan clade, as the D1 blocks join the Altai Denisovan ancestor more recently and hence will be more likely to coalesce with the Altai Denisovan genome before other events. A similar argument applies to the (H,(X,(D,N))) topology – here, proportionally more of the D1 blocks have coalesced with the Altai Denisovan already by the time the Altai Denisovan and Altai Neanderthal share common ancestry, such that (H,(X,(D,N))) should be under-represented in D1 versus D2. Again, this is observed. The increased prevalence of (N,(H,(D,X))) in D2 is not clearly expected and may be associated with detection bias; the greater distance of D2 chunks from the Altai Denisovan render them harder to differentiate from Neanderthal introgression, such that there is a tendency to retrieve topologies in which the Altai Neanderthal is placed deeper in the tree.

The (H,(N,(D,X))) topology is consistent with a series of coalescent events that mirror the proposed population model – with an introgressing population showing greatest affinity to the Altai Denisovan, then the Altai Neanderthal, and then humans. As such, the hypothesis that D1 and D2 represent two divergent populations on the Denisovan clade makes clear predictions about the frequency of different mutation motifs, conditioned on a block showing this topology. Specifically, recalling the mutation motif notation [H,N,D,X], the counts should substantially reflect the branch lengths, as the block- and Altai Denisovan-specific branches are shorter for D1, the motifs 0001 and 0010 are expected to be rarer in D1 than in D2. Given an identical position of the Denisovan/Neanderthal common ancestor for D1 and D2, any shortening of these branches will lead to a corresponding lengthening of the common ancestral branch of D1 and the Altai Denisovan, such that the 0011 motif is expected to be more common in D1. We sought to confirm this pattern. Without conditioning on topologies, the Denisovan-specific motif 0010 was 23% more common in D2 than D1, and the introgressing block-specific motif 0001 was 33% more common. Conditioning on the (H,(N,(D,X))) topology using the mismatch order criteria, for example, places these motifs as 28% and 33% more common in D2 than D1, respectively. Providing further evidence of consistency, the Denisovan/block shared motif 0011 is appropriately shortened when conditioning on the (H,(N,(D,X))), such that it is 20% less common in D2. The frequency of this motif is approximately equal when we do not condition on the topology (3% more frequent in D2).

To summarize, the topology proportions in D1 and D2 blocks do not support the idea that D1/D2 mismatch differences are driven by H. erectus introgression into the Altai Denisovan, into D1 or D2 blocks independently, or both. The prevalence of the (H,(N,(D,X))) topology, and the topology differences between D1 and D2, are consistent with Denisovan-like introgression in Papuans originating from two populations on the Denisovan clade.

The analysis above is useful in teasing apart the proportions and qualities of coalescent topologies represented in D1 and D2, and can help to rule out specific causes of their mismatch distributions. However, especially given the complexity of block identification, we emphasize that the topologies show consistency with a demographic scenario of interest rather than discounting all other scenarios. The topology differences between D1 and D2 could be consistent with other scenarios, including introgression from a Neanderthal/Denisovan sister-clade or extremely complex bottlenecks on the Denisovan clade.

Having established the likely cause of D1 and D2 mismatches as complexity in archaic hominin introgression, we sought to further explore differences between the two block sets. We proceed by assessing evidence for different introgression dates based on block lengths, and by asking whether a model with introgression from two deeply divergent Denisovan-clade populations is supported by simulations.