Post by Admin on Jan 6, 2022 4:43:42 GMT
f. D1/D2 signal is not caused by selection
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.
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).
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 (McVicker et al., 2009). 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.
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.
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.
g. The topology of D1/D2 blocks
The network of interacting hominin populations in the Middle Pleistocene is becoming increasingly complex. One phenomenon that is included in some models (Lipson and Reich, 2017, Mallick et al., 2016, McColl et al., 2018, Prüfer et al., 2014, Skoglund et al., 2016), but not others (such as the main model in the Malaspinas et al., 2016 study) is a usually small Homo erectus component in the Altai Denisovan. Our approach does not allow us to identify genomic regions derived from H. erectus that may have introgressed into modern humans only, but not into the Altai Denisovan (but see STAR Methods S12). However, if the genetic contact between H. erectus and Denisovans occurred before the divergence of the Altai Denisovan and the Denisovan population that introgressed into humans, our D1 and/or D2 blocks could include regions with H. erectus ancestry, which were introduced into modern humans by a Denisovan population that was already pre-admixed with H. erectus. If so, there would be two categories of blocks identified as introgressed in the modern human – those derived from the Denisovan clade and those with an H. erectus origin – which could have different mismatch distributions and create a bimodal signal. While this phenomenon is likely to be rare if the proportion of H. erectus in Altai Denisovan is small (2–14%), some models incorporate a surprisingly large contribution (e.g., 66% in Figure 3 of Skoglund et al., 2016).
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 (Wall, 2017).
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.
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.
h. D1/D2 Denisovan lineages introgressed at different dates
Given the likely demographic origin of the D1 and D2 haplotypes, the question of whether D1 and D2 have different introgression dates is of particular interest. We therefore sought to estimate introgression dates based on haplotype lengths, which are expected to follow an exponential distribution with its decay parameter depending on the introgression date and the amount of introgression (Equation 1 in Gravel, 2012). The accuracy of this approach depends on the power of our methods to detect haplotypes of different lengths. Two factors are relevant. First, shorter haplotypes are expected to be harder to detect as signals of introgression may be mistaken for noise. Second, haplotypes are expected to be broken up by phasing errors; while the incidence of switch point errors is low in our data (see STAR Methods S4), and while our use of archaic genomes in phasing is expected to substantially reduce errors in introgressed regions, some errors are probable. The haplotypes that we assign to D1 and D2 are extremely long (> 180 Kb), such that the signal of introgression is clear, and the power of our methods is expected to be consistently high. Nevertheless, we caution that date estimates derived from this method may be better considered as relative rather than absolute dates.
We sought to use simulations to profile the potential of date estimation through the distribution of block lengths using a set of long introgressing blocks, and to confirm that fitting dates using longer blocks only does not lead to substantial biases. We first note that deviations from the exponential distribution are known to occur under certain combinations of introgression parameters (especially smaller population sizes and larger admixture proportions; see Figure 3 in Liang and Nielsen, 2014). We therefore assessed the correspondence between simulations and the exponential expectation for our parameter range of interest.
The results of 200 replicates of a Wright-Fisher forward-time simulation of 5 Mb chunks, with recombination rate 1.27 × 10−8 (average rate from the HapMap combined genetic map) and the chromosome discretized into 1 Kb segments for computational simplicity, using an introgression proportion of 0.04, haploid Ne of 8334 (Australian Ne in Malaspinas et al., 2016), and introgression times from 50 to 2950 generations in steps of 50 generations, show that the exponential fit is close but not exact (Figures S4A and S4B), even when all individuals in the population are sampled. We attempted to fit each simulation both by maximum likelihood (using the scipy.stats.expon.fit function), using either all blocks, those with minimum length 50 Kb, or those with minimum length 180 Kb. We were able to retrieve accurate introgression dates (Figure S4C), although there is a tendency to underestimate the introgression date for more ancient introgression times. In the regime of interest (introgression times < 1800 generations), the deviation is at most 10%–15%. These fittings confirm that it is possible to fit dates using longer block lengths only.
Additional challenges in inferring introgression dates arise from errors in blocks length estimation. We profiled these using the forward-time simulations and fittings as above, but now modified the introgressing block lengths to Berror = B + Laplace(μ, b) on sampling, where B is the error-free simulated block length and the Laplace distribution location parameter μ = 0, while scale parameter b = 20000. Choosing the Laplace distribution as our error model assumes equal probability to over and under-estimate block length with a constant probability of error per base pair. If Berror < 0, the block was discarded, capturing the difficulty in correctly identifying short blocks. The Laplace(0,20000) distribution has a cumulative density function at 2.5% and 97.5% of −59915 and 59915 bp, respectively, capturing very substantial errors in block sizes. Under these models, we are still able to achieve accurate fitting of dates (Figures S4B and S4D), although bigger biases arise when including smaller blocks in the fitting. Again, there is a slight tendency to underestimate introgression dates by 10%–15% when using longer blocks.
Introgression may well have occurred over many generations rather than as a single event, and we considered it probable that fitting using block lengths would emphasize the most recent introgression time. For example, if very weak introgression were (hypothetically) to occur up to the present and even a small number of very long introgressing blocks were sampled, this could lead to a very recent inferred introgression date, likely depending on the fitting procedure. We therefore repeated the forward-time simulations, this time simulating introgression over 520 generations (15080 years with a generation time of 29 years) at a rate 0.04/520 = 7.69 × 10−5. Note that the effective introgression rate is only marginally reduced compared to the single-event introgression model due to replacement of haplotypes that are already partially introgressed under this model, and the expected correction when there is relatively low total introgression (as in our case) is minimal. Measuring the simulated introgression date as the mid-point in the introgression process (e.g., weak introgression from 1090 to 1610 generations ago corresponds to a mid-point of 1350 generations), there is again a slight bias to infer more recent introgression dates (Figure S4E). For the introgression times of interest, this bias is again no more than 10%–15%.
Finally, we fitted the output of coalescent simulations generated in STAR Methods S10i, which builds on the Malaspinas et al., 2016 model with Denisovan introgression at 1353 generations ago. A slight bias toward inferring more recent introgression dates was observed, of approximately 5%. Adding errors to the sampled block lengths as above did not change the inference when using long blocks > 180 Kb.
The simulations above suggest that introgression date fitting based on block lengths is effective given our demographic parameters and our use of larger block lengths, even under a strong block length estimation error model or introgression occurring over many generations rather than as a single event. Nevertheless, in each case there is a slight bias toward recent dates, that is greatest when introgression is more ancient. The downward bias in date estimation is limited to 10%–15% for the times most likely corresponding to archaic introgression in humans, and is likely closer to 5%. Based on these simulations, we consider the dates we report to be probable lower bounds on introgression dates, with true dates up to 15% more ancient than our fitting suggest.
We proceeded to perform exponential fittings on the lengths of blocks in the D1 and D2 sets (Figures S4F and S4G) using the Python package statsmodels v.0.8.0 (Seabold and Perktold, 2010) and maximum likelihood fitting, assuming either a constant recombination rate or the combined HapMap genetic map (Frazer et al., 2007). We confirmed 95% CIs through a block-bootstrap procedure, whereby the genome was divided into 2 Mb chunks, consecutive chunks were combined if blocks spanned boundaries, and artificial samples were generated by sampling the same chunks over all individuals with probability proportional to chunk length (usually 2 Mb, but sometimes 4 Mb or more) until the total observed number of blocks corresponded to that expected from the data. When calculating the date of introgression, we assumed an introgression proportion of 2% for each of D1 and D2, such that half of a proposed total of 4% Denisovan introgression (Malaspinas et al., 2016) entered from each ancestry into Papuans, and a generation time of 29 years.
The results of the block length fittings are consistent with relatively recent introgression times. Under a constant recombination rate, we estimate dates of D1 introgression as 17.9 kya (95%CI 8.7–29.4) and D2 as 32.9 kya (95%CI 22.9–44.2). While there is no Papuan recombination map, we also sought to incorporate local recombination rates into the fitting by scaling block lengths by the average combined HapMap genetic map recombination rate over all blocks in the D1 or D2 sets, respectively. This maintains the approximately exponential distribution of block lengths. Under this model, we retrieved date estimates of D1 introgression at 29.8 (95%CI 14.4–50.4) and D2 at 45.7 (95%CI 31.9–60.7). The weight of probability supports younger dates within this range (Main Text Figures 4B and 5E). On average, D2 introgression is relatively older than D1 introgression, by 1.84 and 1.53 times for the two recombination models above, respectively.
To assess the robustness of this finding, we repeated the fitting after removing replicates of haplotypes observed in multiple individuals. In this way, we are seeking to observe recombination histories that are independent, and to limit the impact of haplotypes at higher frequency due to selection. Under a constant recombination rate and using unique haplotypes (Figures S4F and S4G), we estimate dates of D1 introgression as 20.2 kya (95%CI 10.4–33.5) and D2 as 29.5 kya (95%CI 23.2–36.1); under the HapMap scaled recombination rate, we estimate dates of D1 introgression as 33.7 kya (95%CI 16.8–57.7) and D2 as 44.3 kya (95%CI 34.4–55.4). The D2 average introgression date is now estimated as 1.46 or 1.31 times as ancient as D1, depending on the recombination model. The various block length fittings all suggest that D1 introgression was relatively more recent than D2 introgression.
Several other lines of evidence are consistent with D1 having a more recent introgression date. First, the variance in the frequencies of non-overlapping D1 blocks is less than that observed for D2 blocks (Fligner’s test for equal variance = 7.1, p = 0.008). After a pulse of introgression, we expect the variance in haplotype frequencies to increase as haplotypes drift away from their initial frequency. Second, there is structure in the geographic distribution of D1 and D2 introgression (Main Text Figures 5A and 5B). The amount of identified sequence in the D1 block set (including blocks with mD < 0.1) is significantly lower in Baining samples as compared to mainland Papuans (1.33 Mb per phased haploid genome versus 1.82 Mb on Papua; Welch’s t test T = –3.9, p = 4 × 10−4), while there is no evidence of a different amount of D2 (including blocks with mD > 0.23) sequence (1.28 Mb versus 1.37 Mb, T = –0.8, p = 0.41). To account for sampling differences between New Guinea (N = 52) and Baining (N = 16) and to estimate CIs around the observed mean, we resampled one million times (with replacement) 16 samples from each of two islands. While the average amounts of D2 per individual in both populations overlap (NG: 2.75, 95%CI 2.36–3.12; Baining 2.37, 95%CI 1.96–2.79), Baining has significantly fewer D1 chunks (NG: 3.64, 95%CI 3.10–4.19; Baining 2.60, 95%CI 2.19–2.99). To additionally assess whether the ratio of D1 between mainland Papuans and Baining is greater than the ratio of D2 between mainland Papuans and Baining, we resampled sets of mainland Papuan (N = 52) and Baining (N = 16) individuals with replacement one million times, recording whether the median pairwise D1 ratio was greater than the median pairwise D2 ratio. D1[Papuaresample]/D1[Bainingresample] was greater than D2[Papuaresample]/D2[Bainingresample] in 95.7% of resampling iterations.
Together, these geographical patterns raise the intriguing possibility that the D1 component is more typical of mainland Papua, and introgression may even have been ongoing in the time-frame of the split between Baining and Papuan populations. Similarly, there is a tendency for populations outside ISEA to show a Denisovan signal more consistent with D2 (Main Text Figure 3C), although the limited amount of Denisovan introgression means that the blocks contributing to the mismatches are shorter, and hence there is less resolution in the mismatch distributions. While it is possible that the reduced D1 signal in Baining samples is caused by weak Asian admixture (given the lack of evidence for the D1 signal in mainland Asia), this would additionally be expected to reduce the Baining D2 signal compared to mainland Papuans (not observed, see above) and generate an excess signal of Asian ancestry in the Baining as compared to mainland Papuans (not apparent in LOTER results, Main Text Figure 1A and Table S2; or in the PCA, Main Text Figure 1B, where the Baining cluster toward Australians rather than with the Asian-admixed Bougainville samples). We additionally note the detailed demographic analyses in Hudjashov et al., 2017, which strongly place the Baining as a recently separated Papuan population that does not harbor additional admixture signals from a wide range of other regional populations.
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.
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).
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 (McVicker et al., 2009). 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.
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.
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.
g. The topology of D1/D2 blocks
The network of interacting hominin populations in the Middle Pleistocene is becoming increasingly complex. One phenomenon that is included in some models (Lipson and Reich, 2017, Mallick et al., 2016, McColl et al., 2018, Prüfer et al., 2014, Skoglund et al., 2016), but not others (such as the main model in the Malaspinas et al., 2016 study) is a usually small Homo erectus component in the Altai Denisovan. Our approach does not allow us to identify genomic regions derived from H. erectus that may have introgressed into modern humans only, but not into the Altai Denisovan (but see STAR Methods S12). However, if the genetic contact between H. erectus and Denisovans occurred before the divergence of the Altai Denisovan and the Denisovan population that introgressed into humans, our D1 and/or D2 blocks could include regions with H. erectus ancestry, which were introduced into modern humans by a Denisovan population that was already pre-admixed with H. erectus. If so, there would be two categories of blocks identified as introgressed in the modern human – those derived from the Denisovan clade and those with an H. erectus origin – which could have different mismatch distributions and create a bimodal signal. While this phenomenon is likely to be rare if the proportion of H. erectus in Altai Denisovan is small (2–14%), some models incorporate a surprisingly large contribution (e.g., 66% in Figure 3 of Skoglund et al., 2016).
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 (Wall, 2017).
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.
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.
h. D1/D2 Denisovan lineages introgressed at different dates
Given the likely demographic origin of the D1 and D2 haplotypes, the question of whether D1 and D2 have different introgression dates is of particular interest. We therefore sought to estimate introgression dates based on haplotype lengths, which are expected to follow an exponential distribution with its decay parameter depending on the introgression date and the amount of introgression (Equation 1 in Gravel, 2012). The accuracy of this approach depends on the power of our methods to detect haplotypes of different lengths. Two factors are relevant. First, shorter haplotypes are expected to be harder to detect as signals of introgression may be mistaken for noise. Second, haplotypes are expected to be broken up by phasing errors; while the incidence of switch point errors is low in our data (see STAR Methods S4), and while our use of archaic genomes in phasing is expected to substantially reduce errors in introgressed regions, some errors are probable. The haplotypes that we assign to D1 and D2 are extremely long (> 180 Kb), such that the signal of introgression is clear, and the power of our methods is expected to be consistently high. Nevertheless, we caution that date estimates derived from this method may be better considered as relative rather than absolute dates.
We sought to use simulations to profile the potential of date estimation through the distribution of block lengths using a set of long introgressing blocks, and to confirm that fitting dates using longer blocks only does not lead to substantial biases. We first note that deviations from the exponential distribution are known to occur under certain combinations of introgression parameters (especially smaller population sizes and larger admixture proportions; see Figure 3 in Liang and Nielsen, 2014). We therefore assessed the correspondence between simulations and the exponential expectation for our parameter range of interest.
The results of 200 replicates of a Wright-Fisher forward-time simulation of 5 Mb chunks, with recombination rate 1.27 × 10−8 (average rate from the HapMap combined genetic map) and the chromosome discretized into 1 Kb segments for computational simplicity, using an introgression proportion of 0.04, haploid Ne of 8334 (Australian Ne in Malaspinas et al., 2016), and introgression times from 50 to 2950 generations in steps of 50 generations, show that the exponential fit is close but not exact (Figures S4A and S4B), even when all individuals in the population are sampled. We attempted to fit each simulation both by maximum likelihood (using the scipy.stats.expon.fit function), using either all blocks, those with minimum length 50 Kb, or those with minimum length 180 Kb. We were able to retrieve accurate introgression dates (Figure S4C), although there is a tendency to underestimate the introgression date for more ancient introgression times. In the regime of interest (introgression times < 1800 generations), the deviation is at most 10%–15%. These fittings confirm that it is possible to fit dates using longer block lengths only.
Additional challenges in inferring introgression dates arise from errors in blocks length estimation. We profiled these using the forward-time simulations and fittings as above, but now modified the introgressing block lengths to Berror = B + Laplace(μ, b) on sampling, where B is the error-free simulated block length and the Laplace distribution location parameter μ = 0, while scale parameter b = 20000. Choosing the Laplace distribution as our error model assumes equal probability to over and under-estimate block length with a constant probability of error per base pair. If Berror < 0, the block was discarded, capturing the difficulty in correctly identifying short blocks. The Laplace(0,20000) distribution has a cumulative density function at 2.5% and 97.5% of −59915 and 59915 bp, respectively, capturing very substantial errors in block sizes. Under these models, we are still able to achieve accurate fitting of dates (Figures S4B and S4D), although bigger biases arise when including smaller blocks in the fitting. Again, there is a slight tendency to underestimate introgression dates by 10%–15% when using longer blocks.
Introgression may well have occurred over many generations rather than as a single event, and we considered it probable that fitting using block lengths would emphasize the most recent introgression time. For example, if very weak introgression were (hypothetically) to occur up to the present and even a small number of very long introgressing blocks were sampled, this could lead to a very recent inferred introgression date, likely depending on the fitting procedure. We therefore repeated the forward-time simulations, this time simulating introgression over 520 generations (15080 years with a generation time of 29 years) at a rate 0.04/520 = 7.69 × 10−5. Note that the effective introgression rate is only marginally reduced compared to the single-event introgression model due to replacement of haplotypes that are already partially introgressed under this model, and the expected correction when there is relatively low total introgression (as in our case) is minimal. Measuring the simulated introgression date as the mid-point in the introgression process (e.g., weak introgression from 1090 to 1610 generations ago corresponds to a mid-point of 1350 generations), there is again a slight bias to infer more recent introgression dates (Figure S4E). For the introgression times of interest, this bias is again no more than 10%–15%.
Finally, we fitted the output of coalescent simulations generated in STAR Methods S10i, which builds on the Malaspinas et al., 2016 model with Denisovan introgression at 1353 generations ago. A slight bias toward inferring more recent introgression dates was observed, of approximately 5%. Adding errors to the sampled block lengths as above did not change the inference when using long blocks > 180 Kb.
The simulations above suggest that introgression date fitting based on block lengths is effective given our demographic parameters and our use of larger block lengths, even under a strong block length estimation error model or introgression occurring over many generations rather than as a single event. Nevertheless, in each case there is a slight bias toward recent dates, that is greatest when introgression is more ancient. The downward bias in date estimation is limited to 10%–15% for the times most likely corresponding to archaic introgression in humans, and is likely closer to 5%. Based on these simulations, we consider the dates we report to be probable lower bounds on introgression dates, with true dates up to 15% more ancient than our fitting suggest.
We proceeded to perform exponential fittings on the lengths of blocks in the D1 and D2 sets (Figures S4F and S4G) using the Python package statsmodels v.0.8.0 (Seabold and Perktold, 2010) and maximum likelihood fitting, assuming either a constant recombination rate or the combined HapMap genetic map (Frazer et al., 2007). We confirmed 95% CIs through a block-bootstrap procedure, whereby the genome was divided into 2 Mb chunks, consecutive chunks were combined if blocks spanned boundaries, and artificial samples were generated by sampling the same chunks over all individuals with probability proportional to chunk length (usually 2 Mb, but sometimes 4 Mb or more) until the total observed number of blocks corresponded to that expected from the data. When calculating the date of introgression, we assumed an introgression proportion of 2% for each of D1 and D2, such that half of a proposed total of 4% Denisovan introgression (Malaspinas et al., 2016) entered from each ancestry into Papuans, and a generation time of 29 years.
The results of the block length fittings are consistent with relatively recent introgression times. Under a constant recombination rate, we estimate dates of D1 introgression as 17.9 kya (95%CI 8.7–29.4) and D2 as 32.9 kya (95%CI 22.9–44.2). While there is no Papuan recombination map, we also sought to incorporate local recombination rates into the fitting by scaling block lengths by the average combined HapMap genetic map recombination rate over all blocks in the D1 or D2 sets, respectively. This maintains the approximately exponential distribution of block lengths. Under this model, we retrieved date estimates of D1 introgression at 29.8 (95%CI 14.4–50.4) and D2 at 45.7 (95%CI 31.9–60.7). The weight of probability supports younger dates within this range (Main Text Figures 4B and 5E). On average, D2 introgression is relatively older than D1 introgression, by 1.84 and 1.53 times for the two recombination models above, respectively.
To assess the robustness of this finding, we repeated the fitting after removing replicates of haplotypes observed in multiple individuals. In this way, we are seeking to observe recombination histories that are independent, and to limit the impact of haplotypes at higher frequency due to selection. Under a constant recombination rate and using unique haplotypes (Figures S4F and S4G), we estimate dates of D1 introgression as 20.2 kya (95%CI 10.4–33.5) and D2 as 29.5 kya (95%CI 23.2–36.1); under the HapMap scaled recombination rate, we estimate dates of D1 introgression as 33.7 kya (95%CI 16.8–57.7) and D2 as 44.3 kya (95%CI 34.4–55.4). The D2 average introgression date is now estimated as 1.46 or 1.31 times as ancient as D1, depending on the recombination model. The various block length fittings all suggest that D1 introgression was relatively more recent than D2 introgression.
Several other lines of evidence are consistent with D1 having a more recent introgression date. First, the variance in the frequencies of non-overlapping D1 blocks is less than that observed for D2 blocks (Fligner’s test for equal variance = 7.1, p = 0.008). After a pulse of introgression, we expect the variance in haplotype frequencies to increase as haplotypes drift away from their initial frequency. Second, there is structure in the geographic distribution of D1 and D2 introgression (Main Text Figures 5A and 5B). The amount of identified sequence in the D1 block set (including blocks with mD < 0.1) is significantly lower in Baining samples as compared to mainland Papuans (1.33 Mb per phased haploid genome versus 1.82 Mb on Papua; Welch’s t test T = –3.9, p = 4 × 10−4), while there is no evidence of a different amount of D2 (including blocks with mD > 0.23) sequence (1.28 Mb versus 1.37 Mb, T = –0.8, p = 0.41). To account for sampling differences between New Guinea (N = 52) and Baining (N = 16) and to estimate CIs around the observed mean, we resampled one million times (with replacement) 16 samples from each of two islands. While the average amounts of D2 per individual in both populations overlap (NG: 2.75, 95%CI 2.36–3.12; Baining 2.37, 95%CI 1.96–2.79), Baining has significantly fewer D1 chunks (NG: 3.64, 95%CI 3.10–4.19; Baining 2.60, 95%CI 2.19–2.99). To additionally assess whether the ratio of D1 between mainland Papuans and Baining is greater than the ratio of D2 between mainland Papuans and Baining, we resampled sets of mainland Papuan (N = 52) and Baining (N = 16) individuals with replacement one million times, recording whether the median pairwise D1 ratio was greater than the median pairwise D2 ratio. D1[Papuaresample]/D1[Bainingresample] was greater than D2[Papuaresample]/D2[Bainingresample] in 95.7% of resampling iterations.
Together, these geographical patterns raise the intriguing possibility that the D1 component is more typical of mainland Papua, and introgression may even have been ongoing in the time-frame of the split between Baining and Papuan populations. Similarly, there is a tendency for populations outside ISEA to show a Denisovan signal more consistent with D2 (Main Text Figure 3C), although the limited amount of Denisovan introgression means that the blocks contributing to the mismatches are shorter, and hence there is less resolution in the mismatch distributions. While it is possible that the reduced D1 signal in Baining samples is caused by weak Asian admixture (given the lack of evidence for the D1 signal in mainland Asia), this would additionally be expected to reduce the Baining D2 signal compared to mainland Papuans (not observed, see above) and generate an excess signal of Asian ancestry in the Baining as compared to mainland Papuans (not apparent in LOTER results, Main Text Figure 1A and Table S2; or in the PCA, Main Text Figure 1B, where the Baining cluster toward Australians rather than with the Asian-admixed Bougainville samples). We additionally note the detailed demographic analyses in Hudjashov et al., 2017, which strongly place the Baining as a recently separated Papuan population that does not harbor additional admixture signals from a wide range of other regional populations.