nature ecology & evolution Article https://doi.org/10.1038/s41559-023-02145-2 Ghost admixture in eastern gorillas Received: 19 December 2022 Harvinder Pawar1, Aigerim Rymbekova 2,3, Sebastian Cuadros-Espinoza1, Xin Huang 2,3, Marc de Manuel1, Tom van der Valk4,5, Irene Lobon 1, Accepted: 30 June 2023 Marina Alvarez-Estape1, Marc Haber6, Olga Dolgova7, Sojung Han 1,2,3, Published online: xx xx xxxx Paula Esteller-Cucala 1, David Juan 1, Qasim Ayub 8,9, Ruben Bautista 8 , Joanna L. Kelley 10, Omar E. Cornejo 10, Oscar Lao1, Aida M. Andrés11, Check for updates Katerina Guschanski 12,13,14, Benard Ssebide15, Mike Cranfield16, Chris Tyler-Smith8, Yali Xue8, Javier Prado-Martinez1,8, Tomas Marques-Bonet 1,17,18,19,20 & Martin Kuhlwilm 1,2,3,20 Archaic admixture has had a substantial impact on human evolution with multiple events across different clades, including from extinct hominins such as Neanderthals and Denisovans into modern humans. In great apes, archaic admixture has been identified in chimpanzees and bonobos but the possibility of such events has not been explored in other species. Here, we address this question using high-coverage whole-genome sequences from all four extant gorilla subspecies, including six newly sequenced eastern gorillas from previously unsampled geographic regions. Using approximate Bayesian computation with neural networks to model the demographic history of gorillas, we find a signature of admixture from an archaic ‘ghost’ lineage into the common ancestor of eastern gorillas but not western gorillas. We infer that up to 3% of the genome of these individuals is introgressed from an archaic lineage that diverged more than 3 million years ago from the common ancestor of all extant gorillas. This introgression event took place before the split of mountain and eastern lowland gorillas, probably more than 40 thousand years ago and may have influenced perception of bitter taste in eastern gorillas. When comparing the introgression landscapes of gorillas, humans and bonobos, we find a consistent depletion of introgressed fragments on the X chromosome across these species. However, depletion in protein-coding content is not detectable in eastern gorillas, possibly as a consequence of stronger genetic drift in this species. Gorillas are a member of the great apes and form a sister clade to Homo The subspecies are distributed across western and eastern Africa (human) and Pan (chimpanzees and bonobos). Extant gorillas consist in a non-continuous manner (Fig. 1a). The current geographic ranges of of four recognized subspecies, which cluster into two species, a west- the different subspecies differ by size, continuity and ecology, impact- ern species of western lowland (Gorilla gorilla gorilla) and Cross River ing connectivity and population sizes5. Western lowland gorillas are (Gorilla gorilla diehli) gorillas and an eastern species of eastern low- endemic to a largely continuous range of considerable size, whereas the land (Gorilla beringei graueri) and mountain gorillas (Gorilla beringei other subspecies have much more fragmented distributions6. Likewise, beringei)1. All gorilla subspecies are either endangered or critically western lowland gorillas exhibit the highest genetic diversity of the endangered under IUCN criteria2–4. subspecies5,7,8, indicative of long-term high effective population sizes, A full list of affiliations appears at the end of the paper. e-mail: tomas.marques@upf.edu; martin.kuhlwilm@univie.ac.at Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 a b c GBB–Bwindi GBB–Virunga Mt Tshiaberimu GBB–Virunga 50 n=1 50 GGG Bwindi 0 PC2 (4.97%) PC4 (3.2%) n=5 0 GGD n=8 Mount Tshiaberimu Virunga n=7 –50 –50 Mount Tshiaberimu GBG GBG GBB–Bwindi GBB –100 –100 GGG GGD –100 –50 0 50 100 –100 –50 0 50 100 PC1 (16.74%) PC3 (3.92%) Fig. 1 | Gorilla samples used in this study. a, Present geographic distribution n = 12); GGG, Gorilla gorilla gorilla (Western lowland gorilla, n = 27); GGD, Gorilla of eastern gorillas, with that of the four gorilla subspecies shown in the inset. gorilla diehli (Cross River gorilla, n = 1). Shape files for the distribution of gorilla The newly sequenced samples are given in black, numbers of previously subspecies were obtained from IUCN. b, PCA with PCs 1 and 2 shown. c, PCA with sequenced eastern gorillas are given in colour. GBG, Gorilla beringei graueri PCs 3 and 4 shown. (Eastern lowland gorilla, n = 9); GBB, Gorilla beringei beringei (Mountain gorilla, while eastern gorilla effective population sizes are smaller9. Mountain eastern lowland gorillas (Extended Data Table 1). We performed a princi- gorillas are currently isolated in two discrete areas, the Virunga Volca- pal component analysis (PCA; Methods) to ascertain whether the newly noes Massif and the Bwindi Impenetrable National Park. The Bwindi sequenced individuals cluster with individuals from the same subspe- National Park is located at a lower elevation than the Virunga Volca- cies. The first PC separates western and eastern gorillas, as previously noes and as such has warmer temperatures10,11. Previous studies of the observed, and the second PC separates mountain gorillas from eastern demographic history of gorillas did not incorporate information from lowland gorillas (Fig. 1b). Since the new individual from the isolated all subspecies and were not conclusive, especially about the divergence Mount Tshiaberimu population clusters within the distribution of the time between the eastern and western clade9,12–15. This might be due other eastern lowland gorillas (Fig. 1b), this individual is, as expected, to gene flow from unsampled lineages, which is probably widespread considered a representative of this subspecies. The third PC reflects but is often insufficiently considered in evolutionary studies16,17. While population stratification within western lowland gorillas, whereas uncovering such hidden introgression events in gorillas is not pos- the fourth PC separates the eastern gorillas, with the two mountain sible from ancient DNA from fossil remains, as has been performed in gorilla populations from Virunga and Bwindi at the extremes (Fig. 1c), humans18, it is possible to address such questions using genomic data explaining 3.2% of the variance. This is well in agreement with previ- from present-day individuals19–22. ous studies7,8. To address this question, we use high-coverage whole-genome sequences of 28 western and 21 eastern gorillas. In addition to previ- ABC modelling favours a ghost lineage in eastern gorillas ously published genomes7,8, we sequenced the genomes of five moun- To infer a demographic model for the four extant gorilla subspecies, we tain gorillas from the Bwindi National Park and one eastern lowland used a neural-network based ABC modelling strategy using windowed gorilla from the isolated population of Mount Tshiaberimu. These summary statistics and extensive simulations (Methods; Extended Data new genomes contribute to a more complete representation of the Fig. 1), based on a previous implementation in the Pan clade22. A main genomic diversity of eastern gorillas. Using this expanded dataset, improvement is the implementation of a broad range of informative representing all four known gorilla subspecies, we explored the summary statistics (Supplementary Table 5), as is common practice demographic history of gorillas and specifically the hypothesis of for ABC studies in modern and ancient humans23. ghost introgression, defined as gene flow from an unsampled archaic We first established a demographic null model of the four popula- lineage. Given its substantial impact in their sister taxa of Pan and tions (Extended Data Fig. 2a and Supplementary Fig. 5), on the basis Homo as well as many other species18,21,22, such ghost introgression of previous studies3,7,8,24–27. Notably, although none of these studies events may explain some of the uncertainties in previous demo- incorporated whole-genome data from all subspecies, our inferred graphic models for gorillas. Using an approximate Bayesian com- parameters are largely coherent with previous work (Supplementary putation (ABC) approach, we find evidence for introgression from Tables 1 and 2). Nevertheless, unaccounted demographic events such an extinct lineage into the common ancestor of eastern gorillas and as ancient population structure or ghost admixture could affect param- characterize some of the functional consequences of this introgressed eter estimates16, particularly given evidence in other great apes22. Initial genetic material. exploratory analyses with f4-statistics and admixture graphs (Supple- mentary Section 2) did not show any asymmetries between the four Results gorilla terminal populations, which would arise if ghost admixture had Eastern gorillas form two population clusters occurred in any of the individual subspecies. However, this does not We newly sequenced six eastern gorillas to high coverage (on average, exclude the possibility of ghost admixture into the common ancestor 28.6×). After reprocessing the sequencing data from previous studies of eastern or western gorillas, which these methods cannot assess. To (Methods), we obtained a dataset of 49 individuals, with 27 western account for this and explicitly test if ghost admixture could improve lowland gorillas, one Cross River gorilla, 12 mountain gorillas and nine the inferred null demographic model (model A), we considered two Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 more complex demographic models, in which we added the possibility we implemented two independent approaches: the S* statistic19,20,29 of ‘ghost’ introgression into the common ancestor of eastern gorillas and the SkovHMM method or hmmix30. The S* statistic detects highly (model B) and western gorillas (model C). We assessed the robustness divergent windows relative to an outgroup, under a given demographic of our ghost models B and C using a wider parameter space (Supple- model, as introgressed sites19,20,29. Hence, the S* approach depends mentary Figs. 8–10; Methods), resulting in coherent posteriors with on the availability of a demographic model. By contrast, hmmix does those observed in models B and C (Supplementary Table 2), albeit not rely on a demographic model to identify putatively introgressed with wider confidence intervals (CIs), as expected given the increased regions but instead uses the density of private mutations in the ingroup model complexity (Supplementary Fig. 8). We performed a formal to partition the genome into ‘internal’ and ‘external’ fractions, walking comparison of these models (Methods), to determine which fits the in small windows of 1,000 base pairs (bp) along the genomes30. Hence, empirical data best. Model B, with archaic gene flow to the common although both S* and hmmix target the same signature of ghost intro- eastern ancestor had the highest posterior model probability of 0.9973, gression, the algorithms are distinct. compared to models A (0.0027) and C (0) and a substantially higher We simulated the expected null distribution of S* scores for east- Bayes factor (374 versus 0.0027 for model A and 0 for model C). In a ern gorillas with posterior parameter estimates from model A, that cross-validation analysis, the model with archaic introgression into is a model without ghost introgression. This yields insights into the eastern gorillas was clearly distinguishable from the model without presence of any outlier windows in our empirical data using the 99% (Supplementary Table 3). We conclude that a model with archaic intro- confidence interval (CI) for expected S* scores, given the mutation gression into the common eastern ancestor best explains the observed density (number of segregating sites) in each window (Supplementary summary statistics in the empirical data (for full posterior distributions Fig. 11; Methods). Indeed, at this threshold we identify an excess of S* see Extended Data Fig. 3). outlier windows, suggestive of introgression from an external source We infer that eastern gorillas experienced bottlenecks and gener- into the common ancestor of eastern gorillas: windows which fall ally had lower effective population sizes than western gorillas, while outside the null expectation constitute, on average, 1.64% of eastern mountain gorillas and eastern lowland gorillas experienced a par- lowland genomes and 2.36% of mountain gorilla genomes, respectively ticularly strong population decrease (Supplementary Tables 1 and 2), (Supplementary Table 9). as described previously8,24. We infer that the eastern subspecies split We assessed the performance of the S* statistic using coales- at 15,000 years ago (ka) (14–16 ka, 95% credible interval (CrI), Supple- cent simulations where we could trace the introgressed fragments mentary Tables 1 and 2). In agreement with previous studies13, we see (Methods). The precision and recall are high, with a 90.96% detec- a population expansion in western lowland gorillas ~40 ka. Our null tion rate of true introgressed fragments for eastern lowland gorillas demographic model infers a large ancestral population size for western (91.06% for mountain gorillas) at the 99% quantile (Fig. 2b, Extended gorillas (effective population size, Ne = 98,135), in comparison to that of Data Fig. 4 and Supplementary Table 7; Methods), comparable to the other gorilla populations, as well as a split between the western gorilla human–Neanderthal scenario31. Since the CrIs of the null demographic subspecies at ~454 ka (448–456 ka 95% CrI). Considering that not all model encompass larger effective population sizes, which would lead summary statistics could be calculated for Cross River gorillas (where to inflated rates of incomplete lineage sorting that might affect the only one sample was available) and gene flow between western gorilla expected distribution of S* scores, we also assessed how these param- subspecies was not modelled, we caution that the confidence in this eters influence our findings. Using the maximum values within the split time might be low. Finally, we infer that gorillas diverged into two 95% CrIs, we find that the recall of the S* statistic remains high, while species ~965 ka (729–1,104 ka 95% CrI), which is within the higher range the precision falls to 55.82% for eastern lowland gorillas (53.33% for of previous estimates9,12,28. mountain gorillas), reflecting an increase in the false discovery rate, as For simplicity, we modelled extant admixture as single migration expected. We conclude that the S* statistic performed well in detecting pulses over one generation, finding a small contribution of gene flow introgressed fragments under our null model, even when assuming from the common eastern ancestor to the western lowland gorillas of misspecification of the null model. 0.80% (0.06–2.14% 95% CrI), as well as from western lowland gorillas Analogous to previous work22, we also used hmmix to detect intro- to the common eastern ancestor of 0.27% (0.22–0.43% 95% CrI). We gressed windows30, which performs well for the given demographic infer a contribution of 2.47% of gene flow from an archaic source into model (Fig. 2b), with precision and recall well above 80% (Supplemen- the common ancestor of eastern gorillas, with a narrow 95% CrI of tary Table 8). Considering the strong support for ghost admixture 2.38–2.49% (Fig. 2c). We infer that this ghost population diverged from into eastern gorillas, we again used western lowland gorillas as the the extant gorilla lineages ~3.4 million years ago (Ma) (2.98–3.8 Ma, 95% outgroup and eastern gorillas as the ingroup. We find that 1.48–2.97% CrI). We estimate the timing of this ghost gene flow to have occurred of the individual eastern gorilla genomes are inferred as external at a 38,281 years ago, although the CrIs for this parameter are wide strict threshold for the mean probability of 0.95, with an estimated (22–108 ka, 95% CrI) (Fig. 2a,c). By contrast, the posterior distribu- introgression time of 37–41 ka. tions for the archaic introgression proportion and the gorilla–ghost While we observe sharing of the putative introgressed regions divergence time have narrow CrIs, indicating a strong support with across the eastern species, sharing is higher within each subspecies, clear peaks for these parameters (Fig. 2c). In contrast, our ABC analysis which again is more pronounced in the mountain gorillas than in the of model C does not confidently infer a contribution of a deeply diver- eastern lowland gorillas (Fig. 3a). This indicates that most of the puta- gent external lineage into the common ancestor of western gorillas. tive introgressed regions are segregating rather than fixed. Pairwise Instead, the best fit of this model suggests a 0.17% (0.09–0.4%, 95% CrI) nucleotide differences are elevated between eastern and western contribution from an external lineage at ~1.1 Ma into the common gorillas in putative introgressed regions in eastern gorillas, compared ancestor of all extant gorillas (Supplementary Table 2). This marginal to random regions (Fig. 3b). Likewise, there is an excess of nucleo- contribution is inferred to originate from an external lineage which tide differences between individuals of the eastern subspecies in the diverged from extant gorillas 1.9 Ma (1.5–3 Ma 95% CrI). putative introgressed regions, indicative of an archaic origin of these regions (Fig. 3b). The ghost introgression landscape in eastern gorillas The overlap of the autosomal hmmix fragments and the S* outliers Having established that a model of ghost introgression into the com- within each individual is, on average, 42% for eastern lowland gorillas mon eastern ancestor provided the best fit to the empirical data, we and 51% for mountain gorillas (Supplementary Table 11). For random aimed to identify the putative introgressed fragments in the genomes genomic regions passing filtering criteria, the observed overlap is, on of eastern gorillas. To explore this landscape of ghost introgression, average, 6% for eastern lowland and 8% for mountain gorillas, suggesting Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 a b 100 100 c 80 80 Precision (%) Recall (%) 60 60 35,381 Density 40 40 3.4 Ma (3.0–3.8 Ma) 20 20 0 0 14,363 ta ix ta ix G ta G G ta G ss m ss EL ss EL ss m rM rM 0 0.5 1.0 1.5 2.0 2.5 hm hm r r 965 ka Archaic introgression (%) 5,325 98,135 Density 454 ka 48,288 2.47% 40 ka (2.39–2.5%) 38 ka 18 218 418 618 818 1,018 (22–108 ka) 0.27% Time of archaic introgression (ka) 34 ka Ghost 0.80% 15 ka 3,009 21,982 Density 9 ka 115 7 ka 243 Present 2,157 20,894 64,748 14,558 Mountain Eastern Western lowland Cross river 1,140 1,640 2,140 2,640 3,140 3,640 lowland Time of gorilla–ghost split (ka) Gorillas Fig. 2 | ABC-based demographic model. a, Model of gorilla population history represent performance using the model presented in a, light bars represent the with archaic admixture from an unsampled ‘ghost’ lineage into the common ‘worst’ model with large Ne, in the case of hmmix to simulate the data to detect ancestor of eastern gorillas. The 95% CrI are shown for the archaic introgression fragments, in the case of S* to obtain the expected distribution of S* scores. proportion, timing of archaic introgression and archaic divergence (purple c, Posterior distributions for the archaic introgression proportion, time of timeframes), inferred under ABC modelling. Numbers on blocks represent archaic introgression and gorilla–ghost split time. The grey line indicates the effective population sizes. b, Precision and recall of hmmix (at the 95% posterior prior distribution. The red line represents the posterior inferred with neural probability cutoff) and S* (at the 99% quantile using sstar31) in simulated data networks. Neural networks reduce the dimensionality of the summary statistics using msprime. Precision (percentage of recovered introgressed fragments) and used and account for possible mismatch between the observed and simulated recall (percentage of true among inferred introgressed fragments) for hmmix summary statistics54. and S* (for ELG, eastern lowland gorilla and MG, mountain gorilla). Dark bars that both methods detect to a large degree the same regions (Fig. 3c of mutational conservation (SIFT, PolyPhen-2 and LINSIGHT scores) the and Supplementary Table 14). We thus consider the regions in the putative introgressed regions of both eastern gorilla subspecies follow intersect of the hmmix outliers and S* outliers as our high-confidence random expectation (Supplementary Fig. 19). We also investigated putative introgressed regions. The overlap between the two methods the distribution of gorilla-defined regulatory element annotations increases to 59% for eastern lowland and 68% for mountain gorillas, from another study33. Here, across categories and populations, we when using more lenient cutoffs for both methods, that is hmmix only observe an excess of strong enhancers (sE) in mountain gorilla fragments of at least 40 kilobases (kb) and 95% CI S* outliers (Supple- introgressed regions, compared to random regions (Supplementary mentary Table 12). Mountain gorillas (with the exception of Turimaso) Fig. 20). These are largely intragenic enhancers (Supplementary consistently exhibit higher proportions of overlapping base pairs of Fig. 21), which agrees with patterns of regulatory architecture observed the two methods than do the eastern lowland gorillas (Supplementary in primate sE33. Tables 11 and 12). Introgressed fragments can carry beneficial alleles and to explore signatures of adaptive introgression within eastern gorillas we applied The interaction of selection and introgression the method VolcanoFinder34. VolcanoFinder scans the genome for In contrast to archaic introgressed regions identified in humans and a signal of a distorted local site frequency spectrum consistent with bonobos, the putative introgressed regions in eastern gorillas are not a selective sweep surrounding an introgressed allele. Outliers of the significantly depleted in genic content compared to random genomic VolcanoFinder approach (95% composite likelihood ratio) within the regions (Fig. 3d). However, we find 127 megabases (Mb) of autosomal putative introgressed regions identified above were considered puta- segments longer than 5 Mb that are depleted for introgressed frag- tive targets of adaptive introgression. We identify seven candidate ments (Fig. 4). Further, we observe a signal of depletion in archaic regions for adaptive introgression (Extended Data Table 2), three of fragments on the X chromosome (Fig. 3f), on a scale comparable to which are shared between eastern lowland and mountain gorillas. The observations in modern humans32 and bonobos22. The putative intro- region with the highest likelihood ratio (LR) in VolcanoFinder (chr. 12: gressed regions of eastern lowland gorillas exhibit a slightly higher 11090005–11324172; maximum LR = 246.2) contains the bitter taste proportion of likely deleterious sites than do mountain gorillas, as esti- receptor TAS2R14, within which we find several protein-coding changes mated by the GERP score (Fig. 3e). However, under alternate measures (Supplementary Table 15). Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 b 500,000 EG–EG EG–WG WG–WG 400,000 a Random regions 50 300,000 Bwiruka Imfura 40 Kaboko Kahungye 30 Katungi Maisha 20 200,000 Nyamunwa Semehe 10 Tuck Turimaso 0 Umurimo 100,000 Zirikana 9732_Mkubwa 100,000 200,000 300,000 400,000 500,000 A929_Kaisi Introgressed regions A967_Victoria Dunia Itebero c d e f Mukokya 2.4 Overlap between methods (%) 50 High impact GERP content (%) 30 Ntabwoba Protein coding content (%) 3.5 Pinga Tumani 40 Autosome:X ratio 2.2 A967_Victoria Itebero Ntabwoba A929_Kaisi Mukokya Dunia Tumani 9732_Mkubwa Pinga Umurimo Kaboko Zirikana Tuck Imfura Maisha Turimaso Katungi Semehe Bwiruka Kahungye Nyamunwa 3.0 20 30 2.0 20 2.5 10 1.8 10 2.0 MG EL MG EL MG EL MG EL Fig. 3 | Characterization of introgressed fragments. a, Sharing of putative of high impact GERP content detected in introgressed regions (red lines) introgressed regions across eastern gorillas for autosomal regions detected and random regions (violin plots) for eastern gorillas. f, Autosome: X ratio of using the S* statistic and hmmix. b, Pairwise nucleotide differences in introgressed fragments inferred using hmmix for eastern gorillas (violin plots), introgressed regions (x axis) and in random regions (y axis) matched for length with reference lines for the equivalent values for bonobos (red line) and humans and proportion of positions with sufficient coverage (avoiding genomic regions (distribution as grey bar). In c–f: MG, mountain gorillas; EL, eastern lowlands. In without callable sites). Colours indicate the comparison: among eastern gorillas c–f, data are presented in violin plots with overlaid boxplots, which represent (EG–EG, green), among western gorillas (WG–WG, orange) and between eastern the median and interquartile range (25th and 75th percentiles). In f, individual and western gorillas (EG–WG, purple). c, Percentage of overlapping base datapoints are additionally plotted as black circles. For c–e, the data in violin pairs in introgressed regions (red lines) and random regions (violin plots) for plots consist of population-wise means for n = 100 iterations of random genomic eastern gorillas. For details of the definition of random regions see Methods. regions; for f, the data consist of hmmix fragments for n = 12 mountain gorillas d, Percentage of protein-coding content detected in introgressed regions and n = 9 eastern lowland gorillas. (red lines) and random regions (violin plots) for eastern gorillas. e, Percentage Discussion temporal sampling from historical specimens24,25 has the potential to Here, we present a demographic model inferred from representatives of be highly informative on variation lost over time. all four extant gorilla subspecies, leveraging the most comprehensive Previous estimates of demographic parameters varied greatly dataset of gorilla genomes available to date and an improved estimate under different models, methods and input data 9,13–15. The ABC for gorilla mutation rate from extended pedigree data35. The newly approach presented here leverages population-wise summary statis- sequenced whole genomes of mountain gorillas from Bwindi National tics. However, since high-coverage, population-level whole genomes Park are genetically close to those from Virunga but form a distinct are not currently available for Cross River gorillas, a subset of the cluster within their subspecies (Fig. 1b,c), confirming earlier results statistics could not be obtained for this subspecies (Methods; Supple- from microsatellite data36. Eastern lowland gorillas, as represented in mentary Table 5) and those calculated may be relatively less informative our dataset, seem to form a genetically less differentiated population, (for example, number of segregating sites). For all other populations, which includes the individual from Mount Tshiaberimu. Nonetheless, multiple individuals were included, yielding a better representation sample size remains a limitation, as high-quality invasive samples are of their diversity in the summary statistics. As such, we have lower highly restricted for endangered species, given ethical and logistical confidence in parameters involving Cross River gorillas, such as the constraints. A more fine-grained analysis of the evolutionary history relatively large divergence time inferred for the western lowland–Cross and population structure of gorillas necessitates denser sampling, River split. This divergence time represents 47% of the inferred eastern– which most likely will only be possible through advances in the use of western species split, compared to 26% estimated in a previous study non-invasive samples. For example, a reconstruction of recent pat- which also inferred a more recent species split time13. We note that this terns of connectivity has been demonstrated from a large panel of difference may be attributed to our inclusion of more western lowlands faecal samples from chimpanzees37. Furthermore, considering the gorillas, known to have high levels of population structure9,38,39. We rapid decline of great ape populations over the past centuries, more also do not include gene flow between western lowlands and Cross Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 13 14 90 Mb 0 Mb 0 Mb 90 Mb 15 12 90 Mb 0 Mb 90 Mb 0M b 0M b 16 b 90 11 M M 90 b 0 M b 17 b M 0 0 b M 18 M b 10 90 0M b 19 0M b 0M Mb 20 b 90 9 0Mb 21 0 Mb 0 Mb 22 90 Mb 0 Mb 8 0 Mb 90 Mb 1 90 Mb Mb 7 180 0M b b 0M 90 b M M 90 b 6 0 M 2 b b 18 M 0 0 M 18 b 90 b Mb 0M 5 0M b b 18 90M 0 Mb 90 Mb 0 Mb 180 Mb 3 4 Fig. 4 | Distribution of introgressed fragments. Outer circle: karyogram of in sliding windows of 2 Mb. Inner circle: long regions depleted of introgression the autosomes based on the human genome (hg19). Second circle from outside: content are shown in orange (length ≥5 Mb) and purple (length ≥8 Mb). Grey: introgression landscape in mountain gorillas (blue), as cumulative amount genomic regions with sufficient data (>20% of 40 kb windows passing threshold). of introgressed material in sliding windows of 2 million base pairs, Mb). Third White: genomic regions without sufficient data. circle from outside: introgression landscape in eastern lowland gorillas (green) River gorillas as a parameter in our modelling, which would reduce bonobos22 and other species41,42, as well as theoretical predictions and divergence estimates. simulations showing an impact of admixture from unsampled lineages The inferred deep divergence time between the two species is at that is probably common rather than exceptional16,17. Using extensive the upper end of previous estimates and conservative for the detection simulations, we find strong support for a model including archaic of putatively introgressed windows under the null model, since larger admixture into eastern gorillas, compared to a null model without such S* scores would be expected to result from an increased number of seg- ghost admixture or a model of such an event in western gorillas. The regating sites22. Indeed, even approximate demographic models with latter may be rather considered similar to a model of deep substructure large divergence times may allow a detection of external gene flow into within gorillas, given the shallower times and small amounts of external a target population31. We demonstrate that the S* statistic performed gene flow inferred. However, we note that further ghost introgression well in detecting introgression under the null model inferred herein events may exist beyond what we describe, for example with regards (model A), even if the true demography was deviating in terms of ances- to much smaller amounts of ghost admixture into gorillas or with shal- tral effective population sizes. Demographic modelling presented lower divergence times of the ghost lineages or in the context of larger here finds the best model for gorilla demography to include archaic effective population sizes in western gorillas. introgression from an unsampled ‘ghost’ lineage into the common Our inference of 2.47% of ghost introgression is associated with ancestor of eastern gorillas. This accords with a growing literature on high confidence as the posterior distribution is well differentiated from the prevalence of introgression from extinct lineages in humans21,40, that of the prior (Fig. 2c). This estimate agrees well with the estimates Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 of genome-wide introgression proportions per individual inferred individuals, we recover a putative 16.4% of the autosomal genome using the S* statistic and hmmix (Supplementary Table 9). We probably of an extinct lineage (Supplementary Table 13), adding to a growing underestimate the timing of archaic introgression, since shorter intro- literature revealing unsampled, now extinct lineages via analysis of gressed fragments are more likely to be missed and another potential variation present in present-day individuals. complication might be relatively high levels of homozygosity in eastern gorillas8, leading to increased haplotype lengths. Our definition of Methods putative introgressed regions as the overlap of outliers inferred with Samples and sequencing both the S* and hmmix methods (Fig. 3c) is conservative and on the Six eastern gorillas were sequenced as part of this study. Five Bwindi order expected for these methods, considering their relatively high mountain gorillas were sampled after death by the Mountain Gorilla false-positive rates31. Nonetheless, these methods are currently the Veterinary Project. One Mount Tshiaberimu individual was sampled only reliable tools available for detecting introgressed fragments in under anaesthesia. Convention on the Trade in Endangered Species of comparably small datasets of non-human species, without the avail- Wild Fauna and Flora (CITES) permits were obtained for all samples. ability of a source genome31. Sequencing was performed on the Illumina HiSeq X platform. Detailed A higher degree of sharing of putative introgressed fragments information on all samples is provided in the Supplementary Materials is observed among mountain gorillas than in eastern lowlands (Fig. (Extended Data Table 1). 3a). This is consistent with smaller effective population sizes of these populations, increasing the impact of drift on introgressed genetic Data processing variation18. High levels of genetic drift and reduced efficacy of natural We integrated the newly sequenced samples alongside previously selection probably also explain the absence of a detectable depletion published, high-coverage genomic data7,8. Raw sequencing reads of genic content in introgressed regions, in contrast to observations were mapped to the human hg19 reference genome, as described in introgressed material of humans and bonobos. Likewise, we do not previously52. Given that the hg19 reference does not belong to any of observe a coherent signature in mutational tolerance in gorilla intro- the gorilla subspecies, equal mapping bias will be exerted across all gressed material across different metrics, possibly due to genetic drift. gorillas in our dataset. This would not have been the case if the gorilla Despite this, we do find some ‘introgression deserts’, that is regions reference genome was used instead, as it might have introduced bias depleted of introgressed material in the population (Fig. 4), possibly in amounts of allele sharing, as observed previously for chimpanzees as a result of purifying selection18 shortly after the introgression took and bonobos52. The final dataset derives from 49 gorillas of known place. Furthermore, we observe a reduction of introgression on the X subspecies: 12 mountain (Gorilla beringei beringei), 9 eastern lowland chromosome, as also seen in humans and bonobos22,30,32. This is prob- (Gorilla beringei graueri), 1 Cross River (Gorilla gorilla diehli) and 27 ably a result of strong purifying selection against introgressed varia- western lowland (Gorilla gorilla gorilla) gorillas. tion, as seen in humans and other species18,30,43, possibly as a result of Processing of data to obtain genotypes followed procedures a combination with multiple factors44. Biased dispersal patterns45 and described in ref. 22. We used bcftools to retain genotypes with a cov- high reproductive skew of gorilla males46 might have led to a further erage larger than fivefold and lower than 101-fold, a mapping quality reduction of the male-haploid X chromosome in introgressed material. over 20, a proportion of MQ0 reads <10% and an allele balance >0.1 Even though the observed patterns are probably a combination of these at heterozygous positions; bedtools and jvarkit53 to filter the data by factors, we can currently not discern their respective contributions. known repeats (RepeatMasker) and mappability (35 k-mer). Following We note that our definition of adaptive introgressed targets is a previous study22, we used the rhesus macaque reference genome highly conservative, as the intersection of the outliers of three differ- (Mmul10) to infer ancestral allele states at each site and generate an ent methods S*, hmmix and VolcanoFinder as putative adaptive intro- ancestral binary genome, as implemented in the freezing-archer reposi- gressed targets. However, in being conservative we aim to minimize tory (https://github.com/bvernot/freezing-archer). Only positions the impact of potential false positives, which is a known caveat of the with genotype information in all individuals after filtering were used VolcanoFinder method34,47. However, at present this is the only method for calculating summary statistics for the demographic model and the available to localize signatures of adaptive introgression without a S* analysis. For hmmix, missing data were allowed, genotypes were fil- source genome. Interestingly, three candidate genes contain putative tered for known repeats and mappability and then an individual-based functional variants segregating in eastern gorillas and fixed ancestral filtering was applied for sequencing coverage (depth 6–100), mapping in western gorillas. One of these genes is TAS2R14, which encodes a quality (20) and retained only biallelic single nucleotide variants. taste receptor implicated in perception of bitter tastes48 and contains six missense variants. Eastern gorillas typically have more herbaceous Demographic modelling diets than the frugivorous western gorillas11, as such taste receptors Null demographic model. To infer a reliable null demographic model are plausible targets of adaptive introgression in eastern gorillas. for the four extant gorilla subspecies, we performed ABC modelling Bitter taste receptors have been suggested as targets of recent positive using the R package abc54 with neural networks, following a previously selection in western lowland gorillas as well, including a region encom- described strategy22. Previous demographic models did not include all passing TAS2R14 (ref. 13). It is possible that different mutations in the of the four extant gorilla subspecies8,13. We first attempted a merging of same region have been under selection in the different species. This these models (Supplementary Table 6) but in simulations this proved a could be interpreted in terms of the essential role of taste receptors poor fit to the empirical data in terms of the distributions of segregat- to avoid toxicity. The gene SEMA5A contains a missense variant and a ing sites, one of the main determinants of S* (Supplementary Fig. 3). splice region variant; this gene has been associated with neural develop- We used ms55 to simulate data and aimed to generate 35,700 coa- ment, with implications in autism spectrum disorder49. However, the lescent simulation replicates, of which 35,543 were successful, whereby functional impact of these variants in gorillas demands further work in per iteration we generated 2,500 windows of length 40 kb, randomly the future. Here, we do not find a contribution of adaptive introgression sampling from wide uniform priors informed by refs. 8,13,35 (Sup- to altitude adaptation, a phenomenon observed in humans and other plementary Table 2). We sampled local mutation rates from a normal species18,50. In mountain gorillas and eastern lowland gorillas at high distribution with mean of 1.235 × 10–8 (mutation rate per generation), altitude, this adaptation is probably driven by different mechanisms, recombination rate from a negative binomial distribution with mean of such as the oral microbiome51. 9.40 × 10−9 and gamma of 0.5 and assumed a generation time of 19 years In conclusion, our work contributes improved resolution to our (ref. 35). We scaled the mean mutation rate to 1.976 (1.235 × 10−8 × win- understanding of the evolutionary history of eastern gorillas. Across dow size of 40 kb × 4 × Ne of 1,000) with a scaled standard deviation of Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 0.460408 (1.976 × 0.233). We also scaled mean recombination rate to Detecting introgressed fragments 1.504 (9.40 × 10−9 × 4 × Ne of 1,000 × window size of 40 kb). Per window Following refs. 20,22,29 we calculated the S* statistic using a customized and per population, we calculated the following summary statistics: version of the package freezing-archer, accommodating non-human mean and standard deviations of heterozygosity, nucleotide diversity samples. We calculated the S* statistic genome-wide in 40 kb windows, (pi) and Tajima’s D, as well as the number of population-wise fixed and sliding every 30 kb, using the following test (i = ingroup) and reference segregating sites, the number of fixed sites per individual and pairwise (o = outgroup) populations: (1) GBG (G. beringei graueri-i and G. gorilla FST (Supplementary Table 5). These measures constitute the input sum- gorilla-o) and (2) GBB (G. beringei beringei-i and G. gorilla gorilla-o). For mary statistics for all ABC analyses performed in this section. Given the S* analysis, 15,181,832 variants were included. that only one diploid sample is available for G. gorilla diehli, we did Identifying outliers for the S* statistic requires a distribution of not use standard deviations of heterozygosity, nucleotide diversity scores for local mutation densities (represented by numbers of seg- and fixed sites per individual, as well as mean nucleotide diversity for regating sites in the dataset) under a demographic scenario without this population. introgression, as the null model. We used the weighted median poste- We calculated the equivalent summary statistics normalized by riors for each parameter value from the above ABC analysis to generate data coverage for the empirical data, which had been prefiltered by simulated data, specifying the number of segregating sites in a stepwise repeats, mappability and sufficiently informative windows (>50% of manner (from 15 to 800 in steps of 5). For each stepwise segregating sites with confident genotype calls in all individuals). We also filtered site (158 in total), we simulated 20,000 windows of length 40 kb, to by sites fixed across all gorillas relative to the human reference genome. which we applied the S* statistic for each of the scenarios (GBG and We accepted parameter values from the prior distribution if they gen- GBB). From this we obtained generalized additive models (GAMs) per erated summary statistics close to those of the empirical data. This scenario for three CIs (95%, 99% and 99.5%) using the R package mgcv, was assessed using a tolerance of 0.005, logit transformation of all following the procedures described in detail in refs. 22,29. From these parameters and 100 neural networks in the ABC analysis. GAMs, we predicted the expected S* distributions under the null model without archaic introgression. Applying the GAMs to the empirical Alternative demographic models. We performed parameter infer- data we inferred whether any windows lay outside the expectation per ence for two further demographic models, in which we allowed gene scenario and per confidence interval, assessing CIs of 95%, 99% and flow from a ‘ghost’ lineage into the common ancestor of (B) eastern 99.5%. As such, the threshold of significance is defined as the 95%, 99% gorillas (G. beringei beringei and G. beringei graueri) (Supplemen- or 99.5% CI from the standard deviation for expected S* scores, given tary Table 2) and (C) western gorillas (G. gorilla gorilla and G. gorilla the mutation density22,29. diehli) (Supplementary Table 2). For each alternate demographic To assess the performance of the S* statistic under our null model model, as above, we performed ABC analysis using 35,700 simulation and its robustness to model misspecifications, we performed valida- replicates, whereby per iteration we generated 2,500 windows of tion analyses following ref. 31, using msprime56,57 simulations with length 40 kb. We fixed parameters with narrow CrIs from model A, to explicit tracking of the introgressed fragments. Briefly, we simulated reduce the complexity of these models. To assess the impact of fixing expected distributions of S* scores for the null model (model A) and well-inferred parameters from the null model on subsequent ghost for a model where the effective population sizes before 40 ka were set parameter inference and explore the ghost parameter space more to the upper end of the 95% CrI (‘worst’ null model, in terms of highest fully we undertook a revised modelling approach (Supplementary expected amount of incomplete lineage sorting). We then simulated Section 3.4). In these revised ghost models, we performed parameter datasets of ten outgroup individuals (western lowland gorillas) and inference sampling all parameters from priors, for ghost gene flow a single ingroup individual (eastern lowland or mountain gorillas) into the common ancestor of (D) eastern gorillas and (E) western for model B and a model where the effective population sizes before gorillas (Supplementary Table 2). We observed a strong correlation 40 ka were set to the upper end of the 95% CI (‘worst’ model B). We between the estimated parameters of the original and the revised then obtained putatively introgressed fragments using the expected ghost models, albeit with wider posterior distributions for the revised scores from either model A or the ‘worst’ null model (Supplementary models due to increased complexity and larger parameter space Table 7). For each model, we performed 100 replicates and calcu- (Supplementary Section 3.4). lated the average precision and recall at different thresholds. For the To compare the three main demographic models—(A) null demog- S* approach, we used the quantiles of the S* statistic as thresholds, raphy, (B) ghost gene flow into the eastern common ancestor and (C) which range from 0 to 0.999. ghost gene flow into the western common ancestor—we simulated In an independent approach to the S* statistic, we applied hmmix30. 10,000 replicates of 250 windows of 40 kb length, fixing the parameters We obtained the input files for this method: weight files, local muta- as the weighted median posteriors for each model. To achieve an equal tion rates and individual observations files using scripts provided simulated timeframe (number of generations) in all models under with the repository for hmmix (https://github.com/LauritsSkov/ comparison, we added a non-interacting ghost population to the null Introgression-detection, as of 2 August 2018), as well as bcftools, bed- demography, with a divergence time between ghost and extant gorillas tools, jvarkit and custom R scripts. The macaque allele (RheMac10 equal to that inferred under Model B above. To determine if the models assembly) was used for polarization of alleles. We then applied the could be differentiated from each other we performed cross-validation method to the eastern gorillas using the following prior parameters: with the function cv4postpr (nval = 1,000, tol = 0.05, method = “neural- starting_probabilities = [0.98, 0.02], transitions = [[0.9995, 0.0005], net”). We calculated the posterior probabilities of each demographic [0.012, 0.988]], emissions = [0.05, 0.5]. We confirmed that using differ- model using the function postpr (tol = 0.05, method = “neuralnet”). ent parameters did not affect the results. We used a recombination rate The resulting confusion matrix is shown in Supplementary Table 3. of 9.40 × 10−9 per site per generation and 19 years generation time with We also performed cross-validation and model comparison for the the median fragment length to estimate introgression time. Decod- five demographic models: (A) null demography, (B) ghost gene flow ing, that is assigning internal and external states to specific genomic into the eastern common ancestor, (C) ghost gene flow into the west- regions, was done with the script provided with the repository. Puta- ern common ancestor, (D) revised model of ghost gene flow into the tive external fragments were filtered for posterior probabilities of 0.9 eastern common ancestor and (E) revised model of ghost gene flow (lenient) or 0.95 (strict) and required to contain at least five private into the western common ancestor, where we still observed model B positions. We also conducted a performance analysis of hmmix on having the highest support (Supplementary Table 4 and Supplemen- introgressed fragments in simulations of either model B or the ‘worst’ tary Section 3.4). model B, with results similar to those for S* (Supplementary Table 8). Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 For performance testing the hmmix approach, we used the posterior Plots were created with ggplot2 (ref. 65), circlize66 and pheatmap probabilities estimated by hmmix as thresholds, which range from (https://github.com/raivokolde/pheatmap). Genomic ranges were 0 to 0.9999. analysed with the GenomicRanges package67. We note that only hmmix could be used to infer archaic intro- gressed fragments on the X chromosome, due to the lack of a gorilla Adaptive introgression demographic model for the sex chromosomes. To explore signatures of adaptive introgression within eastern goril- las, we applied the genome-wide scan VolcanoFinder34. To do so, we Exploring introgressed regions polarized the data to two outgroups. First, we polarized the human To obtain a consensus set of putative introgressed regions, we over- reference allele using the rhesus macaque allele and subsequently lapped the autosomal outlier regions inferred under the two methods polarized the gorilla genotypes by this polarized allele representing within each eastern gorilla. For this overlap, we calculated the percent- the ancestral state. To obtain the allele frequency input files per chro- age of overlapping base pairs, considering in turn each S* confidence mosome, we then filtered our data to only eastern gorilla genotypes interval (95%, 99% and 99.5%) and with and without a 40 kb length cutoff at biallelic sites and also filtered out sites with multiple ancestral for hmmix regions identified under the strict threshold. Imposing alleles (where polarization would be uncertain) and sites of reference a 40 kb length cutoff retains 76.7% of the total strict hmmix regions homozygotes. The second input file required is an empirical unnormal- (Supplementary Tables 9 and 11). We consider the intersect of the S* ized site frequency spectrum (SFS), which we generated by obtaining 99% outliers with the strict hmmix autosomal outliers, as our putative the unfolded SFS, normalizing so all site categories sum to 1 and then introgressed regions of high confidence. To determine whether the filtering out the first category (the 0 entry). We called VolcanoFinder overlap obtained differed from random expectation we generated specifying ‘-big 1000, D = −1, P = 1, Model = 1’. For computational effi- intersections of random regions, of equivalent distribution to the ciency, we performed the VolcanoFinder scan in blocks, whereby each empirical data, for 100 iterations. chromosome was split into blocks of approximately equal numbers As a proxy for gene density we calculated the proportion of of base pairs. We placed a test site every 1,000 bp (-big 1000). We set protein-coding base pairs within these regions of high confidence. D to −1, so VolcanoFinder iteratively tested a grid of values for genetic As above, we compared this to the proportion of protein-coding base distance internally and selected the value that maximizes the likeli- pairs within 100 iterations of random genomic regions, of equal length hood ratio34. We set P to 1 as our input data were polarized. We used distribution as the putative introgressed regions within each eastern Model = 1, following procedures applied to human data34, as well as gorilla. We calculated pairwise nucleotide differences between indi- non-human species68,69. viduals in the putative introgressed regions and in random genomic We took the 95% outliers of composite likelihood ratio scores regions of equal length distribution and sufficient callable sites. This calculated from VolcanoFinder and intersected these regions with our was conducted for three comparisons: (1) among eastern gorillas, (2) putative introgressed regions (identified above), to obtain putative among western gorillas and (3) between eastern and western gorillas. adaptive introgressed targets. To explore potential functional con- To assess mutational tolerance, we used GERP, SIFT, PolyPhen-2 sequences, we assessed which genes and which mutations fall within and LINSIGHT scores58–61. We calculated the proportion of high impact the putative adaptive introgressed targets, using the Variant Effect sites for GERP, SIFT and PolyPhen-2 scores and the mean LINSIGHT Predictor annotation (v.83)70. score within putative introgressed regions and random regions of equal length distribution and sufficient callable sites. To explore the impact Reporting summary of introgression on regulatory elements, we calculated the proportion Further information on research design is available in the Nature Port- of regulatory base pairs using gorilla-defined regulatory element folio Reporting Summary linked to this article. annotations33, within putative introgressed and random regions of equivalent length and callability. This was assessed globally and per Data availability regulatory element type, considering poised, strong and weak, enhanc- The six newly sequenced eastern gorilla samples are publicly available ers and promoters. in the European Nucleotide Archive (ENA) under the project number: We further explored our putative introgressed regions of high PRJEB12821. ENA accession numbers for all samples used in this study confidence using PCA (Supplementary Fig. 13). This was generated are given in Extended Data Table 1. The human reference genome (hg19) using the biallelic sites in our putative introgressed regions. For com- and the rhesus macaque reference genome (Mmul10/rheMac10) were parison, we also generated PCAs of one random set of random regions, downloaded from https://hgdownload.soe.ucsc.edu/goldenPath/. with equal length distribution of random regions as the putative intro- Precalculated GERP scores for hg19 were accessed from http://mendel. gressed regions per eastern gorilla. The PCAs in Fig. 1 were generated stanford.edu/SidowLab/downloads/gerp/ and LINSIGHT scores for using biallelic SNPs of random genomic regions of equivalent length hg19 from https://rdrr.io/github/rcastelo/GenomicScores/src/inst/ distribution to the putative introgressed regions of GBB Bwiruka. This scripts/make-data_linsight.UCSC.hg19.R. sample of random genomic regions is representative of the whole genome. All PCAs were generated with the R package adegenet62. We Code availability generated phylogenetic trees of our putative introgressed regions and Scripts used for data analysis are available on Github under https:// one random replicate (Supplementary Fig. 14), using the ‘K80’ model github.com/h-pawar/gor_ghost_introg. of nucleotide substitution, using the adegenet package63. Haplotype networks were drawn using pegas64. References We localized introgression deserts by screening 1 Mb 1. Grubb, P. et al. Assessment of the diversity of African primates. non-overlapping windows (bins) spanning the genome. We filtered Int. J. Primatol. 24, 1301–1357 (2003). out bins overlapping centromeres and those at the end of each chromo- 2. Gray, M. et al. Genetic census reveals increased but uneven some which were <1 Mb in size. Per bin we calculated the frequency of growth of a critically endangered mountain gorilla population. putative introgressed regions falling within the bin, for each eastern Biol. Conserv. 158, 230–238 (2013). gorilla. We also calculated data coverage of the bins and filtered by 3. Plumptre, A. J. et al. Catastrophic decline of world’s largest mean callable proportion >0.5. Deserts hence constitute bins where primate: 80% loss of Grauer’s gorilla (Gorilla beringei graueri) no eastern gorilla carried a putative introgressed region and which population justifies critically endangered status. PLoS ONE 11, had a reasonable number of callable sites. e0162697 (2016). Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 4. Maisels, F., Williamson, E. A. & Bergl, R. IUCN Red List of 26. Tocheri, M. W. et al. The evolutionary origin and population history Threatened Species: Gorilla gorilla (IUCN, 2016); https://www. of the Grauer gorilla. Am. J. Phys. Anthropol. 159, S4–S18 (2016). iucnredlist.org/species/9404/136250858 27. Roy, J., Gray, M., Stoinski, T., Robbins, M. M. & Vigilant, L. 5. Fünfstück, T. & Vigilant, L. The geographic distribution of genetic Fine-scale genetic structure analyses suggest further male than diversity within gorillas. Am. J. Primatol. 77, 974–985 (2015). female dispersal in mountain gorillas. BMC Ecol. 14, 21 (2014). 6. Bergl, R. A. & Vigilant, L. Genetic analysis reveals population 28. Scally, A. et al. Insights into hominid evolution from the gorilla structure and recent migration within the highly fragmented genome sequence. Nature 483, 169–175 (2012). range of the Cross River gorilla (Gorilla gorilla diehli). Mol. Ecol. 16, 29. Vernot, B. et al. Excavating Neandertal and Denisovan DNA 501–516 (2007). from the genomes of melanesian individuals. Science 352, 7. Prado-Martinez, J. et al. Great ape genetic diversity and 235–239 (2016). population history. Nature 499, 471–475 (2013). 30. Skov, L. et al. Detecting archaic introgression using an unadmixed 8. Xue, Y. et al. Mountain gorilla genomes reveal the impact of outgroup. PLoS Genet. 14, e1007641 (2018). long-term population decline and inbreeding. Science 348, 31. Huang, X., Kruisz, P. & Kuhlwilm, M. sstar: a Python package for 242–245 (2015). detecting archaic introgression from population genetic data 9. Thalmann, O., Fischer, A., Lankester, F., Pääbo, S. & Vigilant, L. The with S*. Mol. Biol. Evol. 39, msac12 (2022). complex evolutionary history of gorillas: insights from genomic 32. Sankararaman, S. et al. The genomic landscape of Neanderthal data. Mol. Biol. Evol. 24, 146–158 (2007). ancestry in present-day humans. Nature 507, 354–357 (2014). 10. Sarmiento, E. E., Butynski, T. M. & Kalina, J. Gorillas of 33. García-Pérez, R. et al. Epigenomic profiling of primate Bwindi-Impenetrable Forest and the Virunga Volcanoes: lymphoblastoid cell lines reveals the evolutionary patterns taxonomic implications of morphological and ecological of epigenetic activities in gene regulatory architectures. Nat. differences. Am. J. Primatol. 40, 1–21 (1996). Commun. 12, 3116 (2021). 11. Robbins, M. M. & Robbins, A. M. Variation in the social 34. Setter, D. et al. VolcanoFinder: genomic scans for adaptive organization of gorillas: life history and socioecological introgression. PLoS Genet. 16, e1008867 (2020). perspectives. Evol. Anthropol. 27, 218–233 (2018). 35. Besenbacher, S., Hvilsom, C., Marques-Bonet, T., Mailund, T. 12. Kuhlwilm, M. et al. Evolution and demography of the great apes. & Schierup, M. H. Direct estimation of mutations in great apes Curr. Opin. Genet. Dev. 41, 124–129 (2016). reconciles phylogenetic dating. Nat. Ecol. Evol. 3, 286–292 (2019). 13. McManus, K. F. et al. Inference of gorilla demographic and 36. Baas, P. et al. Population-level assessment of genetic diversity and selective history from whole-genome sequence data. Mol. Biol. habitat fragmentation in critically endangered Grauer’s gorillas. Evol. 32, 600–612 (2015). Am. J. Phys. Anthropol. 165, 565–575 (2018). 14. Mailund, T. et al. A new isolation with migration model along comp­ 37. Fontsere, C. et al. Population dynamics and genetic connectivity lete genomes infers very different divergence processes among in recent chimpanzee history. Cell Genom. 2, 100133 (2022). closely related great ape species. PLoS Genet. 8, e1003125 (2012). 38. Anthony, N. M. et al. The role of Pleistocene refugia and rivers in 15. Becquet, C. & Przeworski, M. A new approach to estimate shaping gorilla genetic diversity in central Africa. Proc. Natl Acad. parameters of speciation models with application to apes. Sci. USA 104, 20432–20436 (2007). Genome Res. 17, 1505–1519 (2007). 39. Das, R. et al. Complete mitochondrial genome sequence of the 16. Tricou, T., Tannier, E. & de Vienne, D. M. Ghost lineages can Eastern gorilla (Gorilla beringei) and implications for African ape invalidate or even reverse findings regarding gene flow. PLoS Biol. biogeography. J. Hered. 105, 752–761 (2014). 20, e3001776 (2022). 40. Durvasula, A. & Sankararaman, S. Recovering signals of 17. Pang, X.-X. & Zhang, D.-Y. Impact of ghost introgression on ghost archaic introgression in African populations. Sci. Adv. 6, coalescent-based species tree inference and estimation of eaax5097 (2020). divergence time. Syst. Biol. 72, 35–49 (2022). 41. Zhang, D. et al. ‘Ghost introgression’ as a cause of deep 18. Fontsere, C., de Manuel, M., Marques-Bonet, T. & Kuhlwilm, M. mitochondrial divergence in a bird species complex. Admixture in mammals and how to understand its functional Mol. Biol. Evol. 36, 2375–2386 (2019). implications: on the abundance of gene flow in mammalian 42. Rocha, J. L. et al. African climate and geomorphology drive species, its impact on the genome, and roads into a functional evolution and ghost introgression in sable antelope. Mol. Ecol. 31, understanding. Bioessays 41, e1900123 (2019). 2968–2984 (2022). 19. Plagnol, V. & Wall, J. D. Possible ancestral structure in human 43. Skov, L. et al. Extraordinary selection on the human X populations. PLoS Genet. 2, e105 (2006). chromosome associated with archaic admixture. Cell Genom. 3, 20. Vernot, B. & Akey, J. M. Resurrecting surviving Neandertal lineages 100274 (2023). from modern human genomes. Science 343, 1017–1021 (2014). 44. Chevy, E. T., Huerta-Sánchez, E. & Ramachandran, S. Integrating 21. Green, R. E. et al. A draft sequence of the Neandertal genome. sex-bias into studies of archaic admixture on chromosome X. Science 328, 710–722 (2010). Preprint at bioRxiv https://doi.org/10.1101/2022.08.30.505789 22. Kuhlwilm, M., Han, S., Sousa, V. C., Excoffier, L. & Marques-Bonet, (2022). T. Ancient admixture from an extinct ape lineage into bonobos. 45. Harcourt, A. H., Stewart, K. S. & Fossey, D. Male emigration Nat. Ecol. Evol. 3, 957–965 (2019). and female transfer in wild mountain gorilla. Nature 263, 23. Cooke, N. P. & Nakagome, S. Fine-tuning of approximate Bayesian 226–227 (1976). computation for human population genomics. Curr. Opin. Genet. 46. Vigilant, L. et al. Reproductive competition and inbreeding Dev. 53, 60–69 (2018). avoidance in a primate species with habitual female dispersal. 24. van der Valk, T., Díez-Del-Molino, D., Marques-Bonet, T., Behav. Ecol. Sociobiol. 69, 1163–1172 (2015). Guschanski, K. & Dalén, L. Historical genomes reveal the genomic 47. Zhang, X., Kim, B., Lohmueller, K. E. & Huerta-Sánchez, E. The consequences of recent population decline in eastern gorillas. impact of recessive deleterious variation on signals of adaptive Curr. Biol. 29, 165–170 (2019). introgression in human populations. Genetics 215, 799–812 (2020). 25. van der Valk, T. et al. Significant loss of mitochondrial diversity 48. Di Pizio, A. & Niv, M. Y. Promiscuity and selectivity of bitter within the last century due to extinction of peripheral populations molecules and their receptors. Bioorg. Med. Chem. 23, in eastern gorillas. Sci. Rep. 8, 6551 (2018). 4082–4091 (2015). Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 49. Carulli, D., de Winter, F. & Verhaagen, J. Semaphorins in adult programme (grant no. 864203), PID2021-126004NB-100 (MINECO/ nervous system plasticity and disease. Front. Synaptic Neurosci. FEDER, UE), Secretaria d’Universitats i Recerca and CERCA Program 13, 672891 (2021). del Departament d’Economia i Coneixement de la Generalitat de 50. Witt, K. E. & Huerta-Sánchez, E. Convergent evolution in human Catalunya (GRC 2021 SGR 00177) to T.M.-B. H.P. was supported by a and domesticate adaptation to high-altitude environments. Formació de Personal Investigador fellowship from Generalitat de Philos. Trans. R. Soc. Lond. B 374, 20180235 (2019). Catalunya (FI_B100131). M.A.-E. was supported by a Formación de 51. Moraitou, M. et al. Ecology, not host phylogeny, shapes the Personal Investigador PRE2018-083966 from Ministerio de Ciencia, oral microbiome in closely related species. Mol. Biol. Evol. 39, Universidades e Investigación. C.T.-S., Y.X. and J.P.-M. were funded msac263 (2022). by Wellcome grant no. 098051. K.G. was supported by Swedish 52. de Manuel, M. et al. Chimpanzee genomic diversity reveals Research Council grant no. 2020-03398. J.L.K. received the María de ancient admixture with bonobos. Science 354, 477–481 (2016). Maeztu Mobility Fellowship. O.D. was supported by a John Templeton 53. Lindenbaum, P. JVarkit: java-based utilities for Bioinformatics. Foundation grant no. ID 62178. A.M.A. received funding from UCL’s Figshare https://doi.org/10.6084/m9.figshare.1425030.v1 (2015). Wellcome Trust ISSF3 award no. 204841/Z/16/Z. Q.A. is supported by 54. Csilléry, K., François, O. & Blum, M. G. B. abc: an R package for strategic funding from Monash University (STG-000114). approximate Bayesian computation (ABC). Methods Ecol. Evol. 3, 475–479 (2012). Author contributions 55. Hudson, R. R. Generating samples under a Wright–Fisher neutral T.M.-B. and M.K. conceived and conceptualized the study. H.P. model of genetic variation. Bioinformatics 18, 337–338 (2002). performed demographic modelling. Q.A. performed experiments. 56. Baumdicker, F. et al. Efficient ancestry and mutation simulation H.P., A.R., M.d.M., T.v.d.V., I.L., M.H., R.B., O.D., S.H., P.E.-C., J.P.-M. and with msprime 1.0. Genetics 220, iyab229 (2022). M.K. analysed data. X.H., O.L. and D.J. provided software. K.G., A.M.A., 57. Kelleher, J., Etheridge, A. M. & McVean, G. Efficient coalescent C.T.-S., Y.X., T.M.-B. and M.K. provided supervision. B.S., M.C., C.T.-S. simulation and genealogical analysis for large sample sizes. and Y.X. acquired samples and their documentation. M.A.-E. visualized PLoS Comput. Biol. 12, e1004842 (2016). results. S.C.-E., J.L.K. and O.E.C. provided comments. H.P., T.M.-B. and 58. Davydov, E. V. et al. Identifying a high fraction of the human M.K. wrote the paper with input from all authors. genome to be under selective constraint using GERP++. PLoS Comput. Biol. 6, e1001025 (2010). Competing interests 59. Kumar, P., Henikoff, S. & Ng, P. C. Predicting the effects of coding The authors declare no competing interests. non-synonymous variants on protein function using the SIFT algorithm. Nat. Protoc. 4, 1073–1081 (2009). Additional information 60. Adzhubei, I. A. et al. A method and server for predicting damaging Extended data is available for this paper at missense mutations. Nat. Methods 7, 248–249 (2010). https://doi.org/10.1038/s41559-023-02145-2. 61. Huang, Y.-F., Gulko, B. & Siepel, A. Fast, scalable prediction of deleterious noncoding variants from functional and population Supplementary information The online version genomic data. Nat. Genet. 49, 618–624 (2017). contains supplementary material available at 62. Jombart, T. & Ahmed, I. adegenet 1.3-1: new tools for the analysis https://doi.org/10.1038/s41559-023-02145-2. of genome-wide SNP data. Bioinformatics 27, 3070–3071 (2011). 63. Jombart, T. adegenet: a R package for the multivariate analysis of Correspondence and requests for materials should be addressed to genetic markers. Bioinformatics 24, 1403–1405 (2008). Tomas Marques-Bonet or Martin Kuhlwilm. 64. Paradis, E. pegas: an R package for population genetics with an integrated-modular approach. Bioinformatics 26, 419–420 (2010). Peer review information Nature Ecology & Evolution thanks Da-Yong 65. Wickham, H. ggplot2: Elegant Graphics for Data Analysis Zhang, Laurits Skov and Colin Brand for their contribution to the peer (Springer, 2016). review of this work. Peer reviewer reports are available. 66. Gu, Z., Gu, L., Eils, R., Schlesner, M. & Brors, B. circlize implements and enhances circular visualization in R. Bioinformatics 30, Reprints and permissions information is available at 2811–2812 (2014). www.nature.com/reprints. 67. Lawrence, M. et al. Software for computing and annotating genomic ranges. PLoS Comput. Biol. 9, e1003118 (2013). Publisher’s note Springer Nature remains neutral with regard to 68. Moest, M. et al. Selective sweeps on novel and introgressed jurisdictional claims in published maps and institutional affiliations. variation shape mimicry loci in a butterfly adaptive radiation. PLoS Biol. 18, e3000597 (2020). Open Access This article is licensed under a Creative Commons 69. Liu, S. et al. Demographic history and natural selection shape Attribution 4.0 International License, which permits use, sharing, patterns of deleterious mutation load and barriers to introgression adaptation, distribution and reproduction in any medium or format, across Populus genome. Mol. Biol. Evol. 39, msac008 (2022). as long as you give appropriate credit to the original author(s) and the 70. McLaren, W. et al. The ensembl variant effect predictor. Genome source, provide a link to the Creative Commons license, and indicate Biol. 17, 122 (2016). if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless Acknowledgements indicated otherwise in a credit line to the material. If material is not We thank D. Setter for valuable guidance in applying VolcanoFinder. included in the article’s Creative Commons license and your intended We thank the Uganda Wildlife Authority for the Gorilla monitoring and use is not permitted by statutory regulation or exceeds the permitted research permission. We are grateful to the Life Science Compute use, you will need to obtain permission directly from the copyright Cluster of the University of Vienna. This project has been funded holder. To view a copy of this license, visit http://creativecommons. by the Vienna Science and Technology Fund (WWTF) (grant no. org/licenses/by/4.0/. 10.47379/VRG20001) to M.K. and the European Research Council under the European Union’s Horizon 2020 research and innovation © The Author(s) 2023 Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 1 Institute of Evolutionary Biology (UPF-CSIC), PRBB, Barcelona, Spain. 2Department of Evolutionary Anthropology, University of Vienna, Vienna, Austria. 3 Human Evolution and Archaeological Sciences (HEAS), University of Vienna, Wien, Austria. 4Department of Bioinformatics and Genetics, Scilifelab, Swedish Museum of Natural History, Stockholm, Sweden. 5Centre for Palaeogenetics, Stockholm, Sweden. 6Institute of Cancer and Genomic Sciences, University of Birmingham, Dubai, United Arab Emirates. 7Integrative Genomics Lab, CIC bioGUNE—Centro de Investigación Cooperativa en Biociencias, Parque Científico Tecnológico de Bizkaia building 801A, Derio, Spain. 8Wellcome Sanger Institute, Hinxton, UK. 9Monash University Malaysia Genomics Facility, School of Science, Monash University Malaysia, Selangor Darul Ehsan, Malaysia. 10Department of Ecology and Evolutionary Biology, University of California, Santa Cruz, CA, USA. 11UCL Genetics Institute, Department of Genetics, Evolution and Environment, University College London, London, UK. 12 Animal Ecology, Department of Ecology and Genetics, Uppsala University, Uppsala, Sweden. 13Institute of Ecology and Evolution, School of Biological Sciences, University of Edinburgh, Edinburgh, UK. 14Science for Life Laboratory, Uppsala, Sweden. 15Gorilla Doctors, Kampala, Uganda. 16Gorilla Doctors, Karen C. Drayer Wildlife Health Center, One Health Institute, University of California Davis, School of Veterinary Medicine, Davis, CA, USA. 17Catalan Institution of Research and Advanced Studies (ICREA), Passeig de Lluís Companys, Barcelona, Spain. 18CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Barcelona, Spain. 19Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, Barcelona, Spain. 20These authors contributed equally: Tomas Marques-Bonet, Martin Kuhlwilm. e-mail: tomas.marques@upf.edu; martin.kuhlwilm@univie.ac.at Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 Gorilla Genotypes n=49: 12 MG, 9 ELG, 27 WLG, 1 CRG ABC-based Demographic Modelling 1. Parameter Inference Localise archaic Introgression Adaptive Introgression A Null Model B Ghost Eastern Ancestor S* hmmix VolcanoFinder C Ghost Western Ancestor 2. Model Comparison (A, B, C) Introgressed Regions Revised Models Adaptive Introgressed 1. Ghost Parameter Inference Functional Regions Sample all parameters from priors Consequences D Ghost Eastern Ancestor 1. Gene density E Ghost Western Ancestor 2. Pairwise nucleotide differences 2. Model Comparison 3. Introgression deserts (A, B, C, D, E) 4. Autosome: X ratio 5. Mutational conservation Validation 6. Regulatory elements 1. Power of S* and hmmix under model A 2. Power of S* and hmmix under model misspecification Extended Data Fig. 1 | Workflow of the main analyses. Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 A 14,363 10,414-23,298 965 kya 729kya-1.1mya 98,135 91,423-99,820 454 kya 5,325 448-456 kya 1,473-23,015 48,288 42,636-49,691 40 kya 1.52% 0.17-2.48% 22-45 kya 34 kya 0.80% 0.06-2.14% 14 kya 9,911 11-18 kya 2,919-19,163 21,982 8,003- 9 kya 24,686 115 7 kya 15-1,478 243 76-944 present MG EL WL CR 2,157 20,894 64,748 14,558 719-7,680 2,308-29,776 23,758-95,986 6,959-19,630 B 79,021 63,561-93,320 1.9 mya 1.5-3.0 mya 1.07 mya 0.17% 1.06-1.07 mya 0.09-0.48% 14,363 10,414-23,298 965 kya 729kya-1.1mya 5,325 98,135 1,473- 91,423-99,820 454 kya 23,015 448-456 kya 48,288 42,636-49,691 40 kya 22-45 kya 2.00% 0.79-2.46% 34 kya Ghost 0.80% 0.06-2.14% 18 kya 9,811 17-18 kya 5,786-17,669 21,982 8,003- 24,686 9 kya 115 7 kya 15-1478 243 76-944 present MG EL WL CR 2,157 20,894 64,748 14,558 719-7,680 2,308-29,776 23,758-95,986 6,959-19,630 Extended Data Fig. 2 | Demographic models A and C. A Null model of gorilla ghost gene flow to the western common ancestor, the posteriors indicate a small population history (only extant populations). 95% credible intervals are shown contribution to the common ancestor of all gorillas (consistent with ancestral for all parameters inferred. B Alternate model allowing the possibility of ghost substructure), rather than a defined pulse to the western common ancestor. In introgression into the common ancestor of western gorillas, resulted in a model darker colours are the parameters inferred under this alternate model with their of ancestral population structure being inferred. We note under a model of 95% credible intervals. Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 Extended Data Fig. 3 | Prior and posterior distributions for model B. gorillas. Red indicates the posterior distribution inferred with neural networks. Parameter distributions for all parameters inferred under the ABC model Black indicates the posterior distribution inferred under a rejection method. The allowing gene flow from a ghost lineage into the common ancestor of eastern dotted grey line indicates the prior distribution. Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 Extended Data Fig. 4 | Performance of S* and hmmix. Precision-recall curves consider the target population as alternately eastern lowland or mountain for the S* statistic as implemented in sstar31 and for hmmix. Main model refers gorillas, eg Main Model EL. Worst mis-specified is where we generate simulated to a model taking the weighted median posteriors from the ABC-based null data under the worst model but run the S* analysis using the ‘quantile’ or outlier demography presented herein (Extended Data Fig. 2A). Worst model refers to values inferred under the main model. Skov=hmmix method, EL=eastern lowland a model taking the maximum value of the 95% credible interval for all ancestral gorillas, M=mountain gorillas. Ne parameters from the ABC-based null demography. For the S* statistic we Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 Extended Data Table 1 | Information for gorillas analysed herein. 49 samples: 12 Gorilla beringei beringei (mountain gorillas), 9 Gorilla beringei graueri (eastern lowland gorillas), 1 Gorilla gorilla diehli (Cross River gorillas), 27 Gorilla gorilla gorilla (western lowland gorillas) Subspecies Individual SRA ID Sex Country or origin Project Mountain gorilla Bwiruka ERR2300765 F Uganda (Bwindi) this paper Mountain gorilla Imfura ERS168207 M Rwanda (Virunga) Xue et al. 2015 Mountain gorilla Kaboko ERS168410 M DRC (Virunga) Xue et al. 2015 Mountain gorilla Kahungye ERR2300762 F Uganda (Bwindi) this paper Mountain gorilla Katungi ERR2300763 F Uganda (Bwindi) this paper Mountain gorilla Maisha ERS525616 F DRC (Virunga) Xue et al. 2015 Mountain gorilla Nyamunwa ERR2300764 F Uganda (Bwindi) this paper Mountain gorilla Semehe ERR2300766 F Uganda (Bwindi) this paper Mountain gorilla Tuck ERS168204 F Rwanda (Virunga) Xue et al. 2015 Mountain gorilla Turimaso ERS525618 F Rwanda (Virunga) Xue et al. 2015 Mountain gorilla Umurimo ERS525617 F Rwanda (Virunga) Xue et al. 2015 Mountain gorilla Zirikana ERS168174 M Rwanda (Virunga) Xue et al. 2015 Eastern lowland gorilla 9732_Mkubwa SRS396825 M DRC (Tulakwa) Prado-Martinez et al. 2013 Eastern lowland gorilla A929_Kaisi SRS396605 M DRC (Walikale) Prado-Martinez et al. 2013 Eastern lowland gorilla A967_Victoria SRS396876 F DRC Prado-Martinez et al. 2013 Eastern lowland gorilla Dunia ERS525621 F DRC (Walikale) Xue et al. 2015 Eastern lowland gorilla Itebero ERS168205 F DRC (Kahuzi-Biega) Xue et al. 2015 Eastern lowland gorilla Mukokya ERR2300767 M DRC (Mount Tshiaberimu) this paper Eastern lowland gorilla Ntabwoba ERS168206 M DRC (Walikale) Xue et al. 2015 Eastern lowland gorilla Pinga ERS525620 F DRC Xue et al. 2015 Eastern lowland gorilla Tumani ERS525619 F DRC (Walikale) Xue et al. 2015 Cross river gorilla B646_Nyango SRS396855 F West Africa Prado-Martinez et al. 2013 Western lowland gorilla 9749_Kowali SRS396819 F unknown Prado-Martinez et al. 2013 Western lowland gorilla 9750_Azizi SRS396821 M Cameroon Prado-Martinez et al. 2013 Western lowland gorilla 9751_Bulera SRS396820 F Cameroon Prado-Martinez et al. 2013 Western lowland gorilla 9752_Suzie SRS394796 F unknown Prado-Martinez et al. 2013 Western lowland gorilla 9753_Kokomo SRS396849 F unknown Prado-Martinez et al. 2013 Western lowland gorilla A930_Sandra SRS396606 F Cameroon Prado-Martinez et al. 2013 Western lowland gorilla A931_Banjo SRS396826 M Cameroon Prado-Martinez et al. 2013 Western lowland gorilla A932_Mimi SRS396827 F Cameroon Prado-Martinez et al. 2013 Western lowland gorilla A933_Dian SRS396828 F Cameroon Prado-Martinez et al. 2013 Western lowland gorilla A934_Delphi SRS396829 F Congo Prado-Martinez et al. 2013 Western lowland gorilla A936_Coco SRS396830 F Equatorial Guinea Prado-Martinez et al. 2013 Western lowland gorilla A937_Kolo SRS396831 F Cameroon Prado-Martinez et al. 2013 Western lowland gorilla A962_Amani SRS396847 F unknown Prado-Martinez et al. 2013 Western lowland gorilla B642_Akiba_Beri SRS396852 F Cameroon Prado-Martinez et al. 2013 Western lowland gorilla B643_Choomba SRS396853 F West Africa Prado-Martinez et al. 2013 Western lowland gorilla B644_Paki SRS396854 F West Africa Prado-Martinez et al. 2013 Western lowland gorilla B647_Anthal SRS396856 F West Africa Prado-Martinez et al. 2013 Western lowland gorilla B650_Katie SRS396857 F West Africa Prado-Martinez et al. 2013 Western lowland gorilla KB3782_Vila SRS396870 F Congo Prado-Martinez et al. 2013 Western lowland gorilla KB3784_Dolly SRS396873 F Congo Prado-Martinez et al. 2013 Western lowland gorilla KB4986_Katie SRS396874 F unknown Prado-Martinez et al. 2013 Western lowland gorilla KB5792_Carolyn SRS396868 F Congo Prado-Martinez et al. 2013 Western lowland gorilla KB5852_Helen SRS396871 F Cameroon Prado-Martinez et al. 2013 Western lowland gorilla KB6039_Oko SRS396872 F unknown Prado-Martinez et al. 2013 Western lowland gorilla KB7973_Porta SRS396869 F unknown Prado-Martinez et al. 2013 Western lowland gorilla X00108_Abe SRS396850 M unknown Prado-Martinez et al. 2013 Western lowland gorilla X00109_Tzambo SRS396851 M unknown Prado-Martinez et al. 2013 SRA ID = Short Read Archive identifier; F = Female; M = Male; DRC = Democratic Republic of the Congo Nature Ecology & Evolution Article https://doi.org/10.1038/s41559-023-02145-2 Extended Data Table 2 | Regions and genes with signatures of putative adaptive introgression Eastern lowland gorillas Region Genes Minimum LR Maximum LR 12:11051000- TAS2R14, 24.43 246.15 11122000 PRH2, TAS2R13 5:9061000-9111000, SEMA5A 26.75 27.97 5:9315000-9503000 1:191807000- RP11-541F9.2 25.92 35.06 191925000 3:124338000- KALRN 24.07 24.07 124414000 10:113899000- GPAM 24.99 30.16 113999000 7:26805000- SKAP2, HOXA1- 24.65 36.19 27214000 10 Mountain gorillas 12:11075000- TAS2R14, 24.43277 246.15161 11121000 PRH2, TAS2R13 5:9038000-9232000, 5:9249000-9414000, SEMA5A 26.75713 27.97717 5:9436000-9503000 4:69349546- TMPRSS11E 24.0434 24.0434 69349547 3:124318000- KALRN 24.07903 24.07903 124414000 LR = Likelihood Ratio Nature Ecology & Evolution nature portfolio | reporting summary Corresponding author(s): Martin Kuhlwilm Last updated by author(s): Jun 1, 2023 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist. Statistics For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section. n/a Confirmed The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one- or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section. A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted Give P values as exact values whenever suitable. For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated Our web collection on statistics for biologists contains articles on many of the points above. Software and code Policy information about availability of computer code Data collection no software was used Data analysis bcftools 1.6/1.12/1.14, bedtools 2.26.0, bedops, gcc 6.3.0, hmmix, jvarkit, ms, msprime, openssl 1.0.2q, perl, python 2.7.11/2.7.14/2.7.17, R 3.2.0/3.4.0/4.0.1, samtools 1.0, sstar, tabix 0.2.6, Variant Effect Predictor 83, vcftools 0.1.15, VolcanoFinder 1.0, xz 5.2.2, zlib 1.2.8, WebGestaltR For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information. Data March 2021 Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: - Accession codes, unique identifiers, or web links for publicly available datasets - A description of any restrictions on data availability - For clinical datasets or third party data, please ensure that the statement adheres to our policy The six newly sequenced eastern gorilla samples are publicly available in the European Nucleotide Archive (ENA) under the project number: PRJEB12821. ENA accession numbers for all samples used in this study are given in Table S1. The human reference genome (hg19) and the rhesus macaque reference genome 1 (Mmul10/rheMac10) were downloaded from https://hgdownload.soe.ucsc.edu/goldenPath/. Pre-calculated GERP scores for hg19 were accessed from (http:// mendel.stanford.edu/SidowLab/downloads/gerp/) and LINSIGHT scores for hg19 from (https://rdrr.io/github/rcastelo/GenomicScores/src/inst/scripts/make- nature portfolio | reporting summary data_linsight.UCSC.hg19.R). Human research participants Policy information about studies involving human research participants and Sex and Gender in Research. Reporting on sex and gender Not applicable Population characteristics Not applicable Recruitment Not applicable Ethics oversight Not applicable Note that full information on the approval of the study protocol must also be provided in the manuscript. Field-specific reporting Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection. Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf Life sciences study design All studies must disclose on these points even when the disclosure is negative. Sample size 49 whole genomes of endangered wild species. No more samples could be obtained. Data exclusions One sample (Nkuhene) was excluded due to very low quality, with 80% of read duplicates and 2X average coverage. Replication Not applicable, as the 49 genomes are independent individuals from the wild, as representatives of their respective populations. Randomization Not applicable, as these are wild individuals, and no experimental conditions apply. Blinding Not applicable, as no experimental conditions apply. Reporting for specific materials, systems and methods We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response. Materials & experimental systems Methods n/a Involved in the study n/a Involved in the study Antibodies ChIP-seq Eukaryotic cell lines Flow cytometry Palaeontology and archaeology MRI-based neuroimaging Animals and other organisms Clinical data March 2021 Dual use research of concern 2