Mapping gene flow between ancient hominins through demography-aware inference of the ancestral recombination graph

The sequencing of Neanderthal and Denisovan genomes has yielded many new insights about interbreeding events between extinct hominins and the ancestors of modern humans. While much attention has been paid to the relatively recent gene flow from Neanderthals and Denisovans into modern humans, other instances of introgression leave more subtle genomic evidence and have received less attention. Here, we present a major extension of the ARGweaver algorithm, called ARGweaver-D, which can infer local genetic relationships under a user-defined demographic model that includes population splits and migration events. This Bayesian algorithm probabilistically samples ancestral recombination graphs (ARGs) that specify not only tree topologies and branch lengths along the genome, but also indicate migrant lineages. The sampled ARGs can therefore be parsed to produce probabilities of introgression along the genome. We show that this method is well powered to detect the archaic migration into modern humans, even with only a few samples. We then show that the method can also detect introgressed regions stemming from older migration events, or from unsampled populations. We apply it to human, Neanderthal, and Denisovan genomes, looking for signatures of older proposed migration events, including ancient humans into Neanderthal, and unknown archaic hominins into Denisovans. We identify 3% of the Neanderthal genome that is putatively introgressed from ancient humans, and estimate that the gene flow occurred between 200-300kya. We find no convincing evidence that negative selection acted against these regions. Finally, we predict that 1% of the Denisovan genome was introgressed from an unsequenced, but highly diverged, archaic hominin ancestor. About 15% of these “super-archaic” regions—comprising at least about 4Mb—were, in turn, introgressed into modern humans and continue to exist in the genomes of people alive today.

Introduction 1 structure of the input seqeunces. It works on unphased genomes and can accommodate 62 changing migration and recombination rates. ARGweaver-D is a Bayesian method, 63 using Markov chain Monte Carlo (MCMC) iterations to remove and "rethread" 64 branches into the local trees; as a result, the output of ARGweaver-D is a series of 65 ARGs that are sampled from the posterior distribution of ARGs conditional on the 66 input data and demographic model. From these, we can extract posterior probabilities 67 of introgression for any lineage at any genomic position. 68 Another recent method, dical-admix [20], is similar to ARGweaver-D in that it is 69 designed to accommodate generic demographic models, and takes the full haplotype 70 structure of the input sequences into account. However, there are some important 71 differences that make our approach more applicable to the complex history of ancient 72 hominins. dical-admix assumes that there are only a few admixed individuals, and 73 that other genomes are "trunk" lineages that help define the haplotype structure of 74 their respective populations. It therefore cannot infer admixture from an unsampled 75 population, nor is it designed to work when all individuals have some degree of admixed 76 ancestry. Additionally, ARGweaver-D can handle unphased genomes, which is 77 important since there are not enough Neanderthal or Denisovan samples to reliably 78 phase these archaic genomes. 79 After introducing ARGweaver-D, we present simulation studies showing it can 80 successfully detect Nea→Hum introgression, even when using a limited number of 81 genomes. We then use simulations to show that it can also detect older migration 82 events, including Hum→Nea, Sup→Den, and Sup→Afr, depending on the underlying 83 demographic parameters. We then apply this method to African humans and ancient 84 hominins, classifying 3% of the Neanderthal genome as introgressed from ancient 85 humans, and 1% of the Denisovan genome as introgressed from a super-archaic hominin. 86 In contrast to Nea→Hum introgression, we do not see any evidence of selection against 87 Hum→Nea introgression. 88

89
ARGweaver-D can estimate genealogies conditional on arbitrary 90 demographic model 91 ARGweaver-D is an extension of ARGweaver [13]that can infer ARGs conditional on a 92 user-defined population model. This model can consist of an arbitrary number of 93 present-day populations that share ancestry in the past, coalescing to a single panmictic 94 population by the most ancestral discrete time point. Population sizes can be specified 95 separately for each time interval in each population. Migration events between 96 populations can also be added; they are assumed to occur instantaneously, with the 97 time and probability defined by the user. 98 Recall that ARGweaver is a MCMC sampler, in which each iteration consists of 99 removing a branch from every local tree in the ARG ("unthreading"), followed by the 100 "threading" step, which resamples the coalescence points for the removed branches. This 101 threading step is the main engine behind ARGweaver, and is accomplished with a The gray horizontal dashed lines represent the discrete time points in the ARGweaver model, when coalescence and recombination events occur; migration and population divergence times are pre-specified by the user and rounded to the nearest "half time-point" (midway between the dashed lines). Migration is assumed to occur instantaneously at a rate p M specified by the uesr. This ARG currently has three haploid samples, indicated by the solid black lines. A fourth sample from the right-hand population is being threaded into the ARG, with the dotted black line representing one possible threading outcome. Each dot on the tree is a potential coalescence point for the new branch, representing a state in the threading HMM. The black dots are states from a population path with no migration, whereas the blue dots are from the migrant population path. Recombination events occur immediately before positions b 2 , b 3 , and b 4 , as indicated by the red X on the trees preceding those positions. The dotted red line shows the recoalescence of the broken branch, which defines the tree at the next site. The recombinations before b 2 and b 4 would be sampled after the threading algorithm, as they occur on the branch being threaded, whereas the recombination before b 3 is part of the ARG before the threading, and therefore not modified at this stage. Here, we only show a single tree in each interval between recombination events; the local tree is identical within each of these intervals. The lineage being threaded enters an introgressed state at position b 2 , and leaves it at b 4 . The transition probabilities of the HMM are calculated between each pair of adjacent states; the probability of migration p M is a factor in the transition observed at b 2 . It is not a factor at b 3 because the new branch is already in a migrant state. The transition probability at b 4 includes a factor of 1 − p M .
version. This is because coalescence is not possible unless two branches are in the same 111 population at the same time, and so the state space of potential coalescence points will 112 be a subset of the original state space. However, as migration events are added, 113 coalescence points in other populations become possible, and some coalescence points 114 may be reachable by multiple population paths (see Fig 1). Therefore, the complexity of 115 the algorithm can quickly increase. Whereas the original threading algorithm had an 116 asymptotic running time of O(Lnk 2 ) (where L is the number of sites, n the number of 117 samples, and k the number of time points), ARGweaver-D is O(Lnk 2 P 2 ), where P is 118 the maximum number of population paths available to any single lineage.

119
One way to improve the efficiency is to allow at most one migration event at any 120 genomic location. Note that this assumption still allows multiple lineages to be 121 introgressed at the same genomic position, if they are descended from a common 122 migrant ancestor. This assumption is reasonable when the number of samples is small 123 and the migration rate is low, and is set as a default in ARGweaver-D that we use 124 throughout this paper. It has two advantageous side-effects: it avoids strange parts of 125 the state space that could cause MCMC mixing problems (such as back-migrations, or 126 population label switching issues). It also means that if we are modelling introgression 127 from a "ghost" population such as a super-archaic hominin (from which we have no 128 samples), there will be at most one (migrant) lineage in the population at any location. 129 Therefore, the population size of ghost populations does not matter as coalescence will 130 not occur within them.

131
After running ARGweaver-D, it is straightforward to identify predicted introgressed 132 regions; they are encoded in each ARG as lineages that follow a migration band. By 133 examining the set of ARGs produced by the MCMC sampler, ARGweaver-D can 134 compute posterior probabilities of introgression across the genome; this can be done in 135 any way that the user would like: as overall probabilities of migration anywhere in the 136 tree, or probabilities of a specific sampled genome having an ancestral lineage that 137 passes through a particular migration band. For a diploid individual, we can look at 138 probabilities of being heterozygously or homozygously introgressed. Throughout this 139 paper we use the cutoff of p ≥ 0.5 to define predicted introgressed regions, and compute 140 total rates of called introgression for a diploid individual as the average amount called 141 across each haploid lineage.

142
More details of the ARGweaver-D algorithm are given in Fig 1 and  We performed a set of simulations to assess the performance of ARGweaver-D for 148 identifying Neanderthal introgression into modern humans. These simulations 149 realistically mimic human and archaic demography, as well as variation in mutation and 150 recombination rates (see Methods). We compared the performance with the CRF  Next, we predicted introgressed regions in two non-African human samples from the 157 Simons Genome Diversity Panel (SGDP), one European (Basque) and a Papuan. The 158 ARGweaver-D model used is illustrated in Fig 3; in this case only the "Recent migration" 159 bands were included. We compared to calls on the same individuals from the CRF. 160 Again, ARGweaver-D used two Africans, whereas the CRF used 43. And while the CRF 161 uses Africans as a control group, ARGweaver-D allows for introgression into any of the 162 human samples. The results are summarized in Fig 4. Overall, the two methods call a 163 large fraction of overlapping regions, but each method also produces a substantial 164 fraction not called by the other method (between 15-40%), and ARGweaver-D generally 165 calls more regions. While ARGweaver-D seems to have greater power in simulations, 166 another factor in the discrepancy may be that the ARGweaver-D segments were called 167 with the inclusion of both the Altai and Vindija Neanderthals, whereas the CRF calls 168 were produced with the Altai Neanderthal only. Both methods show a strong depletion 169 of introgression on the X chromosome, especially in the Basque individual.  African)/F4(Altai, chimp; Vindija, African) [11]. Both ARGweaver-D and CRF 183 predicted fewer elements (1.95% and 1.56%, respectively), compared to the F4 ratio 184 statistic (2.31%). Looking across the chromosomes, there is a higher correlation between 185 coverage predicted by ARGweaver-D and the expectation from the F4 ratio 186 (Spearman's ρ = 0.75), than between CRF and the F4 ratio (ρ = 0.51) ( Fig S1).

187
ARGweaver-D can detect older introgression events 188 We next did a series of simulations to assess ARGweaver-D's power to detect other We analyzed these data sets with ARGweaver-D using the model depicted in Fig 3,206 with only the "old migration" bands. As we do not have good prior estimates for the 207 migration time (t mig ) or super-archaic divergence time (t div ), we tried four values of 208 t mig (50kya, 150kya, 250kya, 350kya) and two values of t div (1Mya, 1.5Mya). We 209 generated data sets under all 8 combinations of t mig and t div , and then analyzed each 210 data set with ARGweaver-D under all 8 models, in order to assess the effects of model 211 misspecification on the inference.

212
The power of ARGweaver-D to detect introgression is summarized in Fig 5. The left 213 side of the plot represents simulations generated with t div = 1Mya, whereas the right 214 side used t div = 1.5Mya. Power to detect super-archaic introgression is clearly much 215 higher when the divergence is higher, but (as expected) does not affect power to detect 216 Hum→Nea introgression. Looking from top to bottom, the plots show the effect of 217 increasing the true time of migration. In the top plot with t mig = 50kya, only results for 218 Sup→Afr are shown because the archaic hominin fossil ages pre-date the migration time. 219 For all events, we see power decrease as the true migration time decreases.

220
For a given simulation set, the effect of the parameters used by ARGweaver-D are 221 generally more subtle. We note that power tends to be better when older migration 222 times are used in the model, even when the true migration time is recent; in particular, 223 the power when t mig = 150kya (red bars) is often much worse than when later times are 224 used, especially for the Hum→Nea event. Similarly, power is often better when t div is 225 set to 1Mya in the ARGweaver-D model, as opposed to 1.5Mya.

226
In summary, ARGweaver-D has reasonably good power to detect super-archaic 227 introgression when the divergence time is old, but power is more limited as the divergence decreases. The power to detect Sup→Afr is always lower than the power to 229 detect Sup→Den, as the African population size is much larger, making introgression 230 more difficult to distinguish from incomplete lineage sorting. For the Hum→Nea event, 231 we have around 50% power if the migration time is 150kya, and around 30% power 232 when it is 250kya.

233
False positive rates are less than 1% when a posterior probability threshold of 0.5 is 234 used (Fig S2). When analyzing the simulated data sets, we included two additional 235 migration bands in the ARGweaver-D model as controls: one from the super-archaic 236 population to Neanderthal (Sup→Nea), and another from humans into Denisova 237 (Hum→Den). The rates of calling these events were also less than 1% for all models.

238
Importantly, the rate of mis-classification is very low for all categories ( Fig S3); in 239 particular, the model can easily tell the difference between Hum→Nea and Sup→Den 240 events, despite both resulting in similar D statistics [7,8].

241
More details about the simulation results are available in the Supplementary Text. 242 One issue to note is that, although the simulated data sets were generated with a 243 human recombination map, the ARGweaver-D model used a simple constant 244 recombination rate. Performance is somewhat better when ARGweaver-D uses the true 245 recombination map, but in practice there are not enough Neanderthal or Denisovan 246 samples to generate a reliable recombination map, and there is no data to infer the 247 recombination map for the super-archaic population. The Supplementary Text also 248 shows results when we simulate more African individuals. We find that performance 249 does not improve as samples are added, so in the main text we focus on analysis with 250 two African samples (four haploid genomes).

252
We next applied the models from the previous section to real modern and archaic 253 human genomes. Our goals were to identify and characterize introgressed regions from 254 previously proposed migration events, as well as to see if we find evidence for other Simulation results. Each shaded box represents a set of simulations generated with a different value of t mig and t div . Each bar gives the basewise true positive rate for a particular migration event and ARGweaver-D model, using a posterior probability threshold of 0.5. The color of each bar represents the value of t mig used in the inference model (orange=50kya, red=150kya, purple=250kya, blue=350kya); shaded bars have t d iv=1.5Mya in the inference model, whereas solid bars use t div = 1.0Mya. Because the archaic hominin fossil ages are older than 50kya, results for t mig = 50kya (top) are only applicable for introgression into humans.
two Africans from the SGDP [21], two Neanderthals [2,8], the Denisovan [4], and a 257 chimpanzee outgroup. We again use the demography illustrated in Fig 3, with old 258 migration events only. We focus on the model with t mig =250kya and t div =1Mya, 259 because this model seemed to have high power in all our simulation scenarios, and 260 because our results suggest that it may be the most realistic (as discussed below). The 261 results using other models are consistent with those presented here, and are described in 262 the Supplementary Text.

263
An overview of the coverage of predicted introgressed regions is depicted in Fig 6; a 264 more detailed summary is given in Fig S4. The most immediate observation is that 265 Hum→Nea regions are called most frequently, at a rate of ∼ 3% in both the Altai and 266 Vindija Neanderthal. This number is almost certainly an underestimate, given that the 267 true positive rate for this model was measured between 30-55%. By contrast, only 268 ∼ 0.37% of regions are classified as Hum→Den. As no previous study has found 269 evidence for Hum→Den migration, this serves as a control, verifying that our false 270 positive rate estimated in simulations is likely fairly accurate, as we estimated a FP rate 271 of 0.41%.

272
Whereas there is a depletion on the X chromosome of archaic introgression into 273 humans, we see high coverage of Hum→Nea on the X for both Altai and Vindija. The 274 fact that it is somewhat higher on the X than the autosomes might be partially explained by increased power on the X; simulations suggest that power will be ∼ 20% 276 higher for this event when population sizes are multiplied by 0.75 ( Fig S5). Overall,

277
there is a lot of variation in detected introgression across the chromosomes, and several 278 autosomal chromosomes have higher predicted coverage than the X, including 1, 6, 21, 279 and 22 (Fig S4).

280
Although the Vindija sample is 70ky younger than the Altai sample [8], there is no 281 apparent depletion of human ancestry on Vindija compared to Altai on the autosomes, 282 suggesting that negative selection did not cause a significant loss of human introgressed 283 regions in the Neanderthal during that time. Several chromosomes do show drops in 284 coverage from Altai to Neanderthal, with the largest drop on the X chromosome (Fig   285   S4).

286
Other migrations are detected at lower levels. We identify 1% of the Denisovan 287 genome as introgressed from a super-archaic hominin, which is double our estimated 288 false positive rate for this event. The fact that we found much less than the ∼ 6% 289 estimated by previous methods [8] might suggest that the super-archaic divergence time 290 is closer to 1Mya, since we would expect to have more power with a higher divergence 291 time. Still, this analysis resulted in 27Mb of sequence that may represent a partial 292 genome sequence from a new archaic hominin. ARGweaver-D also predicted a small 293 fraction of the Neanderthal genomes as introgressed from a super-archaic hominin 294 (0.75% for Altai and 0.70% for Vindija). These amounts are only slightly above the 295 estimated false positive rates (0.65%), and the Sup→Nea event has not been previously 296 hypothesized.

297
One interesting aspect of Sup→Den and Sup→Nea regions is that, to the extent that 298 these predictions are accurate, there is the potential that this super-archaic sequence 299 was passed to modern humans through subsequent Den→Hum and Nea→Hum 300 migrations. We explored these regions further by intersecting them with introgression 301 predictions across the full SGDP data set. This analysis is detailed in the 302 Supplementary Text. It first confirms that most Sup→Den and Sup→Nea regions have 303 higher-than expected divergence to the Denisovans and Neanderthals (respectively) 304 across all humans, and not just the two African humans used by ARGweaver-D. 15% of 305 the Sup→Den regions overlap with sequence introgressed into Asian and Oceanian 306 individuals from Denisovans, and many of these regions also contain a high number of 307 variants consistent with super-archaic introgression. We also see that 35% of the 308 Sup→Nea regions are introgressed in at least one modern-day non-African human. We 309 also identified one region of hg19 (chr6:8450001-8563749) which appears to be 310 Neanderthal-introgressed and overlaps a Sup→Nea region. We compiled a list of 311 Sup→Den and Sup→Nea regions that overlap human introgressed regions, and the 312 genes that fall in these regions. These are given in Tables S1 and S2. 313 We examined lengths of all our sets of predicted regions, as they might be 314 informative about the time of migration. However, we find that there is strong 315 ascertainment bias towards finding longer regions, so that the length distributions are  Instead, we looked at the frequency spectrum of introgressed regions to gain insight 318 into the times of migration events. The older the migration, the more likely that an 319 introgressed region has drifted to high frequency and is shared across the sampled 320 individuals. For the Hum→Nea event, we observed 37% of our regions are inferred as 321 "doubly homozygous" (that is, introgressed across all four Neanderthal lineages). This is 322 very close to what we observe in regions predicted from our simulations with migration 323 at 250kya (38%), whereas simulations with migration at 150kya and 350kya had 324 doubly-homozygous rates of 10% and 55%, respectively. To further narrow down the 325 range of times, we did additional simulations with t mig =200, 225, 275, and 300kya, and 326 compared the frequency spectrum of introgressed regions after ascertainment with 327 ARGweaver-D. We find that the observed frequency spectrum is consistent with 200kya 328 < t mig < 300kya (Fig S6). The same approach suggests that t mig > 225kya for the for 329 the Sup→Den event ( Fig S7). against the Hum→Nea regions. We sought to look for other signals that might hint at 344 possible functional consequences of this event. 345 We first looked at deserts of introgression that were detected in [5]. They noted four 346 10Mb deserts in which the rate of both Nea→Hum and Den→Hum introgression is 347 < 1/1000. The coverage of Hum→Nea introgression within these deserts is shown in 348 Table 1; the fairly high coverages suggest that these deserts are unidirectional. For two 349 of the deserts, the Hum→Nea coverage is very high, especially in the Altai Neanderthal. 350 The third region is interesting as it overlaps the FOXP2 gene, which contains two 351 human-chimp substitutions that have been implicated in human speech [23,24], 352 although the Hum→Nea introgressed region is upstream of these substitutions (Fig S9).  The columns "cov Altai" and "cov Vindija" show the coverage of Hum→Nea within the given region on each Neanderthal; "num Altai" and "num Vindija" show the number of introgressed regions > 50kb. The genome-wide average coverage for Altai and Vindija is 0.034 and 0.033, respectively.
We next looked at all deserts of Nea→Hum ancestry, to see if this larger set of When chimpanzee alignments are available, the non-chimp allele is colored; otherwise the minor allele is colored. Lack of a color may mean that the haplotype has the chimpanzee or major allele, or that it has missing data. The phasing of the variants represents the final phase sampled by the ARGweaver-D algorithm. Here, the Denisovan is usually homozygous and shares variants with Africans and Neanderthals outside of the introgressed region; but within it, the Denisova 2 haplotype has many singleton variants, whereas Denisova 1 continues to share many variants with Neanderthals and Africans. randomly chosen genomic regions matched for size (Fig S10). population size, mutation and recombination rates) affect the power to detect regions. 365 While the overall levels of enrichment are therefore difficult to interpret, it is interesting 366 to note that the enrichment of functional regions (such as CDS, promoters, and UTRs) 367 tends to be higher in the Altai than the Vindija Neanderthal, which is the opposite  Neanderthals, we confirm that previously identified Nea→Hum deserts are not depleted 407 for introgression in the opposite direction. We do see a slight decrease in Hum→Nea 408 introgression on the X chromosome in the Vindija Neanderthal compared to the Altai, 409 which could be explained by weak negative selection removing some introgressed regions 410 in the ∼70ky that separate these fossils. An interesting question is whether this lack of 411 selection is because human introgression introduced healthy variation into the 412 Neanderthal genome, or because the Neanderthal population was too small for natural 413 selection to act against anything but the most harmful variants. However, without more 414 archaic samples, these questions will be challenging to answer.

415
ARGweaver-D also identified 1% of the Denisovan genome as introgressed from a 416 super-archaic hominin. Previous studies have estimated the total amount of Sup→Den 417 as roughly 6% [8], but this is the first study to be able to identify specific introgressed 418 regions. The fact that we only find a small fraction of the total amount suggests that 419 the introgressing population was not too highly diverged from other hominins; this low 420 power is much more consistent with a divergence time of 1Mya than 1.5Mya. Still, we 421 report 27Mb of putative super-archaic sequence from this previously-unsequenced 422 hominin, and we note that 15% of these regions have been passed on to modern 423 humans through Den→Hum introgression. It may be possible to obtain more of this 424 super-archaic sequence by applying ARGweaver-D to the set of 161 Oceanian genomes 425 recently sequenced [25], looking for super-archaic segments passed through the 426 Denisovans.

427
There have been several studies suggesting super-archaic introgression into various 428 African poulations [9,10,26]. However, ARGweaver-D only detected a small amount of 429 Sup→Afr introgression, which was somewhat lower than our estimated false positive 430 rate. One aspect to note here is that the power to identify introgression from an 431 unsequenced population is highly dependent on the population size of the recipient 432 population. The larger the population, the deeper the coalescences are within that 433 population, making it more difficult to discern which long branches might be explained 434 by super-archaic introgression. In the case of Africans, we used a population size of 435 23,700, which was our best estimate from previous runs of GPhoCS [7,12]. If we had 436 used a smaller population size, ARGweaver-D would have produced more Sup→Afr 437 predictions, but most of these would be false positives unless that smaller population 438 size is closer to the truth. Overall, we caution that the problem of detecting 439 super-archaic introgression into a large and structured population such as Africas is very 440 difficult, and that claims of such introgression need to be robust to the demographic 441 model used in analysis. It may not be possible to address the question of ancient 442 introgression into Africans without directly sequencing fossils from the introgressing 443 population. 444 We also explored some introgression events that do not have any support from 445 previous literature; namely the Hum→Den and Sup→Nea events. A priori, we expected 446 that levels predicted for these events would likely serve to confirm our false positive 447 rates in real data. However, it is also possible that there is some amount of these types 448 of gene flow, which has not been detected previously because it goes against the net 449 direction of gene flow. For the Hum→Den event, we predicted a slightly smaller fraction 450 (0.37%) than our predicted false positive rate from simulations (0.41%). For Sup→Nea, 451 we predicted 0.75% of the genome introgressed, which is slightly higher than our For all ARGweaver-D runs in this paper, the MCMC chain was run for 2000 iterations, 462 with the first 500 discarded as burnin, and ARGs sampled every 20 steps thereafter.

463
Except where otherwise noted, phase was randomized for all individuals and the phase 464 integration feature of ARGweaver-D was used (--unphased). We used site compression 465 throughout (--compress 10). We also used --start-mig 100, which disallows 466 migrations for the first 100 iterations of the sampler, enabling ARGweaver-D to 467 establish an ARG with a good general structure before exploring the migration space.

468
Recombination rate. Rather than use a recombination map calculated from modern 469 humans, which may not be accurate for ancient hominins, we used a constant  Mutation rate. For real data analysis, the mutation rate map was based on primate 477 divergence levels in 100kb sliding windows, using genome-wide alignments of human, 478 chimp, gorilla, orangutan, and gibbon sequences (see Supplementary Text for details), 479 and scaled to an average rate of 1.45e-8/generation/site. Simulated data sets were 480 generated by sampling rates from this map, and the same map was used for analysis. divergence times used were taken from [8], and population sizes from [7] (which were 484 based on estimates from GPHoCS [12]). When analyzing chrX, population sizes were 485 scaled by a factor of 0.75. When analyzing non-African humans, we only included the 486 "recent" migration bands from Neanderthals and Denisovans into humans, whereas when 487 looking for older introgression events, we excluded the "recent" bands as well as 488 non-African humans.

494
Migration events occurred at half-time points including 50, 150, 250, and 350kya. Note 495 that on this time scale, the European/African split is very recent, so that we did not  The coverage of introgressed regions for an individual is computed as one-half times 511 the coverage of heterozygous regions, plus the coverage of homozygous regions. In 512 theory, this fails to account for sites that switch between the heterozygous and 513 homozygous states without reaching the threshold for either, but in practice this occurs 514 at a negligible fraction of sites.

515
Analysis of hominin data 516 Data preparation 517 We ran a series of ARGweaver-D analyses on freely available hominin data, described in 518 Table 2. The panTro4 chimpanzee sequence was used as a haploid outgroup. The chimp 519 alignment to hg19 was extracted from the alignments of 99 vertebrates with human 520 available on the UCSC Genome Browser 521 (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/multiz100way). Any region 522 which did not have an alignment for chimp is masked in the chimp sequence. For each individual, we masked genotypes with quality scores less than 20 or sequencing 525 depths outside the range [20,80]. For each ancient individual, we also used the filters 526 recommended by [8] and provided here: 527 hg19.wgEncodeDukeMapabilityRegionsExcludable; ∼ 9% of the genome for which 532 SGDP genotype calls were not provided (85% of this set overlapped previously 533 mentioned filters). ARGweaver-D was run in 2.2Mb windows, but we excluded any 534 window for which any of the ancient filters, or the combined site filters, exceeds 50% of 535 bases. In total we analyzed 1,166 autosomal windows and 52 windows on the X 536 chromosome, covering 2.56Gb of hg19.

CRF calls
538 Introgression calls from [5] were downloaded from https://sriramlab.cass.idre. 539 ucla.edu/public/sankararaman.curbio.2016/summaries.tgz. As recommended by 540 the README contained therein, "set1" calls were used for Neanderthal ancestry in the 541 Basque individual, whereas "set2" calls were used for Denisovan ancestry in both the Neanderthal ancestry according to [11]). For this analysis we masked all sites that did 550 not have a filter level (FL) field of 9 in the SGDP individuals. For the Neanderthals we 551 used the same filters described previously.

552
Simulated data sets (deep introgression) 553 We performed a series of simulations to assess ARGweaver-D's ability to detect older 554 migration events. Each simulated data set consists of a 2Mb region with 5 unphased 555 diploid individuals and one haploid outgroup, mimicking the demographic histories and 556 sampling dates of the individuals from the real data analysis. All simulations were 557 produced with the software msprime [28]. 558 The population tree used in the simulations is identical to the one depicted in Fig 3, 559 and sampling dates correspond to the sample ages in Table 2. The human population 560 size history also corresponds to the one in Fig 3. For the archaic hominins, we simulated 561 a more detailed model of population size change, using piecewise-constant estimates 562 produced by PSMC [29] and published in [8]. For the Neanderthal population history, 563 we averaged the histories produced separately for the Altai and Vindija individuals, for 564 the time periods when they overlap. Similarly, we averaged the Denisova and 565 Neanderthal population size estimates during the time frame of their common ancestral 566 population (415-575 kya).

567
For each data set, a random 2Mb region of the autosomal genome was chosen as a 568 template region from which we chose recombination rates and mutation rates used to 569 generate the simulated data. We used the recombination map estimated from 570 African-American samples [30]. For the mutation map, we used the same map as in the 571 real data analysis (based on primate divergence levels). Missing data patterns were also 572 taken from the template region; we applied the same ancient genome masks and 573 mapability/blacklist masks to the simulated data. (We did not mimic the sequencing 574 depth or quality score masks, which affected a relatively small fraction of sites).

575
Overall we produced several sets of simulations, each consisting of 100 2Mb regions. 576 One set served as a control and contained no migration events. All other sets each had 577 three types of migration (Hum→Nea at a rate of 8%, Sup→Den at 4%, and Sup→Afr at 578 0.5%). The rates of each event were chosen so as to have enough events per data set to 579 be able to assess power, while still being less common than the non-migrant state. They 580 were also chosen (by trial and error) to produce roughly similar levels of predicted 581 introgression as observed in the real data. The simulated data sets varied in the 582 demographic parameters used (migration time and super-archaic divergence time). A 583 smaller set of additional simulations was produced with population sizes scaled by 0.75 584 to see how power might change on the X chromosome (see Fig S5).  Simulated data sets (Nea→Hum introgression) 591 We also did a smaller simulation study to assess performance on the Nea→Hum event 592 and compare performance to the CRF. Most of the settings were the same as above, 593 except that we sampled 86 haploid African lineages and 4 Europeans, along with the 594 two diploid Neanderthals and a haploid chimpanzee outgroup. Demographic parameters 595 were the same as above, except that a European population diverged from the African 596 population 100kya and had a initial size of 2100; at 42kya it experienced exponential 597 growth at a rate of 0.002, for a present-day population size of 37236. (These parameters 598 were roughly adapted from [31], but modified to reflect current smaller estimates of the 599 mutation rate in humans.) We then added 2% migration from Neanderthal into 600 Europeans at 50kya. In some supplementary analysis we also included 5% migration 601 from human to Neanderthal 250kya.

602
For this analysis only, we used true haplotype phases, in order to have a fair 603 comparison with CRF, which assumes phased samples.     Each row shows a true migration category, and each column shows the fraction of bases predicted in the category indicated at the foot of the column. The color of each bar represents the true parameters used in simulation, as indicated in the legend, with darker colors used for the older super-archaic divergence time. Multiple bars of the same color show results on the same data set, using an ARGweaver-D model with a different t mig . The value of t mig used by ARGweaver-D is not indicated in the plot, but increases from left-to-right: t mig = 50, 150, 250, 350kya, with 50kya only shown for Sup→Afr. All the models used t div = 1Mya; the plot with t div = 1.5Mya is nearly identical.    Fig S6; here we look at putative Sup→Den regions. Because there is only one Denisovan individual, there are only two categories: heterozygous or homozygous. Note that while we expect rates of heterozygosity to decrease with migration time, the confidence intervals here are wide, and there may be conflicting ascertainment effects that cause the apparent increase in heterozygous segments for the simulated data sets with oldest t mig values.  Fig S8. UCSC Genome Browser shot of a region with predicted homozygous Hum→Nea introgression in Vindija. This region on chromosome 10 has a high-probability introgressed region in both Vindija (but neither Altai) haplotypes. The top green bar indicates a predicted Hum→Nea region in Vindija, and below this is the posterior probability of introgression across the region in both Neanderthals. The variant track is similar to Fig 8. Here, we see almost identical haplotypes between Vindija and the Africans, whereas Altai shares haplotypes with the Denisovan.  We compared the distribution of various statistics across all non-overlapping 15Mb windows in the genome (black), to the distribution within deserts of Neanderthal introgression in humans of at least 10Mb (red). We excluded any window that crosses a telomere or centromere, or where ≥ 50% of the window does not pass our filters. In the bottom-right corner of each plot is shown the Kolmogorov-Smirnov statistic p-value, indicating that there is no significant difference between the black and red distributions. The statistics shown are indicated on the x-axis label. "Hum→Nea coverage" is average fraction of the window that contains any Hum→Nea region. "Mean Hum→Nea frequency" is the average number of introgressed haploid lineages of Hum→Nea across the window (where a frequency of zero indicates no introgression, and a frequency of 4 indicates homozygous introgression in Altai and Vindija). "Mean frequency of Hum→Nea regions" is the mean frequency, among regions with Hum→Nea calls. "Altai -Vindija coverage" is difference in mean coverage between the Altai and Vindija within each window.