Report Historical Genomes Reveal the Genomic Consequences of Recent Population Decline in Eastern Gorillas Highlights Authors d Gorilla genomes from museum samples enabled Tom van der Valk, measurement of temporal genomic erosion David Dı́ez-del-Molino, Tomas Marques-Bonet, d Genetic diversity declined and inbreeding increased in n Katerina Guschanski, Love Dale Grauer’s gorillas d The frequency of deleterious mutations increased Correspondence significantly in Grauer’s gorillas katerina.guschanski@ebc.uu.se (K.G.), love.dalen@nrm.se (L.D.) d No temporal changes in genomic parameters were observed in the mountain gorilla In Brief Van der Valk et al. investigate the genomic consequences of recent population decline using 100-year-old museum specimens. They find a reduction in genetic diversity, increase in inbreeding, and increase in deleterious mutations in the critically endangered Grauer’s gorillas. van der Valk et al., 2019, Current Biology 29, 165–170 January 7, 2019 ª 2018 Elsevier Ltd. https://doi.org/10.1016/j.cub.2018.11.055 Current Biology Report Historical Genomes Reveal the Genomic Consequences of Recent Population Decline in Eastern Gorillas Tom van der Valk,1 David Dı́ez-del-Molino,2 Tomas Marques-Bonet,3,4,5,6 Katerina Guschanski,1,7,8,9,* n2,7,8,* and Love Dale 1Department of Ecology and Genetics (Animal ecology), Uppsala University, 752 36 Uppsala, Sweden 2Department of Bioinformatics and Genetics, Swedish Museum of Natural History, Box 50007, 10405 Stockholm, Sweden 3Institute of Evolutionary Biology (UPF-CSIC), PRBB, Dr. Aiguader 88, 08003 Barcelona, Spain 4Catalan Institution of Research and Advanced Studies (ICREA), Passeig de Lluı́s Companys, 23, 08010, Barcelona, Spain 5CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Baldiri i Reixac 4, 08028 Barcelona, Spain 6Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, c/ Columnes s/n, 08193 Cerdanyola del Vallès, Barcelona, Spain 7Senior author 8These authors contributed equally 9Lead Contact *Correspondence: katerina.guschanski@ebc.uu.se (K.G.), love.dalen@nrm.se (L.D.) https://doi.org/10.1016/j.cub.2018.11.055 SUMMARY may increase the probability that deleterious alleles drift to fixa- tion, further lowering the fitness of the population [6]. Maintaining Many endangered species have experienced severe genetic diversity is thus of considerable importance for species population declines within the last centuries [1, 2]. conservation [7]. However, the genomic consequences of rapid However, despite concerns about negative fitness population decline, as faced by many species today, remain effects resulting from increased genetic drift and largely unexplored due to limited empirical data [9]. inbreeding, there is a lack of empirical data on In this study, we directly quantify genomic changes within the last century in the two eastern gorillas: the critically endangered genomic changes in conjunction with such declines Grauer’s (Gorilla beringei graueri) and the endangered mountain [3–7]. Here, we use whole genomes recovered from gorillas (Gorilla beringei beringei) [10] by sequencing whole ge- century-old historical museum specimens to quan- nomes from historical specimens collected up to ca. 100 years tify the genomic consequences of small population ago and comparing them to present-day (modern) genomes size in the critically endangered Grauer’s and endan- [11, 12] (Figure 1A). After low-depth sequencing of 59 historical gered mountain gorillas. We find a reduction of gorilla specimens (see Data S1), we selected samples from unre- genetic diversity and increase in inbreeding and ge- lated adult individuals displaying high endogenous DNA content netic load in the Grauer’s gorilla, which experienced and DNA quantity for uracil-DNA glycosylase (UDG) treatment severe population declines in recent decades. In and subsequent deep genome sequencing. Following whole- contrast, the small but relatively stable mountain genome sequencing, we collapsed the obtained forward and gorilla population has experienced little genomic reverse high-quality sequencing reads generated on the Illumina HiSeq X platform and aligned them to the gorilla reference change during the last century. These results sug- genome [13]. For seven Grauer’s and four mountain gorilla sam- gest that species histories as well as the rate of de- ples collected between 1910 and 1962 (median: 1923), we ob- mographic change may influence how population tained adequate coverage (3.1–10.8 X) to infer genomic changes declines affect genome diversity. through time (see Data S1). We also included published modern genomes from eight Grauer’s, seven mountain, and 17 western RESULTS AND DISCUSSION lowland gorillas (Gorilla gorilla gorilla) (see Data S1) [11], the last to be used as a comparative dataset. In the last centuries, continuously increasing anthropogenic To study changes in genetic diversity over time, we called ge- pressures have dramatically accelerated wild animal population notype likelihoods and variants across all historical and modern losses, especially affecting large-bodied mammals [1, 2]. These samples. We filtered low-quality sites, recalibrated quality population declines can have severe consequences for long- scores accounting for post-mortem DNA damage, and removed term species survival due to environmental and demographic cytosine-phosphate-guanine (CpG) sites to limit possible biases stochastic events [3], the negative effects of reduced genetic di- from remaining DNA damage (Figures S1A and S1B) [14]. versity on fertility [4], resistance to infectious diseases [8], and Using principal-component and admixture analyses based on adaptability to changing environments [5]. Moreover, natural se- genotype likelihoods, we found indications for the presence of lection is thought to be less efficient in small populations, which population structure in historical Grauer’s gorillas, corroborating Current Biology 29, 165–170, January 7, 2019 ª 2018 Elsevier Ltd. 165 Figure 1. Geography and Genetic Structure (A) Sampling location of historical and previously published modern samples. Two modern Grauer’s gorilla individuals are not shown on the map due to un- certainty about their exact geographic origin. y depicts samples from extinct Grauer’s gorilla populations. (B) Principal-component analysis of eastern gorilla genomes. See also Figure S2. previous results based on modern data (Figures 1B and S2A– is driven by the proportion of the genome contained in very S2E) [12, 15]. Geography could not explain the apparent struc- long ROH (2.5–10 Mb, Figures 2B–2C and S3), which increased ture, as individuals from the extremes of the species distribution by 24% in modern Grauer’s gorilla genomes (Figures 2B and 2C). range were found to belong to the same genomic subcluster Such long tracts are likely to arise from mating between closely (Figures S2A and S2B). Next, we estimated the number of differ- related individuals less than ten generations ago [16]. In regions ences between chromosome pairs to assess intra-individual ge- of the genome outside of ROH, heterozygosity levels were similar netic diversity (heterozygosity) after downsampling modern between historical and modern Grauer’s gorillas (Figure 2A), indi- samples to a similar coverage as the historical genomes to cating that recent inbreeding is the main cause of the observed control for coverage biases (Figure S1C). We found a decline in reduction in genome-wide genetic diversity. The proportion of autosomal heterozygosity in Grauer’s gorillas, with modern indi- the genome contained in long ROH increased by 6.9% in modern viduals containing on average 20% fewer heterozygous sites mountain gorillas; however, the difference with the historical than historical individuals (Figures 2A and S1D). We also samples was not statistically significant (Figures 2B and S3). observed lower genetic diversity in modern mountain gorillas; To test whether population decline and small population size however, the difference compared to the historical individuals led to an accumulation of deleterious mutations in eastern was small (3%) and not statistically significant (Figures 2A and gorillas, we used three independent and complementary ap- S1D). We obtained similar results independent of the method proaches to characterize genetic load on genome-wide, gene, used to call genomic variants (Figure S1D), and we can exclude and protein level. First, we estimated genetic load for each indi- population substructure, difference in sample sizes, and retained vidual by measuring the proportion of homozygous-derived mu- post-mortem DNA damage as possible explanations (Figures tations in regions of the genome that are conserved across a S1C and S2F–S2H). Estimating diversity at population level whole-genome alignment of 20 vertebrates and thus can be as the total number of variable sites per population also sup- assumed to be mostly deleterious [17]. Although the differences ported the observed reduction in genetic diversity in the in genetic load between modern and historical eastern gorilla in- modern Grauer’s gorillas compared to the historical samples dividuals were not significant, a larger proportion of modern ge- (Figure S1E). nomes of both Grauer’s and mountain gorillas fell into the upper Next, we inferred the effects of the recent population decline range of the distribution (Figure 3A). Second, we assessed the on inbreeding levels by identifying tracts of within-individual effects of DNA sequence variants on protein-coding genes and chromosomal sequence sharing (runs of homozygosity [ROH]) found a significant increase in frequency of both missense and using a Hidden Markov Model. Whereas an average of 35.1% loss-of-function (LoF) variants in modern Grauer’s gorillas and 33.9% of the historical Grauer’s and mountain gorilla ge- compared to the historical individuals, consistent with higher ge- nomes, respectively, consisted of ROH above 100 kb, this pro- netic load in the modern population (Figure 3B). In contrast, our portion increased to 39.2% in modern Grauer’s and to 36.3% results suggest that the frequency of missense and LoF variants in modern mountain gorillas (Figures 2B and S3). This change has remained stable in the mountain gorillas during the last 166 Current Biology 29, 165–170, January 7, 2019 Figure 2. Temporal Genetic Diversity and Inbreeding (A) Genome-wide heterozygosity of gorilla in- dividuals (solid outlines) and heterozygosity excluding runs of homozygosity (ROH, dotted outlines). Each black dot represents one individ- ual. **p < 0.01; NS, not significant. (B) Fraction of the genomes in ROH above 100 kb (open bars), and long ROH between 2.5 and 10 Mb (solid bars). Error bars, ±1 SD. *p < 0.05, **p < 0.01; NS, not significant. (C) Cumulative proportion of the genomes con- tained in ROH below the length displayed on the x axis. Thin lines represent individual genomes, whereas thick lines show population averages. Trajectories of modern and historical genomes start diverging at long ROH (>1 Mb). See also Figures S1–S3. tween species are dominated by long- term demographic processes. The observed decrease in genome-wide di- versity and increase in putatively delete- rious variants in eastern gorillas (Figures 2 and 3) is striking given the short time period spanned by our study, corre- sponding to 4–5 gorilla generations. Gra- uer’s gorillas have been more severely affected than mountain gorillas, which we hypothesize can be attributed to their contrasting demographic histories. Eastern gorillas experienced contin- uous population decline for the last 100,000 years [12], before diverging into century (Figure 3B). Finally, we predicted the possible delete- Grauer’s and mountain gorillas, ca. 10,000 years ago [21]. Sub- rious impact of each amino acid substitution on the protein sequently, Grauer’s gorillas went through a period of population function using physical considerations [18]. We observed a sig- growth and range expansion 5,000–10,000 years ago [22, 23]. nificant increase in the frequency of amino acid substitutions This demographic expansion may have led not only to a classified as ‘‘probably damaging’’ in the modern Grauer’s gorilla historically higher genetic diversity in Grauer’s compared to individuals, whereas no significant change was observed be- mountain gorillas (Figure 2) but also to a higher number of tween historical and modern mountain gorillas (Figure 3C). low-frequency deleterious mutations in the Grauer’s gorilla Genes affected by LoF mutations with the highest frequency in- population [24, 25]. As Grauer’s gorillas went through a severe crease in the modern Grauer’s gorillas were associated with func- population decline of 80% in the last 20 years to less than tions related to immunity and methylation (see Data S2 and S3), 4,000 individuals today [26], these deleterious mutations appear possibly indicating a decline in pathogen resistance during recent to have increased in frequency, likely due to increased drift and decades. LoF variants with the highest frequency increase in inbreeding (Figure 2), thus leading to a more pronounced change modern mountain gorillas were enriched for genes related to in genetic load in Grauer’s compared to mountain gorillas (Fig- Sertoli cell differentiation, thus possibly affecting male fertility ure 3). In contrast, the population size of mountain gorillas from (see Data S2 and S3). In addition, both Grauer’s and mountain the Virunga Massif population has likely remained small since gorillas show several high-frequency LoF variants that are related their divergence from the Grauer’s gorillas [22, 23]. Numbering to malformation of fingers and toes, corroborating the reports of fewer than 1,000 individuals at least since the late 1950s [27], syndactyly in eastern gorillas (see Data S3) [19]. mountain gorillas experienced a population low of 250 individ- By analyzing complete genomes from museum-preserved his- uals in the 1980s [28] but recovered to 450 individuals in 2013 torical specimens, we have identified and quantified the genomic [29]. However, we do not detect significant differences in genetic consequences of population decline in eastern gorillas during diversity, inbreeding, and genetic load between historical and the last century. Comparisons to the modern western lowland modern mountain gorilla samples. This might be due to both gorillas (estimated population size: 360,000 individuals), which purging of genetic load as a result of the continuously small pop- diverged from the eastern gorillas around 150,000 years ago ulation size in mountain gorillas during the last 10,000 years as [12, 20], demonstrate that differences in genetic diversity be- well as the rapid population recovery in recent decades thanks Current Biology 29, 165–170, January 7, 2019 167 Figure 3. Functional Consequences of Pop- ulation Decline in Eastern Gorillas (A) Relative genetic load measured as the fraction of homozygous-derived sites in evolutionary- constrained regions of the genome. (B) Ratio of derived alleles between historical and modern genomes at LoF and missense sites. Error bars, ±1 SD. This ratio is significantly different from 1 in Grauer’s gorillas. (C) The likelihood of an amino acid substitution being damaging for homozygous variants within eastern gorilla populations. Solid lines show pop- ulation averages for all alleles classified as possible and likely damaging (PolyPhen2 score > 0.45). **p < 0.01, ***p < 0.001; NS, not significant. See also Figures S2 and S3. to implemented conservation measures and/or a less pro- ACKNOWLEDGMENTS nounced population decline than previously reported [28, 30]. We thank Emmanuel Gilissen, Patrick Semal, Daniela C. Kalthoff, Gontran So- Many species have experienced severe population declines net, Isabel Ordaz Nemeth, and Zoltan T. Nagy for their support during the study in the last centuries [2] and may thus face similar genomic con- of natural history collections and Matthias Meyer for helpful discussions. sequences as reported here for Grauer’s gorillas. In species Sequencing was performed by the Swedish National Genomics Infrastructure with shorter generation time, such consequences may be (NGI) at the Science for Life Laboratory, which is supported by the Swedish even more severe due to faster genomic change. Nevertheless, Research Council and the Knut and Alice Wallenberg Foundation. The authors if population declines occur more gradually or bottlenecks are acknowledge support from the Uppsala Multidisciplinary Centre for Advanced Computational Science for assistance with massively parallel sequencing and short-lived and followed by population recovery within only a access to the UPPMAX computational infrastructure. This research was sup- few generations, species are likely to be less impacted by ported by BFU2017-86471-P (MINECO/FEDER, UE), U01 MH106874 grant, negative genomic consequences, as exemplified by the moun- Howard Hughes International Early Career, Obra Social ‘‘La Caixa’’ and Secre- tain gorillas. taria d’Universitats i Recerca and CERCA Programme del Departament d’Economia i Coneixement de la Generalitat de Catalunya (GRC 2017 SGR STAR+METHODS 880) to T.M.-B., the Human Frontier Science Program Postdoctoral Fellowship (LT000800/2011-L), the Jan Löfqvist and the Nilsson-Ehle Endowments of the Royal Physiographic Society in Lund, FORMAS (2016-00835) to K.G., and Detailed methods are provided in the online version of this paper FORMAS (2015-676) to L.D. and include the following: AUTHOR CONTRIBUTIONS d KEY RESOURCES TABLE d CONTACT FOR REAGENT AND RESOURCE SHARING Conceptualization, T.v.d.V., L.D., and K.G.; Methodology, T.v.d.V., D.D.-d.-M., d EXPERIMENTAL MODEL AND SUBJECT DETAILS T.M.-B., K.G., and L.D.; Investigation, T.v.d.V., D.D.-d.-M.; Writing, T.v.d.V. d METHOD DETAILS and K.G. with input from all co-authors; Funding Acquisition, L.D. and K.G.; Re- B DNA extraction, library preparation and sequencing sources, L.D. and K.G.; Supervision, L.D. and K.G. B Read mapping, post-mortem DNA damage, SNP call- DECLARATION OF INTERESTS ing and filtering B Post mortem DNA damage in historical genomes The authors declare no competing interests. B Genotype likelihood and SNP calling B PCA and Admixture Received: September 12, 2018 B Heterozygosity and nucleotide diversity Revised: November 12, 2018 B Runs of homozygosity Accepted: November 21, 2018 Published: December 27, 2018 B Genetic load d DATA AND SOFTWARE AVAILABILITY REFERENCES SUPPLEMENTAL INFORMATION 1. Waters, C.N., Zalasiewicz, J., Summerhayes, C., Barnosky, A.D., Poirier, C., Ga1uszka, A., Cearreta, A., Edgeworth, M., Ellis, E.C., Ellis, M., et al. Supplemental Information includes three figures and three data files and can (2016). The Anthropocene is functionally and stratigraphically distinct be found with this article online at https://doi.org/10.1016/j.cub.2018.11.055. from the Holocene. Science 351, aad2622. 168 Current Biology 29, 165–170, January 7, 2019 2. Pimm, S.L., Jenkins, C.N., Abell, R., Brooks, T.M., Gittleman, J.L., Joppa, 22. Tocheri, M.W., Dommain, R., McFarlin, S.C., Burnett, S.E., Troy Case, D., L.N., Raven, P.H., Roberts, C.M., and Sexton, J.O. (2014). The biodiversity Orr, C.M., Roach, N.T., Villmoare, B., Eriksen, A.B., Kalthoff, D.C., et al. of species and their rates of extinction, distribution, and protection. (2016). The evolutionary origin and population history of the grauer gorilla. Science 344, 1246752. Am. J. Phys. Anthropol. 159 (Suppl 61 ), S4–S18. 3. Lande, R. (1993). Risks of Population Extinction from Demographic and 23. van der Valk, T., Sandoval-Castellanos, E., Caillaud, D., Ngobobo, U., Environmental Stochasticity and Random Catastrophes. Am. Nat. 142, Binyinyi, E., Nishuli, R., Stoinski, T., Gilissen, E., Sonet, G., Semal, P., 911–927. et al. (2018). Significant loss of mitochondrial diversity within the last cen- 4. Reed, D.H., and Frankham, R. (2003). Correlation between fitness and ge- tury due to extinction of peripheral populations in eastern gorillas. Sci. netic diversity. Conserv. Biol. 17, 230–237. Rep. 8, 6551. 5. Lande, R., and Shannon, S. (1996). The Role of Genetic Variation in 24. Gazave, E., Chang, D., Clark, A.G., and Keinan, A. (2013). Population Adaptation and Population Persistence in a Changing Environment. growth inflates the per-individual number of deleterious mutations and re- Evolution 50, 434–437. duces their mean effect. Genetics 195, 969–978. 6. Keller, L.F., and Waller, D.M. (2002). Inbreeding effects in wild populations. 25. Peischl, S., and Excoffier, L. (2015). Expansion load: recessive mutations Trends Ecol. Evol. 17, 230–241. and the role of standing genetic variation. Mol. Ecol. 24, 2084–2094. 7. Lande, R. (1988). Genetics and demography in biological conservation. 26. Plumptre, A.J., Nixon, S., Kujirakwinja, D.K., Vieilledent, G., Critchlow, R., Science 241, 1455–1460. Williamson, E.A., Nishuli, R., Kirkby, A.E., and Hall, J.S. (2016). Catastrophic Decline of World’s Largest Primate: 80% Loss of Grauer’s 8. Smith, K.F., Sax, D.F., and Lafferty, K.D. (2006). Evidence for the role of in- Gorilla (Gorilla beringei graueri) Population Justifies Critically Endangered fectious disease in species extinction and endangerment. Conserv. Biol. Status. PLoS ONE 11, e0162697. 20, 1349–1357. 27. Harcourt, A.H., and Fossey, D. (1981). The Virunga gorillas: decline of an 9. Dı́ez-del-Molino, D., Sánchez-Barreiro, F., Barnes, I., Gilbert, M.T.P., and n, L. (2018). Quantifying Temporal Genomic Erosion in Endangered ‘island’ population. Afr. J. Ecol. 19, 83–97. Dale Species. Trends Ecol. Evol. 33, 176–185. 28. Kalpers, J., Williamson, E.A., Robbins, M.M., McNeilage, A., Nzamurambaho, A., Lola, N., and Mugiri, G. (2003). Gorillas in the cross- 10. IUCN (2017). The IUCN Red List of Threatened Species. Version 2017.3. fire: population dynamics of the Virunga mountain gorillas over the past 11. Prado-Martinez, J., Sudmant, P.H., Kidd, J.M., Li, H., Kelley, J.L., Lorente- three decades. Oryx 37, 326–337. Galdos, B., Veeramah, K.R., Woerner, A.E., O’Connor, T.D., Santpere, G., et al. (2013). Great ape genetic diversity and population history. Nature 29. Gray, M., Roy, J., Vigilant, L., Fawcett, K., Basabose, A., Cranfield, M., 499, 471–475. Uwingeli, P., Mburanumwe, I., Kagoda, E., and Robbins, M.M. (2013). Genetic census reveals increased but uneven growth of a critically endan- 12. Xue, Y., Prado-Martinez, J., Sudmant, P.H., Narasimhan, V., Ayub, Q., gered mountain gorilla population. Biol. Conserv. 158, 230–238. Szpak, M., Frandsen, P., Chen, Y., Yngvadottir, B., Cooper, D.N., et al. (2015). Mountain gorilla genomes reveal the impact of long-term popula- 30. Robbins, M.M., Gray, M., Fawcett, K.A., Nutter, F.B., Uwingeli, P., tion decline and inbreeding. Science 348, 242–245. Mburanumwe, I., Kagoda, E., Basabose, A., Stoinski, T.S., Cranfield, M.R., et al. (2011). Extreme conservation leads to recovery of the 13. Gordon, D., Huddleston, J., Chaisson, M.J.P., Hill, C.M., Kronenberg, Virunga mountain gorillas. PLoS ONE 6, e19788. Z.N., Munson, K.M., Malig, M., Raja, A., Fiddes, I., Hillier, L.W., et al. (2016). Long-read sequence assembly of the gorilla genome. Science 31. Rohland, N., and Reich, D. (2012). Cost-effective, high-throughput DNA 352, aae0344. sequencing libraries for multiplexed target capture. Genome Res. 22, 939–946. 14. Sawyer, S., Krause, J., Guschanski, K., Savolainen, V., and Pa €a €bo, S. (2012). Temporal patterns of nucleotide misincorporations and DNA frag- 32. Meyer, M., and Kircher, M. (2010). Illumina sequencing library preparation mentation in ancient DNA. PLoS ONE 7, e34131. for highly multiplexed target capture and sequencing. Cold Spring Harb. Protoc. https://doi.org/10.1101/pdb.prot5448. 15. Korneliussen, T.S., Albrechtsen, A., and Nielsen, R. (2014). ANGSD: Analysis of Next Generation Sequencing Data. BMC Bioinformatics 33. Schubert, M., Lindgreen, S., and Orlando, L. (2016). AdapterRemoval 15, 356. v2: rapid adapter trimming, identification, and read merging. BMC Res. Notes 9, 88. 16. Kardos, M., Åkesson, M., Fountain, T., Flagstad, Ø., Liberg, O., Olason, P., Sand, H., Wabakken, P., Wikenros, C., and Ellegren, H. (2018). Genomic 34. Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with consequences of intensive inbreeding in an isolated wolf population. Nat Burrows-Wheeler transform. Bioinformatics 25, 1754–1760. Ecol Evol 2, 124–131. 35. McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., 17. Pollard, K.S., Hubisz, M.J., Rosenbloom, K.R., and Siepel, A. (2010). Kernytsky, A., Garimella, K., Altshuler, D., Gabriel, S., Daly, M., and Detection of nonneutral substitution rates on mammalian phylogenies. DePristo, M.A. (2010). The Genome Analysis Toolkit: a MapReduce frame- Genome Res. 20, 110–121. work for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303. 18. Adzhubei, I., Jordan, D.M., and Sunyaev, S.R. (2013). Predicting functional effect of human missense mutations using PolyPhen-2. Curr. Protoc. Hum. 36. Danecek, P., Auton, A., Abecasis, G., Albers, C.A., Banks, E., DePristo, Genet. Chapter 7: Unit 7.20. M.A., Handsaker, R.E., Lunter, G., Marth, G.T., Sherry, S.T., et al.; 1000 Genomes Project Analysis Group (2011). The variant call format and 19. Mudakikwa, A., Cranfield, M.R., Sleeman, J.M., and Eilenberger, U. (2005). VCFtools. Bioinformatics 27, 2156–2158. Clinical medicine, preventive health care and research on mountain go- rillas in the Virunga Volcanoes region. In Mountain Gorillas : Three decades 37. Jónsson, H., Ginolhac, A., Schubert, M., Johnson, P.L.F., and Orlando, L. of research at Karisoke, M.M. Robbins, P. Sicotte, and K.J. Stewart, eds. (2013). mapDamage2.0: fast approximate Bayesian estimates of ancient (Cambridge: Cambridge University Press), pp. 342–360. DNA damage parameters. Bioinformatics 29, 1682–1684. 20. Strindberg, S., Maisels, F., Williamson, E.A., Blake, S., Stokes, E.J., Aba’a, 38. Quinlan, A.R., and Hall, I.M. (2010). BEDTools: a flexible suite of utilities for R., Abitsi, G., Agbor, A., Ambahe, R.D., Bakabana, P.C., et al. (2018). comparing genomic features. Bioinformatics 26, 841–842. Guns, germs, and trees determine density and distribution of gorillas 39. Bolger, A.M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible and chimpanzees in Western Equatorial Africa. Sci. Adv. 4, eaar2964. trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. 21. Roy, J., Arandjelovic, M., Bradley, B.J., Guschanski, K., Stephens, C.R., 40. Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, Bucknell, D., Cirhuza, H., Kusamba, C., Kyungu, J.C., Smith, V., et al. G., Abecasis, G., and Durbin, R.; 1000 Genome Project Data Processing (2014). Recent divergences and size decreases of eastern gorilla popula- Subgroup (2009). The Sequence Alignment/Map format and SAMtools. tions. Biol. Lett. 10, 20140811. Bioinformatics 25, 2078–2079. Current Biology 29, 165–170, January 7, 2019 169 , J., Marco Sola, S., Knowles, D.G., Raineri, E., Guigó, 41. Derrien, T., Estelle 50. Nielsen, R., Paul, J.S., Albrechtsen, A., and Song, Y.S. (2011). Genotype R., and Ribeca, P. (2012). Fast computation and applications of genome and SNP calling from next-generation sequencing data. Nat. Rev. mappability. PLoS ONE 7, e30377. Genet. 12, 443–451. 42. Narasimhan, V., Danecek, P., Scally, A., Xue, Y., Tyler-Smith, C., and Durbin, 51. DePristo, M.A., Banks, E., Poplin, R., Garimella, K.V., Maguire, J.R., Hartl, R. (2016). BCFtools/RoH: a hidden Markov model approach for detecting C., Philippakis, A.A., del Angel, G., Rivas, M.A., Hanna, M., et al. (2011). autozygosity from next-generation sequencing data. Bioinformatics 32, A framework for variation discovery and genotyping using next-generation 1749–1751. DNA sequencing data. Nat. Genet. 43, 491–498. 43. McLaren, W., Gil, L., Hunt, S.E., Riat, H.S., Ritchie, G.R.S., Thormann, A., 52. Li, H. (2014). Toward better understanding of artifacts in variant calling Flicek, P., and Cunningham, F. (2016). The Ensembl Variant Effect from high-coverage samples. Bioinformatics 30, 2843–2851. Predictor. Genome Biol. 17, 122. 53. Fumagalli, M., Vieira, F.G., Korneliussen, T.S., Linderoth, T., Huerta- Sánchez, E., Albrechtsen, A., and Nielsen, R. (2013). Quantifying popula- 44. Eden, E., Navon, R., Steinfeld, I., Lipson, D., and Yakhini, Z. (2009). GOrilla: tion genetic differentiation from next-generation sequencing data. a tool for discovery and visualization of enriched GO terms in ranked gene Genetics 195, 979–992. lists. BMC Bioinformatics 10, 48. 54. Skotte, L., Korneliussen, T.S., and Albrechtsen, A. (2013). Estimating indi- 45. Dabney, J., Knapp, M., Glocke, I., Gansauge, M.-T., Weihmann, A., Nickel, vidual admixture proportions from next generation sequencing data. €a B., Valdiosera, C., Garcı́a, N., Pa €bo, S., Arsuaga, J.-L., and Meyer, M. Genetics 195, 693–702. (2013). Complete mitochondrial genome sequence of a Middle 55. Nielsen, R., Korneliussen, T., Albrechtsen, A., Li, Y., and Wang, J. (2012). Pleistocene cave bear reconstructed from ultrashort DNA fragments. SNP calling, genotype calling, and sample allele frequency estimation Proc. Natl. Acad. Sci. USA 110, 15758–15763. from New-Generation Sequencing data. PLoS ONE 7, e37558. 46. Ersmark, E., Orlando, L., Sandoval-Castellanos, E., Barnes, I., Barnett, R., 56. Thompson, E.A. (2013). Identity by descent: variation in meiosis, across n, L. (2015). Population Demography and Stuart, A., Lister, A., and Dale genomes, and in populations. Genetics 194, 301–326. Genetic Diversity in the Pleistocene Cave Lion. Open Quat. 1, 1–15. , L.R., Bustamante, C.D., Clark, A.G., and Gravel, S. 57. Henn, B.M., Botigue n, L., and Guschanski, K. (2017). 47. van der Valk, T., Lona Durazo, F., Dale (2015). Estimating the mutation load in human genomes. Nat. Rev. Whole mitochondrial genome capture from faecal samples and Genet. 16, 333–343. museum-preserved specimens. Mol. Ecol. Resour. 17, e111–e121. 58. Cooper, G.M., and Shendure, J. (2011). Needles in stacks of needles: 48. Li, H. (2013). Aligning sequence reads, clone sequences and assembly finding disease-causal variants in a wealth of genomic data. Nat. Rev. contigs with BWA-MEM. arXiv, arXiv:1303.3997, https://arxiv.org/abs/ Genet. 12, 628–640. 1303.3997. 59. Librado, P., Gamba, C., Gaunitz, C., Der Sarkissian, C., Pruvost, M., €a 49. Briggs, A.W., Stenzel, U., Meyer, M., Krause, J., Kircher, M., and Pa €bo, S. Albrechtsen, A., Fages, A., Khan, N., Schubert, M., Jagannathan, V., (2010). Removal of deaminated cytosines and detection of in vivo methyl- et al. (2017). Ancient genomic changes associated with domestication of ation in ancient DNA. Nucleic Acids Res. 38, e87. the horse. Science 356, 442–445. 170 Current Biology 29, 165–170, January 7, 2019 STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Chemicals, Peptides, and Recombinant Proteins Tango Buffer Thermofisher Scientific Cat#: BY5 dNTP’s Thermofisher Scientific Cat#: R1121 ATP Thermofisher Scientific Cat#: R0441 T4 PNK Thermofisher Scientific Cat#: EK0031 USER enzyme NEB Cat#: M5505S T4 DNA ligase Buffer Thermofisher Scientific Cat#: EL0014 T4 DNA Polymerase Thermofisher Scientific Cat#: EP0061 PEG-4000 Merck Cat#: 25322-68-3 T4 DNA ligase Thermofisher Scientific Cat#: EL0014 Isothermal Amplification buffer New England Biolabs Cat#: B0537S Bst polymerase 2.0, Large Fragment New England Biolabs Cat#: M0275S Pfu Turbo Cx Hotstart DNA Buffer Agilent technologies Cat#: 600410 Pfu Turbo Cx Hotstart DNA Polymerase Agilent technologies Cat#: 600410 TritonX Merck Cat#: 600410 Proteinase K Sigma Cat#: 9002-93-1 EDTA Thermofisher Scientific Cat#: 15575020 Guanidine hydrochloride Sigma Cat#: G3272-500G Tween-20 Sigma Cat#: P9416-50ML Isopropanol Sigma Cat#: 650447 Tris-HCl Thermofisher Scientific Cat#: AM9855G Critical Commercial Assays High Pure Viral Nucleic Acid Large Volume Kit Roche Cat#: 05114403001 MinElute PCR Purification Kit QIAGEN Cat#: 28004 AMPure XP Beads Agencourt Cat#: A63881 Deposited Data Raw and analyzed data This paper ENA: PRJEB29503 Gorilla reference genome [13] ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/ 151/905/GCF_000151905.2_gorGor4 Gorilla whole genome resequencing data [11] https://www.ncbi.nlm.nih.gov/bioproject/189439 Gorilla whole genome resequencing data [12] https://www.ebi.ac.uk/ena/data/view/PRJEB3220 Base-wise phyloP scores [17] ftp://hgdownload.cse.ucsc.edu/goldenPath/mm9/ phyloP30way/placentalMammals/ hg38ToGorGor4 chainfile N/A http://hgdownload.cse.ucsc.edu/gbdb/hg38/liftOver/ Experiment models: Organism 38 Gorilla beringei graueri museum specimens Royal Museum for Central Africa See Data S1 (teeth and skin) (Tervuren, Belgium) and Royal Belgian Institute of Natural Sciences (Brussels, Belgium) 21 Gorilla beringei beringei museum specimens Swedish Museum of Natural See Data S1 (teeth) History (Stockholm, Sweden) Oligonucleotides Barcode-primer1: CytoGene [31] 50 -> 30 CTTTCCCTACACGACGCTCTTCCGA TCTatcgatt (Continued on next page) Current Biology 29, 165–170.e1–e6, January 7, 2019 e1 Continued REAGENT or RESOURCE SOURCE IDENTIFIER Barcode-primer2: CytoGene [31] 50 -> 30 CTTTCCCTACACGACGCTCTTCCG ATCTaacgcta Barcode-primer3: CytoGene [31] 50 -> 30 CTTTCCCTACACGACGCTCTTCCG ATCTggtatcg Barcode-primer4: CytoGene [31] 50 -> 30 CTTTCCCTACACGACGCTCTTCCG ATCTtgactgg Barcode-primer5: CytoGene [31] 50 -> 30 CTTTCCCTACACGACGCTCTTCCG ATCTcaattgc Barcode-primer6: CytoGene [31] 50 -> 30 CTTTCCCTACACGACGCTCTTCCG ATCTtggccat IS5_reamp.P5: CytoGene [32] 50 -> 30 AATGATACGGCGACCACCGA IS6_reamp.P7: CytoGene [32] 50 -> 30 CAAGCAGAAGACGGCATACGA Pre_Hyb_F: CytoGene [31] 50 -> 30 CTTTCCCTACACGACGCTCTTC Pre_Hyb_R: CytoGene [31] 50 -> 30 GTGACTGGAGTTCAGACGTGTGCT P7 index and sequencing-primer: CytoGene [32] 50 -> 30 AAGCAGAAGACGGCATACGAGAT gaagattGTGACTGGAGTTCAGACGTGT P5 sequencing-primer: CytoGene [32] 50 -> 30 AAGCAGAAGACGGCATACGAGATG TGACTGGAGTTCAGACGTGT Software and Algorithms AdapterRemovalV2 [33] https://github.com/MikkelSchubert/adapterremoval/; RRID: SCR_011834 BWA [34] http://bio-bwa.sourceforge.net/; RRID: SCR_010910 GATK v3.8 [35] https://software.broadinstitute.org/gatk/; RRID: SCR_001876 Picard N/A https://broadinstitute.github.io/picard/; RRID: SCR_006525 Vcftools [36] https://vcftools.github.io/index.html; RRID: SCR_001235 MapDamage2.0 [37] ginolhac.github.io/mapDamage/; RRID: SCR_001240 BEDTools [38] https://bedtools.readthedocs.io/en/latest/; RRID: SCR_006646 Trimmomatic [39] http://www.usadellab.org/cms/?page=trimmomatic; RRID: SCR_011848 SAMtools [40] http://samtools.sourceforge.net/; RRID: SCR_002105 Python 2.7 N/A https://www.python.org/downloads/release/python- 2715/; RRID: SCR_008394 R N/A https://www.r-project.org/; RRID: SCR_001905 ANGSD [15] https://github.com/ANGSD/angsd GEM [41] http://gemlibrary.sourceforge.net htsbox N/A https://github.com/lh3/htsbox (Continued on next page) e2 Current Biology 29, 165–170.e1–e6, January 7, 2019 Continued REAGENT or RESOURCE SOURCE IDENTIFIER bcftools [42] http://www.htslib.org/doc/bcftools.html; RRID: SCR_005227 liftOver N/A https://genome.sph.umich.edu/wiki/LiftOver; RRID: SCR_001173 Ensembl Variant Effect Predictor [43] https://www.ensembl.org/info/docs/tools/vep/ index.html; RRID: SCR_007931 PolyPhen-2 [18] http://genetics.bwh.harvard.edu/pph2/; RRID: SCR_013189 GOrilla [44] http://cbl-gorilla.cs.technion.ac.il/; RRID: SCR_006848 CONTACT FOR REAGENT AND RESOURCE SHARING Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Katerina Guschanski (katerina.guschanski@ebc.uu.se). EXPERIMENTAL MODEL AND SUBJECT DETAILS Historical Grauer’s gorilla samples (n = 38) with age ranging from 108 – 32 years (median year of collection = 1949) were collected at the Royal Belgian Institute of Natural Sciences (RBINS) in Brussels and the Royal Museum for Central Africa (RMCA) in Tervuren, Belgium. Historical mountain gorilla samples from the Virunga Massif population (n = 21) with age ranging from 104 – 61 years (median year of collection = 1921) were obtained from the Swedish Museum of Natural History (NRM) in Stockholm. METHOD DETAILS DNA extraction, library preparation and sequencing DNA from historical samples, consisting of teeth and dried soft-tissue samples, was extracted in dedicated ancient DNA facilities at Uppsala University, Uppsala, and the Swedish Museum of Natural History, Stockholm. Surface contamination from teeth was removed by exposing the roots to UV light (245 nm) for 10 min, abrading the surface with an engraving cutter using a Dremel 8100 drill, and again subjecting the cleaned surfaces to UV light for 10 min. No surface decontamination was performed for soft-tissue samples. The Grauer’s gorilla samples were extracted using an extraction protocol optimized for retention of short fragments in ancient DNA samples [45]. Mountain gorilla samples were extracted with Vivapsin filters as in Ersmark et. al. [46]. For screening, we prepared pooled Illumina sequencing libraries from all mountain and Grauer’s gorilla DNA extracts following the strategy outlined in [31] and [47]. Briefly, 20 mL of DNA extract was used in a 40 mL blunting reaction (final concentrations: 1 3 buffer Tango, 100 mM each dNTP, 1 mM ATP, 25 U T4 polynucleotide kinase (Thermo Scientific), 5U T4 DNA polymerase (Thermo Scientific)). Samples were incubated at 25 C for 15 min followed by 5 min at 12 C. DNA fragments within each sample were then ligated to a combination of incomplete, partially double-stranded P5- and P7-adapters (10 mM each), each containing a unique seven base pair sequence, referred to as barcode. Adapter ligation was performed in a 40 mL reaction volume using 20 mL of blunted DNA and 1 mL of unique P5 and P7 barcoded adapters per sample (final concentrations: 1 3 T4 DNA ligase buffer, 5% PEG-4000, 5U T4 DNA ligase (Thermo Scientific)). Samples were incubated for 30 min at room temperature and cleaned using MinElute spin columns following the man- ufacturer’s protocol. Adapter fill-in was performed in 40 mL final volume using 20 mL adapter ligated DNA (final concentrations: 1 3 T4 DNA ligase buffer, 5% PEG-4000, 5 U T4 DNA ligase (Thermo Scientific)), incubated at 37 C for 20 min, heat-inactivated at 80 C for 20 min, and cleaned using MinElute spin columns, as above. For the purpose of shallow shotgun sequencing, indexing PCR was performed for 10 cycles in 25 mL reaction volume (3 mL of adapter-ligated library as input) using a unique P7 indexing primer for each sample, as in Meyer & Kircher (2010) [32] (final concentrations: 1x AccuPrime reaction mix, 0.3 mM IS4 primer, 0.3 mM P7 index- ing primer, 7 U AccuPrime Pfx (Thermo Scientific), cycling protocol: 95 C for 2 min, 10 cycles at 95 C for 30 s, 55 C for 30 s and 72 C for 1 min and a final extension at 72 C for 5 min). Sample libraries were cleaned using MinElute spin columns and subsequently frag- ment length distribution and concentration were measured on the Bioanalyzer. Samples were pooled in equimolar amounts into two final library pools and size selected with two rounds of AMPure XP bead clean-up, using 0.5X and 1.8X bead:DNA ratio, respectively. Both library pools with final concentration of 10nM were subjected to 125 bp paired-end shotgun sequencing on the Illumina HiSeq 2500 platform (one lane each, High Output Mode) at the SciLifeLab sequencing facility in Stockholm. Samples were demultiplexed based on their unique indices using bcl2fastq v2.17.1 with default settings (Illumina Inc.). We removed adapter sequences, including barcodes, and used AdapterRemovalV2 [33] with default parameters to collapse and remove reads shorter than 15 base pairs. Collapsed reads were mapped against the gorilla reference genome (gorGor4) [13] using bwa-mem with default parameters [48]. We then estimated the endogenous content for each sample as the fraction of reads mapping to the reference genome with a quality Current Biology 29, 165–170.e1–e6, January 7, 2019 e3 score > 20 (Phred-scale) (see Data S1). Note that we did not remove duplicate reads at this stage, since duplicate removal tools (sam- tools rmdup and Picard MarkDuplicates) only remove duplicates among the mapped reads, biasing the endogenous content esti- mates. Post-mortem DNA damage patterns were identified using mapDamage2.0 [37]. We also assessed human contamination as in van der Valk et al. (2017) [47]. Briefly, we identified fixed nucleotide difference between a panel of mitochondrial genomes from nine human with diverse geographic origin [11] and 13 published eastern gorilla mitochondrial genomes [12]. For each of these diagnostic sites covered by at least three sequencing reads, we calculated the proportion of reads containing the human variant and averaged the estimate over all diagnostic sites (see Data S1). After screening we selected seven Grauer’s and four mountain gorilla samples displaying high endogenous content and low levels of human contamination for UDG treatment to remove post-mortem DNA damage and subsequent deep sequencing. Libraries were then prepared as above with slight modifications: During blunt end repair, 20 mL of DNA extract was used in a 50 mL reaction together with USER enzyme treatment to remove uracil bases resulting from post-mortem damage (final concentrations: 1 3 buffer Tango, 100 mM each dNTP, 1 mM ATP, 25 U T4 polynucleotide kinase (Thermo Scientific), 3U USER enzyme (NEB)). Samples were incubated for 3 h at 37 C, followed by the addition of 1 mL T4 DNA polymerase (Thermo Scientific) and incubation at 25 C for 15 min and 12 C for 5 min. Indexing PCR for all adapter ligated libraries was performed as above for 10 cycles, but in 125 mL volume (18 mL of adapter ligated input DNA) and each of these samples was sequenced on one HiSeq X lane 150bp paired-end. Samples with insufficient coverage (n = 7) after deep-sequencing were pooled and sequenced on an additional five HiSeq X lanes. Read mapping, post-mortem DNA damage, SNP calling and filtering Pooled libraries were demultiplexed based on their unique 7bp barcodes (see Data S1) and we subsequently removed adapter and barcode sequences and collapsed fastq-reads by individual using AdapterRemovalV2 [33]. Collapsed reads were mapped against the western lowland gorilla reference genome with bwa-mem on default settings (including readgroups and marking shorter split hits as secondary for Picard compatibility) [48]. After mapping, alignments around indels were improved using GATKv3.8 indelRealigner [35]. We then removed duplicate reads using Picard (https://broadinstitute.github.io/picard/). We estimated remaining levels of post- mortem DNA damage (note that the libraries were UDG-treated) and rescaled mapping quality scores to correct for any remaining damage using Mapdamage2.0 on default parameters [37]. After quality score rescaling, we filtered out all reads below mapping qual- ity of 20 (Phred-scale) and removed sex chromosomes from further analyses. Raw fastq files from published genomes of wild-born gorillas (17 western lowland, 7 mountain, and 8 Grauer’s gorillas) were obtained from the European nucleotide archive (PRJNA189439 and PRJEB3220) [11, 12]. We removed sequencing adapters and subsequently bases at the start and end of the reads were trimmed if falling below the quality score of 15 (Phred33-scale) using Trim- momatic [39]. Read mapping, indel-realignment and duplicate removal was performed as above. Post mortem DNA damage in historical genomes All historical libraries were UDG-treated, removing the majority of post-mortem DNA damage [49]. We confirmed (using mapDamage2.0) that remaining typical ancient DNA post-mortem damage (C-to-T substitutions at read ends) in historical samples was present only at very low rates (Figure S1A). However UDG-treatment has lower efficiency of removing uracil nucleotides at CpG sites [49] and therefore remaining post-mortem DNA damage at these positions may cause biases in downstream analyses. We estimated the extent of damage at CpG sites as the fraction of mapped reads containing a thymine nucleotide at a cytosine position that is part of a CpG site. To this end, we first identified all CpG sites in the gorilla reference genome defined as sites where C is imme- diately followed by a G and used samtools mpileup to obtain all mapped bases with Q > 20 at each site per genome. We then quan- tified the fraction of thymine nucleotides at the identified CpG sites for each genome and averaged this fraction over all identified CpG sites within the genome (Figure S1B). Genotype likelihood and SNP calling Most of our analyses were conducted in a genotype likelihood framework implemented in the ANGSD software, allowing for the incor- poration of statistical uncertainty in low-coverage data [15, 50]. Genotype likelihoods for all gorilla samples were calculated adjusting mapping quality for regions with excessive mismatches (-C 50), removing reads with multiple hits (-uniqueOnly 1), base quality score below 20 (-minQ 20), and downgrading quality scores around indels (-baq 1). We then only considered genotypes with a likelihood ratio test statistic of minimum 24 (-SNP_pval 2e-6) using the samtools genotype model (-GL 1). SNPs were identified by creating a genomic VCF file for each individual using GATK UnifiedGenotyper with default parameters [35]. Subsequently, we used the GATK hard-filtering best-practice guideline for SNP-filtering [51]. Vcftools [36] was used to remove all sites with coverage depth below 3x, minimal allele count below 3, sites for which more than 75% of samples had a non-reference allele in a heterozygous state, indels, sites with higher than 3-fold average coverage across all individuals, and positions with SNP quality < 30 (Phred-scale), as recommended in [52]. Additionally, we removed all SNP within repetitive regions as identified with the GEM module [41]. Finally, we removed all SNPs at CpG sites, since these showed patterns of remaining post-mortem DNA dam- age (Figure S1B). We extensively assessed the effects of differential coverage between historical and modern genomes on popula- tion-level summary statistics by downsampling all modern genomes to similar coverage as the historical genomes (4X) using sam- tools V1.5 [40] and obtaining genotype-likelihoods for the downsampled genomes as described above. We observed slight biases on diversity estimates related to differential coverage (Figure S1C), and therefore all genomic comparisons were subsequently based on the modern genomes downsampled to 4X coverage. e4 Current Biology 29, 165–170.e1–e6, January 7, 2019 PCA and Admixture PCA plots were constructed using PCAngsd with default parameters and 200 EM iterations for computing the population allele fre- quencies [53]. To further explore genetic affinities and substructure in our data we used NGSadmix to estimate individual genome admixture proportions [54]. NGSadmix was run with the filtered genotype likelihoods and additional filtering of alleles below a fre- quency of 0.05 (-minMaf 0.05) (Figures S2A–S2E). Heterozygosity and nucleotide diversity We measured genome-wide heterozygosity using two different approaches. First, we estimated heterozygosity from the site fre- quency spectrum, obtained from the genotype-likelihoods with realSFS as implemented in ANGSD. We excluded transitions to eliminate possible biases from post-mortem DNA damage, such as those at CpG sites (Figures S1A and S1B), and performed 1000 bootstrap replicates for each sample to obtain confidence intervals (Figure S1C) [15, 55]. Second, we estimated genome- wide heterozygosity directly from the genomic VCF file (with down-sampled modern genomes and excluding CpG sites) by counting the proportion of heterozygous sites out of total sites covered by at least five reads in the filtered genomic VCF data (Figure S1D). In addition to estimating heterozygosity for each individual, we obtained population-level estimates of genetic diversity by calculating the number of divergent sites between individuals belonging to the same temporal population (treating historical and modern ge- nomes of different gorilla taxa separately). To this end, we randomly sampled a single allele at each covered site excluding low quality bases (htsbox v1.0 -R -q 30 -Q 30 -l 35 -s 1, https://github.com/lh3/htsbox) within each individual genome. We then calculated the number of within population divergent sites in non-overlapping sliding windows of 10.000 base pairs, subsampling four random ge- nomes for each window to eliminate the possibility that results are driven by differences in sample size between populations (Figure S1E). We performed a permutation test to control for the possibility that the inferred differences in levels of genome-wide heterozygosity between modern and historical eastern gorillas could be driven by differences in sample size. For each dataset (i.e., VCF files of his- torical and modern genomes), we repeatedly (n = 1000) randomly subsampled four genomes (the smallest sample size in our analysis, equal to the number of historical mountain gorilla genomes), calculated genome-wide heterozygosity as described above, and assessed the significance of population level differences using pairwise t-tests in Scipy python package (v1.13.3). In 96.6% of randomization tests, the historical Grauer’s gorilla populations shows significantly higher (p < 0.05) population level heterozygosity (Figure S2F), whereas none of the repeated 1000 comparisons were significant in mountain gorillas (p value > 0.05). Population struc- ture in Grauer’s gorillas (Figures 1B, S2A, and S2B) [12] could potentially explain the observed differences in heterozygosity between modern and historical samples. Although population structure was found in both modern and historical Grauer’s gorillas, we per- formed an additional test to verify that it does not explain the observed decline in heterozygosity in this species. To this end, we esti- mated the average genome-wide heterozygosity based on the genomic VCF files within the larger Grauer’s gorilla sub-cluster (popB) (Figure S2D), as it provided a large enough sample size for this analysis. Confirming our results at population level, heterozygosity was significantly higher (p < 0.05) in historical individuals in the sub-cluster popB compared to modern individuals from the same sub- cluster (Figure S2G). We then used the same approach as above to randomly subsample three historical and three modern individuals from popB (note three individuals is the lowest number required for significance testing) and found none of the tests resulted in higher heterozygosity in the modern dataset (Figure S2H). Runs of homozygosity To identify regions of the genome in complete homozygosity (runs of homozygosity, ROH), which are indicative of inbreeding, we applied a Hidden Markov Model (HMM) first described by Xue et al. (2015) [12] to the filtered SNP dataset assuming a constant recombination rate of 1 3 108 per base pair (Figure S3A). The emission probabilities in the HMM correspond to a Hardy-Weinberg model with inbreeding. Transition probabilities in the HMM incorporate the likelihood of a recombination event since the last site (see [42], for model details). The HMM model employed to infer ROH assigns confidence score to each of the genomic regions inferred to be in complete autozygosity. Even at relatively short ROH length (100.000bp), the majority of such regions have a high confidence score (Figure S3B). This method thus allows for identification of genome stretches in complete homozygosity despite medium to low coverage. However, applying the model to low coverage genomes leads to a slight over-estimation of the number and length of homozygous tracts (Figure S3C), particularly for long ROH (> 2.5Mbp) (Figure S3D). Therefore, we controlled for differential coverage by down-sampling the modern genomes to 4x coverage as described above, before carrying out the comparisons of the distribution of ROH between historical and modern data (Figures 2B and 2C). Runs of homozygosity arise when identical haplotypes are inherited from both parents. Over time these runs are broken apart by recombination events. Therefore, long runs of homozgyosity represent recent inbreeding events, whereas shorter runs stem from more ancient processes or mating between distant relatives. The largest difference between historical and modern genomes stems from ROH that are 5-6 Mb in length (Figure 2C). We estimate the mean number of generations back to the common ancestor of these homologous sequences using the formula from [56]; g = 100=2$ROHlength where g is the number of generations (exponentially distributed) and ROHlength the ROH length in centimorgans. We use the physical length of ROH as approximation for the genetic length and estimate that ROH of 5-6 Mb trace back to less < 10 generations ago. Current Biology 29, 165–170.e1–e6, January 7, 2019 e5 Genetic load Estimating an individual’s genetic load based on genomic data is challenging [57]. However, using evolutionary constraints as pre- dictor for the fitness consequences of a mutation can serve as an approximation [58]. We therefore measured genetic load in each individual as the number of homozygous derived alleles in sites that are under strict evolutionary constraints divided by the total num- ber of homozygous alleles, following the approach by Librado et al. (2017) [59] (Figure 3A). Base-wise phyloP scores, pre-calculated for a whole-genome alignment of 20 mammals were obtained from Ensembl (Release 90). These scores show high accuracy in measuring regions of the genome under evolutionary constraints [17]. We converted the human genomic coordinates of the phyloP scores to the gorilla reference genome using the command-line version of the liftOver tool (https://genome-store.ucsc.edu/) and the hg38ToGorGor4 chainfile available from UCSC (http://hgdownload.cse.ucsc.edu/ goldenPath/hg38/liftOver/). All sites with phyloP scores < 1, indicating low mammalian sequence conservation, were excluded. We considered only homozygous sites, since estimating the genetic load at heterozygous positions requires additional assumptions about the dominance coefficient of mutations. Using the genotype likelihood data, we then defined the genetic load (phyloPi) as the phyloP score at each given position multiplied by the probability of observing a derived homozygous allele (Pi) divided by the total number of homozygous sites in a given individual (hom). To identify loss-of-function (LoF) and missense mutations in each gorilla individual (Figure 3B), we used the Ensembl Variant Effect Predictor (VEP) tool [43], which predicts the consequences of genomic variants on the protein sequence. Annotation of genes and transcripts for the gorilla reference genome (gorGor4) assembly were retrieved from the Ensembl VEP database (release 91, ftp://ftp. ensembl.org/pub/current_variation/VEP/). We used VEP to obtain classification for each of the SNPs in the filtered SNP dataset and subsequently classified SNPs into LoF, missense, and synonymous variants as in Xue et al. (2015): d LoF variants: transcript_ablation, splice_donor_variant, splice_acceptor_variant, stop_gained, frameshift_variant, inframe_in- sertion, inframe_deletion, splice_region_variant d Missense variants: missense_variant d Synonymous variants: synonymous_variant We compared the frequency of LoF variants between historical and modern data following the strategy outlined in Xue et al. (2015). Briefly, for each category of variants, we calculated at each site i the observed allele frequency in Population A as fAi = dAi / nAi, where nAi, is the total number of alleles called in population A and dAi is the total number of called derived alleles. Similarly, we define fBi for population B. Then, if C is a category of protein-coding sites: X   LA;B ðCÞ = fiA 1  fiB i˛C We then calculate the ratio RA/B(C) = LA,B(C) / LB,A(C), which is a measure of the relative number of derived alleles found more often in population A compared to population B. We estimated the variance in RA/B(C) by running a Jackknife approach in blocks of 1000 from the set of sites in C. We compared the frequency of variants in historical and modern data as measured by RA/B(C) in two categories: 1) missense changes in the protein coding sequence, and 2) variants which are likely to cause loss of function. To predict the likelihood of a variant being deleterious we used PolyPhen-2 (Figure 3C) [18]. This tool predicts the effects of amino acid substitutions on the structure and function of human proteins using a set of comparative and structural criteria. SNPs identified in the gorilla genome were lifted over to the human reference (hg38) using the command-line version of the liftOver tool (https:// genome-store.ucsc.edu/) and the gorGor4ToHg38 chainfile available from UCSC (http://hgdownload.cse.ucsc.edu/goldenPath/ gorGor4/liftOver/). Next, we used the PolyPhen-2 implementation of the Variant Effect Predictor tool to obtain PolyPhen-2 scores, with higher scores indicating a higher likelihood of an allele being damaging (scores above 0.45 indicate possibly damaging muta- tions). Finally, to obtain an estimate of the number of putatively damaging alleles in the historical and modern populations, we randomly subsampled an allele for each site among all individual from the respective population and retrieved the associated PolyPhen-2 score for each of these alleles. For all identified LoF variants (see above) we obtained the historical and modern allele frequencies from the VCF-files and identified all LoF alleles with an absolute allele frequency difference between the modern and historical populations above 0.2 (see Data S2). For each of these LoF alleles we obtained the associated gene based on the human genome annotation by lifting over the variant to the human genome using the command-line version of the liftOver tool (https://genome-store.ucsc.edu/) and the gorGor4ToHg38 chainfile available from UCSC. We then assessed the functional enrichment of LoF variants with the highest difference in allele fre- quency between modern and historical populations in GOrilla [44] (see Data S3). DATA AND SOFTWARE AVAILABILITY The accession number for the FASTQ files reported in this paper is ENA: PRJEB29503. e6 Current Biology 29, 165–170.e1–e6, January 7, 2019