Multiple episodes of interbreeding between Neanderthal and modern humans.

Neanderthals and anatomically modern humans overlapped geographically for a period of over 30,000 years following human migration out of Africa. During this period, Neanderthals and humans interbred, as evidenced by Neanderthal portions of the genome carried by non-African individuals today. A key observation is that the proportion of Neanderthal ancestry is ~12-20% higher in East Asian individuals relative to European individuals. Here, we explore various demographic models that could explain this observation. These include distinguishing between a single admixture event and multiple Neanderthal contributions to either population, and the hypothesis that reduced Neanderthal ancestry in modern Europeans resulted from more recent admixture with a ghost population that lacked a Neanderthal ancestry component (the 'dilution' hypothesis). To summarize the asymmetric pattern of Neanderthal allele frequencies, we compiled the joint fragment frequency spectrum of European and East Asian Neanderthal fragments and compared it with both analytical theory and data simulated under various models of admixture. Using maximum-likelihood and machine learning, we found that a simple model of a single admixture did not fit the empirical data, and instead favour a model of multiple episodes of gene flow into both European and East Asian populations. These findings indicate a longer-term, more complex interaction between humans and Neanderthals than was previously appreciated.

Introduction

When anatomically modern humans dispersed out of Africa, they encountered and hybridized with Neandertals [6]. The Neandertal component of the modern human genome is ubiquitous in non-African populations, and yet is quantitatively small, representing on average only ~2% of those genomes [6, 22]. This pattern of Neandertal ancestry in modern human genomes was initially interpreted as evidence of a single period of admixture, occurring shortly after the out-of-Africa bottleneck [6, 26]. However, subsequent research showed that Neandertal ancestry is higher by ~12–20% in modern East Asian individuals relative to modern European individuals [22, 18, 36].

Neandertals occupied a vast area of Asia and Europe at the time AMH dispersed outside of Africa (~75,000 BP [12]), and later Europe and Asia (~47–55,000 BP [20, 32]). Moreover, the breakdown of Neandertal segments in modern human genomes is indicative of a time-frame for admixture of 50,000–60,000 BP [26, 32] prior to the diversification of East Asian and European lineages. The genome of Ust’-Ishim, an ancient individual of equidistant relation to modern East Asians and Europeans, has similar levels of Neandertal ancestry as modern Eurasians, but found in longer haplotypes, consistent with an admixture episode occurring ~52,000–58,000 BP [5]. Given the extensive support for a single, shared admixture among Eurasians, there is extensive debate surrounding the observation of increased Neandertal ancestry in East Asians.

There are several hypotheses that may explain the discrepancy in Neandertal ancestry between Europeans and East Asians. It is possible that admixture occurred in a single episode, or ‘pulse’, of gene flow, but demographic and/or selective forces shifted the remaining Neandertal alleles into the frequencies we see in modern populations. Among these explanations are differential strength of purifying selection across Eurasia [27] and that modern Europeans lost part of their Neandertal ancestry through ‘dilution’ by a ghost population which was unadmixed [34, 16]. It is also possible that admixture occurred multiple times; the first pulse of Neandertal gene flow into the population ancestral to East Asians and Europeans was supplemented by additional pulses after both populations had diverged [34, 35].

Sankararaman et al. [27] proposed that differences in the level of Neandertal ancestry in East Asian individuals could be explained by their lower ancestral effective population size relative to Europeans, which would reduce the efficacy of purifying selection against deleterious Neandertal alleles [7]. However, Kim and Lohmueller [14] found that differences in the strength of purifying selection and population size are unlikely to explain the enrichment of Neandertal ancestry in East Asian individuals. This conclusion was further strengthened by Juric et al. [10].

Another hypothesis consistent with a single episode of gene flow is that Neandertal ancestry in modern Europeans was diluted by one of the populations that mixed to create modern Europeans [15, 16]. This population, dubbed ‘Basal Eurasian’, possibly migrated out of Africa separately from the population receiving the pulse of Neandertal gene flow, and thus had little to no Neandertal ancestry.

On the other hand, admixture may have occurred multiple times; the first pulse of Neandertal gene flow into the population ancestral to East Asians and Europeans was supplemented by additional episodes after both populations had diverged [34, 35]. The finding of an individual from Peștera cu Oase, Romania with a recent Neandertal ancestor provides direct evidence of additional episodes of interbreeding, although this individual is unlikely to have contributed to modern-day diversity [5]. However, Neandertal ancestry has remained relatively constant across tens of thousands of years of Eurasian history [19], suggesting that any additional admixture events must have been smaller scale than the initial episodes of interbreeding.

Here, we study the asymmetry in the pattern of Neandertal introgression in modern human genomes between individuals of East Asian and European ancestry. We summarize the asymmetric distribution of Neandertal ancestry tracts in the East Asian and European individuals in the 1000 Genomes Project panel in a joint fragment frequency spectrum (FFS) matrix. We first fit analytical models using maximum likelihood to explain the distribution of fragments in European and Asian individuals marginally. We then compare the joint FFS to the output of genomic data simulated under specific models of admixture between Neandertals and AMH to achieve a higher resolution picture of the interplay of different demographic forces. Our results support a complex model of admixture, with early admixture occurring before the diversification of European and East Asian lineages, and secondary episodes of gene flow into both populations independently.

Results

We constructed the joint fragment frequency spectrum by analyzing published datasets of Neandertal fragment calls in 1000 Genomes Project individuals [27, 33]. To avoid complications due to partially overlapping fragments and difficulties calling the edges of fragments, we computed fragment frequencies by sampling a single site every 100kb and asking how many haplotypes were introgressed with confidence above a certain cutoff (Figure 1). Our main results make use of fragments called with posterior probability of 0.45 in the Steinrücken dataset, although we verified robustness across a range of cutoffs and between datasets (Supplementary Material). The observed average proportion of Neandertal ancestry in European individuals was 0.0137 and 0.0164 in East Asian individuals, corresponding to an average enrichment of 19.6% in East Asian individuals (Supplementary Figure 5 shows how this quantity changes across cutoffs).

Figure 1:

Graphic representation of the process to make an FFS based on the posterior probability calls from Steinrücken et al. [33]. At a single position every 100 kb, introgression is assigned for individuals presenting a posterior probability above the global cut-off. The count of individuals between East Asian and European samples determines the cell in the FFS where that site is counted

We first developed analytic theory to understand what the FFS would look like in each population separately under different demographic models. To our surprise, we found that when looking only at the marginal distribution of introgressed fragment frequencies, the one pulse model and the dilution model are not statistically identifiable (Supplementary Material). On the other hand, the two pulse model is identifiable. Moreover, the analytic theory reveals that population size history only impacts the FFS within each population as a function of effective population size; intuitively, this arises because once fragments enter the population, their frequency dynamics only depend on the effective population size, rather than the specifics of the population size history. With this in mind, we developed a maximum likelihood procedure to fit a one pulse and a two pulse model to the European and East Asian marginal spectra (Methods), and found strong support for the two pulse model in both cases (Λ = 193.91 in East Asians, nominal p = 7 × 10−43; Λ = 212.64 in Europeans, nominal p = 6 × 10−47. Figure 2a and b). A subsequent goodness-of-fit test strongly rejected the fit of the one pulse model (p = 2 × 10−26 in East Asians, p = 0.0 in Europeans; χ2 goodness of fit test) but could not reject the fit of the two pulse model in either population (p = 1 in East Asians, p = 0.95 in Europe; χ2 goodness of fit test); see also Supplementary Figure 1, which shows the residuals of each fit. Thus, we concluded from analyzing each population in isolation that that the history of admixture was complex, and involved multiple matings with Neandertals.

Figure 2:

Individual and joint Fragment Frequency Spectra. a,b) Marginal FFS for the East Asian and European populations (first 20 bins, excluding the 0 bin). The lines represent the best fitted one pulse and two pulse model for each population. c) FFS of the Steinrücken et al. [33] introgression data. d)FFS of the Steinrücken et al. [33] introgression data projected down to 64×64 bins, as used to train the FCNN.

Nonetheless, looking at each population individually, we did not have power to estimate the relative contribution of dilution and multiple admixtures in shaping the patterns of Neandertal fragments seen between Europe and Asia. To gain a more global picture of the history of human-Neandertal interbreeding, we developed a supervised machine learning approach. A difficulty when simulating Neandertal admixture is the large number of free parameters associated with modeling multiple populations from which we have incomplete demographic information. Supervised machine learning applied to genomic datasets is becoming a popular solution for inference [for examples see: 28, 25, 30]. Of particular interest to this study, supervised machine learning has demonstrated the capacity for optimizing the predictive accuracy of an algorithm in datasets that cannot be adequately modeled with a reasonable number of parameters [29]. In practice, this results in the ability to describe natural processes even based on incomplete or imprecise models [29]. Supervised machine learning implementing hidden layers, or deep learning, is particularly effective in population genetic inference and learning informative features of data [30]. A definitive advantage of deep learning is how it makes full use of datasets to learn the mapping of data to parameters, allowing inference from sparse data sets [29]. Comparable likelihood-free inference methods, such as ABC, typically use a rejection algorithm, resulting in most simulations being thrown away. This necessitates a very large number of simulations for accurate inference [30, 1]. Deep learning methods also have the potential to generalize in non-local ways, allowing them to make predictions for data not covered by the training set [1, 29].

We simulated Neandertal admixture by specifying five demographic models with different numbers of admixture events (Figure 3), and produce FFS under a wide range of parameters. We used the simulated FFS to train a fully-connected neural network (FCNN). The trained network classified models successfully ~58% of the time, well above the 20% expected by chance, and was not overfit to the training data (Supplementary Figure 2). We then examined how the precision of the prediction changed when we required different levels of support for the chosen model (Figure 4a). Crucially, we see that when the classifier has high confidence in a prediction, it is very often correct, and that multiple pulse models are not often confused with the dilution model (Supplementary Figure 3).

Figure 3:

Representation of the five different demographic models simulated in msprime. Lines within each population phylogeny indicate different paths why which alleles could enter into modern populations, with red lines indicating Neandertal alleles and blue lines indicating non-Neandertal alleles. NEAN: Neandertal, ASN: East Asians, EUR: Europeans, BE: Basal Eurasians

Figure 4:

Results from the FCNN classifier. (a) Posterior probability that a chosen model is correct (precision), for all models under different levels of support for the chosen model. The x axis shows the probability cutoff that we used to classify models, and the y axis shows the precision. Each line corresponds to a different model. Simulated datasets in which no model surpassed the cutoff were deemed unclassified. (b) Posterior probability of the empirical introgression data matching each of the five demographic models, determined by the FCNN classifier.

Finally, we applied the trained FCNN to our empirical joint FFS (Figure 4b). Strikingly, we found that the FCNN supported our two most complicated demographic models, favoring a model with 3 pulses of admixture (Posterior probability ~ 0.55), and with a lower probability, a model with 3 pulses of admixture and dilution (Posterior probability ~ 0.44). These results are consistent across a range of cutoffs for calling introgressed fragments (Supplementary Figure 4), are robust to errors in fragment calling (Supplementary Figures 9 and 10), and dovetail with our maximum likelihood results showing that the best fit model must include multiple episodes of human-Neandertal interbreeding.

Discussion

Despite initial indications of a simple history of admixture between humans and Neandertals, more detailed analyses suggested that there might be additional, population specific episodes of admixture. By analyzing the joint fragment frequency spectrum of introgressed Neandertal haplotypes in modern Europeans and Asians, we found strong support for a model of multiple admixture events. Specifically, our results support a model in which the original pulse of introgression into the ancestral Eurasian population is supplemented with additional pulses to both European and East Asian populations after those populations diverge, resulting in elevated Neandertal ancestry in East Asians relative to Europeans. This is similar to a model recently proposed by Vernot et al. [35] for explaining differential levels of Neandertal ancestry across Europe, Asia, and Melanesia. Importantly, our results exclude a demographic model where the difference in Neandertal ancestry between Europeans and East Asians is driven primarily through dilution of Neandertal ancestry in Europe due to recent admixture with Basal Eurasians, a population lacking Neandertal ancestry. Nonetheless, we cannot exclude dilution as playing a role in the differences in Neandertal ancestry between Europe and East Asia; a model which includes multiple pulses of Neandertal introgression and dilution through Basal Eurasians was the second likeliest model in the five model comparison. Given the evidence that Basal Eurasians contributed to the modern European gene pool [16], we suspect that dilution does play a role in shaping the pattern of Neandertal ancestry across Eurasia. However, a large amount of dilution would be necessary if it were the only factor explaining the ~19.6% difference in Neandertal ancestry between Europe and East Asia, in contrast with recent work that inferred a smaller (~9.4%) contribution of Basal Eurasians to modern European individuals [11].

Several confounding factors could impact our inference. Although it is unlikely that differential purifying selection is responsible for the discrepancy between European and East Asian Neandertal ancestry [14, 34, 10], some Neandertal ancestry was likely deleterious [7, 10] and our models assume neutrality. However, the strength of selection against introgressed fragments is likely to be small compared to the demographic forces at work; moreover, there is relatively little evidence of strong differences in the strength of selection between different non-African populations [31, 4]. To explore the impact of selection, we obtained the FFS from simulations of deleterious Neandertal ancestry by Petr et al. [19] and asked if we classified their scenarios with selection as a two pulse model using maximum likelihood. We found that we rejected a one pulse model at the 5% level in only 1 out of 15 different simulations with selection, suggesting that we are not likely to misclassify selection against Neandertal ancestry as a two pulse model.

Of additional concern is power to detect fragments in each population. To address this, we implemented a model of fragment calling errors (Supplementary Material). Based on simulations done by Steinrücken et al. [33], we expect false positive rates of approximately 0.1% and false negative rates of approximate 1%; such rates do not cause substantial shifts in the FFS (Supplementary Figure 6). Moreover, after extensive simulations, we found that the neural network trained with errors is robust to false positive fragment calls at a rate of 0.2%, and produce consistent results when applied to the real data (Supplementary Material). Finally, in an attempt to see inside the “black-box” of the fully connected network, we examined how the weights propagate from each entry of the JFFS to the final assignments (Supplementary Material). In doing so, we found that moderate frequency haplotypes are most important to distinguish between models (Supplementary Figure 7), whereas errors in calling fragments are most likely to impact low frequency haplotypes. This, combined with the fact that our results are robust across two different datasets and a range of cutoffs for determining archaic ancestry, convince us that our results are robust to errors in fragment calling.

In addition, it is possible that some of the Neandertal ancestry in East Asia has been misclassified, and in fact originated from Denisovan introgression. The misclassified archaic fragments could then mimic the signal of additional pulses of Neandertal introgression. To address this concern, we inferred the position of Denisovan fragments based on data from Browning et al. [2] and masked 1.6% of positions across the genome (Methods). The masking resulted in 0.49% of sites called as Neandertal introgression in the Steinrücken et al. [33] data to be removed from the introgression data we used in all analyses. Though we do not believe we removed all misclassified Denisova ancestry, we think it is unlikely that a substantial enough proportion remains as to mimic the signal of additional Neandertal pulses. These problems are likely to be further resolved as our ability to make accurate introgression calls for the various ancient human populations improves in the future.

Our work provides additional evidence for the ubiquity of archaic admixture in recent human history, consistent with recent work showing that humans interbred with Denisovans multiple times [2]. Though we find that additional pulses of admixture in both East Asians and Europeans are necessary to explain the distribution of Neandertal ancestry in Eurasia, we are unable to settle why East Asians have elevated Neandertal ancestry. Interestingly, in contrast to Denisovans, there does not seem to be evidence of Neandertal population structure within introgressed fragments [2]. Combined with our results, this indicates that the Neandertal population or populations that admixed with Eurasians must have been relatively closely related. This is consistent with the established inference of a long-term small effective size across Neandertals [21, 22], which has held up to scrutiny despite some claims of a larger Neandertal effective size [23, 24, 17]. Thus, we believe that a likely explanation for our results is that gene flow between humans and Neandertals was intermittent and ongoing, but in a somewhat geographically restricted region. Differential levels of admixture between different Eurasian groups may primarily reflect how long those populations coexisted with Neandertals in that region.

Methods

Data

We obtained the joint fragment frequency spectrum by first downloading publicly available Neandertal introgression calls from two sources [27, 33]. The Sankararaman data consists of the location of introgressed fragments along the genome in each phased haplotype for the 1000 Genomes Project populations. The Steinrücken data consists of the probability of Neandertal origin in 500 bp windows of the genome across each phased haplotype for the Central European (CEU) and Han Chinese (CHB) and Southern Han Chinese (CHS) from the 1000 Genomes Project. We computed fragment frequencies by sampling a single site every 100kb from both sources of data. To compute the Joint Fragment Frequency Spectrum (FFS), we counted how many haplotypes were introgressed at each position. For the Steinrücken data, we called a site introgressed if it has a posterior probability of being introgressed above 45%. We then applied the 1000 Genomes accessibility mask (downloaded from ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/accessible_genome_masks/20141020.pilot_mask.whole_genome.bed). In the Supplement, we show that our results are robust to the cutoff and consistent between both datasets. We also masked Denisova fragments that were falsely called as Neandertal by downloading the S’ fragment calls from Browning et al. [2] and masking any fragment that matched Denisova > 35% and Neandertal < 25%, resutling in removal of 1.6% of the genome overall and 0.49% of Neandertal fragments. Given that the Denisovan and Neandertal populations diverged relatively early following the divergence of the Neandertal/Denisova lineage and the modern human lineage, it is unsurprising that a relatively small fraction of Denisova haplotypes were falsely assigned as Neandertal introgression. Finally, we masked the (0,0), (0,1), (1,0), and (1,1) position of the FFS matrix, in order to reduce the impact of false negative and false positive fragment calls.

Analytical Model

We model introgression of intensity f as injection of alleles at frequency f into the population at the time of introgression. In the one pulse model, this results in an exact expression for the expected fragment frequency spectrum under the Wright-Fisher diffusion model (Supplementary Material). Multiple pulse models could be solved analytically using a dynamic programming algorithm as in Kamm et al. [11], but we instead approximate the expected frequency spectrum by making the approximation that the probability of sampling k introgressed haplotypes in a sample of size n + 1 is the same as sampling k haplotypes in a sample of size n for large n (c.f. Jouganous et al. [9]). This results in closed-form expressions for the expected frequency spectrum under both the two pulse and the dilution model (Supplementary Material). With an expected frequency spectrum given parameters θ and model (k out of n haplotypes are introgressed|θ,), we compute the likelihood where x is the number of fragments found in k out of n individuals. We optimized the likelihood using scipy, and compared models using the likelihood ratio statistic where θ and correspond to the i-pulse model. Under the null, Λ should be χ2 distributed with 2 degrees of freedom (since there are 2 additional parameters in the two pulse model). Simulations in the Supplementary Material suggest that p-values under this model are well calibrated, despitre the impact of linkage; thus, we opted to report nominal p-values from our likelihood ratio test.

Simulations

We used msprime [13], to simulate Neandertal introgression into two modern populations with multiple potential admixture episodes and dilution from Basal Eurasians (Figure 3). For each replicate, we simulated the complete genomes for 170 European individuals, and 394 East Asian individuals, matching the sampling available from the 1000 Genomes Project panel. We used the human recombination map (downloaded from http://www.well.ox.ac.uk/~anjali/AAmap/; Hinch et al. [8]). In each simulation we mimicked our sampling scheme on the real data by sampling 1 site every 100kb and calling a Neandertal fragment by asking which individuals coalesced with the Neandertal sample more recently than the human-Neandertal population split time.

For each simulation, we drew demographic parameters, including effective population sizes and divergence times, from uniform distributions. For effective population sizes we used 500–5000 individuals for Neandertals, 5000–50000 for Eurasians, 5000–100000 individuals for the European and East Asian populations. For divergence times, we used 12000–26000 generations for Neandertals and humans, and 1300–2000 generations for the Eurasian split. The divergence between Basal Eurasians and Eurasians was fixed at 3000 generations. Lastly, we drew introgression times between 1500–3000 generations for gene flow into Eurasians, 800–2000 generations for gene flow into the European and East Asian populations. The time for the introgression event between Basal Eurasians and Europeans (dilution) was drawn from an uniform distribution of 200–2000 generations.

In order to ensure that our simulations focused on the correct parameter space, we constrained the resulting amount of Neandertal introgression in the modern European and East Asian genomes. The average Neandertal ancestry a was drawn from an uniform distribution between 0.01–0.03, and the difference in ancestry d between the East Asian and European populations was drawn from an uniform distribution between 0–0.01. We then determined the introgression intensity given a and d (Supplementary Material).

Machine Learning (FCNN)

Using the resulting joint FFS, we trained a simple fully-connected neural network (FCNN) to categorize a joint FFS into one of five demographic models. The network was implemented in Keras [3] using a TensorFlow back-end. The network used a simple architecture of three Dense layers (from 1024 nodes, to 512 nodes, to 64 nodes), each followed by a dropout layer (0.20).

Data availability

No novel datasets were generated or analysed during the current study.

Code availability

Our simulation pipeline is available at https://github.com/Villanea/Neandertal_admix/blob/master/n_admix_10.py, which generates the simulated genomes, identifies Neandertal introgression in 100kb windows, applies the 1000 Genomes accessibility mask, and outputs a FFS and accompanying parameters for each replicate.

The code used to implement the fully-connected neural network can be found at https://github.com/Villanea/Neandertal_admix/blob/master/Fully_connected_network.py.

Code for performing maximum likelihood estimation of FFS parameters can be found at https://github.com/Villanea/Neanderthal_admix/blob/master/sym_stat_theory.py.