Article Ancient and modern genomes unravel the evolutionary history of the rhinoceros family Graphical abstract Authors Shanlin Liu, Michael V. Westbury, Nicolas Dussex, ..., Pierre-Olivier Antoine, Love Dalén, M. Thomas P. Gilbert Correspondence shanlin.liu@cau.edu.cn (S.L.), love.dalen@nrm.se (L.D.), tgilbert@sund.ku.dk (M.T.P.G.) In brief The comparison of de novo genomes from the white, black, Sumatran, and greater one-horned rhinoceroses with the genomes of a historic Javan rhinoceros and three extinct Pleistocene species resolves the evolutionary relationships within the Rhinocerotidae family and reveals that low genetic diversity is a long-term feature of rhinoceroses. Highlights d Analysis of genomes from all five extant and three extinct rhinoceros species d Strong phylogenomic support for the geographical hypothesis of rhinoceros evolution d Basal split between African and Eurasian lineages in the early Miocene (16 mya) d While all rhinoceroses have low genome diversity, it is lowest in modern-day ones Liu et al., 2021, Cell 184, 1–12 September 16, 2021 ª 2021 The Authors. Published by Elsevier Inc. https://doi.org/10.1016/j.cell.2021.07.032 ll Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article Ancient and modern genomes unravel the evolutionary history of the rhinoceros family Shanlin Liu,1,2,37,* Michael V. Westbury,2 Nicolas Dussex,3,4,5 Kieren J. Mitchell,6 Mikkel-Holger S. Sinding,2 Peter D. Heintzman,7 David A. Duchêne,2 Joshua D. Kapp,8 Johanna von Seth,3,4,5 Holly Heiniger,6 Fátima Sánchez-Barreiro,2 Ashot Margaryan,2 Remi André-Olsen,9 Binia De Cahsan,2 (Author list continued on next page) 1Department of Entomology, College of Plant Protection, China Agricultural University, Beijing 100193, China 2The GLOBE Institute, Faculty of Health and Medical Sciences, University of Copenhagen, 1353 Copenhagen, Denmark 3Centre for Palaeogenetics, Svante Arrhenius vag 20C, Stockholm 10691, Sweden 4Department of Bioinformatics and Genetics, Swedish Museum of Natural History, Stockholm 10405, Sweden 5Department of Zoology, Stockholm University, Stockholm 10691, Sweden 6Australian Centre for Ancient DNA, School of Biological Sciences, University of Adelaide, Adelaide 5005, Australia 7The Arctic University Museum of Norway, UiT The Arctic University of Norway, Tromsø 9037, Norway 8Department of Ecology and Evolutionary Biology, University of California, Santa Cruz, Santa Cruz, CA 95064, USA 9Science for Life Laboratory, Department of Biochemistry and Biophysics, Stockholm University, 17121 Solna, Sweden 10China National Genebank, BGI Shenzhen, Shenzhen 518083, China 11Center for Ecological and Environmental Sciences, Northwestern Polytechnical University, Xi’an 710072, China 12Department of Cell and Molecular Biology, National Bioinformatics Infrastructure Sweden, Science for Life Laboratory, Uppsala University, Uppsala, Sweden 13Department of Zoology, University of Venda, Thohoyandou 0950, Republic of South Africa 14Editor of the Rhino Resource Center, Utrecht, the Netherlands 15School of Biosciences, Sir Martin Evans Building, Cardiff University, Cardiff CF10 3AX, UK 16Sustainable Places Research Institute, Cardiff University, Cardiff CF10 3BA, UK 17San Diego Zoo Wildlife Alliance, Beckman Center for Conservation Research, San Diego, CA 92027, USA 18Rotterdam Zoo, Rotterdam, the Netherlands 19N.A. Shilo North-East Interdisciplinary Scientific Research Institute, Far East Branch, Russian Academy of Sciences (NEISRI FEB RAS), Magadan 685000, Russia (Affiliations continued on next page) SUMMARY Only five species of the once-diverse Rhinocerotidae remain, making the reconstruction of their evolutionary history a challenge to biologists since Darwin. We sequenced genomes from five rhinoceros species (three extinct and two living), which we compared to existing data from the remaining three living species and a range of outgroups. We identify an early divergence between extant African and Eurasian lineages, resolving a key debate regarding the phylogeny of extant rhinoceroses. This early Miocene (16 million years ago [mya]) split post-dates the land bridge formation between the Afro-Arabian and Eurasian landmasses. Our analyses also show that while rhinoceros genomes in general exhibit low levels of genome-wide diversity, heterozygosity is lowest and inbreeding is highest in the modern species. These results suggest that while low genetic diversity is a long-term feature of the family, it has been particularly exacerbated recently, likely reflecting recent anthropogenic-driven population declines. INTRODUCTION all of which are highly endangered and global priorities for con- servation. Rhinocerotoidea, the clade including the rhinoceros Understanding the relationships among rhinoceros species and family (Rhinocerotidae), diverged from tapirs 55–60 million years when they diverged has been a question addressed by evolu- ago (mya) in either Eurasia or North America (Bai et al., 2020). The tionary biologists since the dawn of the field. Darwin himself dis- family subsequently radiated into at least 100 species distributed cussed the topic in 1842 as one of a handful of examples in his across Africa, Eurasia, North, and Central America (Cerdeño, short treatise on evolution that preceded On the Origin of Spe- 1998) and included some of the largest land mammals that cies in 1859 (Darwin, 1909). Although rhinoceroses were once ever lived. Most rhinocerotids went extinct prior to the Pleisto- a diverse clade, extant rhinoceroses comprise only five species, cene, with just nine species surviving into the Late Pleistocene, Cell 184, 1–12, September 16, 2021 ª 2021 The Authors. Published by Elsevier Inc. 1 This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article Guanliang Meng,10 Chentao Yang,10 Lei Chen,11 Tom van der Valk,12 Yoshan Moodley,13 Kees Rookmaaker,14 Michael W. Bruford,15,16 Oliver Ryder,17 Cynthia Steiner,17 Linda G.R. Bruins-van Sonsbeek,18 Sergey Vartanyan,19 Chunxue Guo,10 Alan Cooper,20 Pavel Kosintsev,21,22 Irina Kirillova,23 Adrian M. Lister,24 Tomas Marques-Bonet,25,26,27,28 Shyam Gopalakrishnan,2 Robert R. Dunn,2,29 Eline D. Lorenzen,2 Beth Shapiro,8,30 Guojie Zhang,10,31,32,33 Pierre-Olivier Antoine,34 Love Dalén,3,4,5,36,* and M. Thomas P. Gilbert2,35,36,* 20South Australian Museum, Adelaide, SA 5000, Australia 21Institute of Plant and Animal Ecology, Ural Branch of the Russian Academy of Sciences, Yekaterinburg, Russia 22Ural Federal University, Yekaterinburg, Russia 23Institute of Geography, Russian Academy of Sciences, Moscow 119017, Russia 24Department of Earth Sciences, Natural History Museum, London, UK 25Institute of Evolutionary Biology (UPF-CSIC), Barcelona, Spain 26Centre Nacional d’Anàlisi Genòmica, Centre for Genomic Regulation (CNAG-CRG), The Barcelona Institute of Science and Technology, Barcelona, Spain 27Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain 28Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Barcelona, Spain 29Department of Applied Ecology, North Carolina State University, Raleigh, NC, USA 30Howard Hughes Medical Institute, University of California, Santa Cruz, Santa Cruz, CA 96050, USA 31Villum Center for Biodiversity Genomics, Section for Ecology and Evolution, Department of Biology, University of Copenhagen, Copenhagen, Denmark 32State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming 650223, China 33Center for Excellence in Animal Evolution and Genetics, Chinese Academy of Sciences, Kunming 650223, China 34Institut des Sciences de l’Évolution, Université Montpellier, CNRS, IRD, EPHE, Montpellier 34095, France 35Norwegian University of Science and Technology (NTNU) University Museum, Trondheim 7012, Norway 36Senior author 37Lead contact *Correspondence: shanlin.liu@cau.edu.cn (S.L.), love.dalen@nrm.se (L.D.), tgilbert@sund.ku.dk (M.T.P.G.) https://doi.org/10.1016/j.cell.2021.07.032 during which additional extinctions occurred. These consist of RESULTS the five extant species, as well as the now-extinct Siberian uni- corn (Elasmotherium sibiricum), Merck’s rhinoceros (Stephano- Resolving the rhinoceros phylogeny rhinus kirchbergensis) and its relative Stephanorhinus hemitoe- Three hypotheses have been proposed to explain the phyloge- chus (not studied here), and the woolly rhinoceros (Coelodonta netic relationships within living Rhinocerotidae: (1) the ‘‘horn hy- antiquitatis). pothesis,’’ which groups the two-horned rhinoceros species Despite decades of study, fundamental questions remain together, specifically placing the Sumatran rhinoceros as sister regarding the evolutionary relationships among the extant rhi- to the African Diceroti (black and white rhinoceroses) and has noceros species and their recently extinct relatives. Further- been supported by morphology (Antoine et al., 2010), genetic more, several rhinoceros species lack available genomic re- (e g. Steiner and Ryder, 2011), and paleoproteomic analyses of sources that would allow applications, including DNA-based dental enamel (Cappellini et al., 2019); (2) the ‘‘geographical hy- monitoring, conservation management, and environmental pothesis,’’ which places the Asian species together, with Suma- DNA studies. To address these questions and needs, we tran as sister to the greater one-horned and Javan rhinoceroses, analyzed a genome dataset representing eight rhinoceros spe- and is based on morpho-anatomical evidence (Antoine et al., cies (Figure 1), including all seven genera that survived into the 2021), biogeographic parsimony, genetic analyses using a Late Pleistocene (Cerdeño, 1998). Our data include the five limited number of loci (Kirillova et al., 2017; Kosintsev et al., extant rhinoceros species represented by four de novo genome 2019; Orlando et al., 2003; Tougard et al., 2001), and paleopro- assemblies of black (Diceros bicornis), white (Ceratotherium si- teomic analysis using collagen sequences (Welker et al., 2017); mum), Sumatran (Dicerorhinus sumatrensis), and greater one- and (3) a hypothesis that the Sumatran rhinoceros is sister to horned (Rhinoceros unicornis, also known as Indian) rhinocer- the clade comprising the four other extant species, which has oses and a resequenced genome of a Javan rhinoceros been supported by a more recent analysis of complete mito- (R. sondaicus). The Javan rhinoceros genome was retrieved chondrial genomes (Margaryan et al., 2020). These conflicting from a museum specimen dating to 1838 and resequenced hypotheses emphasize the limitations of using lower-resolution to high coverage (253). In addition, we sequenced the ge- markers in reconstructing evolutionary relationships within Rhi- nomes of three extinct rhinoceros species from Late Pleisto- nocerotidae and highlight the potential of applying phyloge- cene fossils that are close to, or beyond, the radiocarbon nomic approaches. Prior studies have also debated the phyloge- dating limit of 50 thousand years ago (kya), specifically a netic placement of the three extinct species included in this Siberian unicorn, a Merck’s rhinoceros, and a woolly rhinoc- study. For example, the relationship of Merck’s and woolly eros, sequenced to 93, 123, and 353 coverage, respectively rhinoceroses to each other, the Sumatran rhinoceros, and the (Table S1). two African Diceroti remained contentious due to contrasting 2 Cell 184, 1–12, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS Figure 1. Ranges of the eight rhinoceros species studied The historical distribution range maps of the five extant species use information published previously (Dinerstein, 2003; Havmøller et al., 2016; Khan and van Strien, 1997; Rookmaaker and Antoine, 2012), while the ranges of the three extinct species were drawn based on their fossil records (Kahlke and Lacombat, 2008; Kosintsev et al., 2019, 2020; Shpansky and Boeskorov, 2018). All ranges are approximations aimed at conveying each species’ general, rather than detailed, distribution. As the greater one-horned rhinoceros sample derives from captive-born zoo stock, a geographic origin is not shown. conclusions drawn from morpho-anatomical versus mitoge- Kosintsev et al., 2019). Within the DCS clade, we found strong nomic and paleoproteomic evidence (Antoine et al., 2021; Cap- support for the Sumatran rhinoceros as sister to the clade that pellini et al., 2019; Kirillova et al., 2017; Kosintsev et al., 2019). includes the extinct Merck’s and woolly rhinoceroses. A remain- To circumvent the problem of reference genome biases (Go- ing challenge is understanding the relationships of more ancient palakrishnan et al., 2017; Heintzman et al., 2017; Figures S1A extinct species, for which DNA remains unrecoverable. As our and S1B), we conducted whole-genome alignments for the eight findings suggest that results based purely on morphology (the rhinoceros species using domestic horse (Equus caballus) as an horn hypothesis) are not supported, attempting to fill in the rest outgroup. We then inferred a genome-wide species tree summa- of the phylogeny based solely on the morphology of extinct rizing the phylogenetic signal from individual gene trees based taxa may prove difficult. on 22,066 100-kb genomic windows. Our phylogenomic analysis identified three major clades within the subfamily Rhinocerotinae Gene flow among species and provided strong support for the geographical hypothesis of While we were able to resolve a genome-wide rhinoceros spe- rhinoceros evolution. A clade comprising the two African species cies tree, we uncovered significant phylogenetic discordance Diceros bicornis and Ceratotherium simum, the Diceroti, is the across the rhinoceros genomes, suggesting gene flow or incom- sister lineage to the remaining rhinoceroses in our dataset (with plete lineage sorting (ILS) among taxa. While this topology repre- the exception of the Siberian unicorn). A second clade includes sents the dominant signal of species relationships across the the Sumatran, Merck’s, and woolly rhinoceroses (hereafter whole genome, analyses of individual chromosomes did not al- referred to as the Dicerorhinus-Coelodonta-Stephanorhinus ways recover the same topology (Figure S1C). Most prominently, [DCS] clade), all of which have current or past geographic ranges we observed that the species-tree position of the DCS clade was that include parts of Asia. The third clade includes the two Rhi- supported by only 45% of the individual sliding-window trees, noceros species (Figure 2A). Thus, the principal divergence which is substantially fewer than recovered for other nodes in the among the rhinoceros lineages is related to the geographical di- phylogeny (Figure 2B). Because a genomic region of 100 kb in vision between species on the African and Eurasian continents. length may contain multiple recombination breakpoints, we Our phylogenomic analyses confirm prior conclusions based also inferred gene trees using 5-kb alignments randomly sub- on morphological and biomolecular evidence (Figure 2A) that sampled from within each 100-kb sliding window. The results place the extinct Siberian unicorn as outgroup to the subfamily corroborated the phylogenetic discordance discovered in the Rhinocerotinae (P.-O.A., unpublished data; Becker et al., 2013; initial dataset (Figure S1D). We then simulated the gene tree Cell 184, 1–12, September 16, 2021 3 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article Figure 2. Phylogeny of the Rhinocerotidae (A) Dated species tree of the Rhinocerotidae based on a consensus of trees generated every 100 kb across the genome using maximum likelihood methods and multiple fossil calibrations as detailed in STAR Methods. Blue horizontal bars show 95% confidence intervals of estimated divergence dates between lineages. Black crosses indicate extinct species. (B and C) Frequency of the three bipartitions induced by branches 4 and 5 (B) in the tree (C). Numbers on the x axis text and header (B) correspond to the branch identity of the tree (C). For example 4,7|8,9 represents the bipartition ((horse and tapir, Elasmotherium) Diceroti) | (Rhinoceros, DCS) induced by branch 5. distribution expected under a hypothesis of genetic drift with ILS logenies (Willerslev et al., 2009; Margaryan et al., 2020; and compared this to the empirical data. We found disagreement Figure 2B). between the simulated and empirical data, suggestive of both To further explore the origins of the phylogenomic discor- ILS and gene flow as the cause of the diminished support for dances, we used D-statistics analyses to differentiate ancient the species tree. We found in the empirical data an excess of gene flow events from ILS (Green et al., 2010). The use of this gene tree topologies supporting the placement of the Sumatran method can be problematic when applied to relatively divergent rhinoceros as sister to the other species of the Rhinocerotinae species due both to possible violations in assumptions such as subfamily, as observed previously in some mitochondrial phy- equal mutation rates and the infinite sites model and biases 4 Cell 184, 1–12, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS A B Figure 3. Comparison of whole-genome heterozygosity estimated in various taxa and ROH distribution in the eight rhinoceros species (A) Heterozygosity estimates of a broad range of animals, mammals, and ruminants. Species with heterozygosity values >1% were not included. See also Table S3. (B) Runs of homozygosity (ROH) size distributions for seven of the species investigated. We did not identify any ROH in the Siberian unicorn; data are not shown here, as we cannot exclude the influ- ence of reference bias. derived from the mapping of resequenced data to a reference Drivers of low genetic diversity genome (Figures S2A and S2B). Despite these potential sources Previous genomic studies on black, white, and Sumatran rhinoc- of error, the D-statistics are congruent with the phylogenomic eroses identified low levels of genetic diversity (Mays et al., 2018; analyses in suggesting that gene flow and ILS occurred between Moodley et al., 2020; Tunstall et al., 2018). These findings are the ancestors of the Diceroti and Rhinoceros species (Figures 2B consistent with the observation that all extant rhinoceros species and 2C). This gene flow may have been enabled by the Eurasian have gone through recent population size declines, even though origin of both African species (Antoine et al., 2021; Geraads, some species (white and greater one-horned rhinoceroses) have 2005, 2020). Our analyses also revealed an excess of shared since recovered (Ellis and Talukdar, 2019; Emslie, 2020; Rook- derived alleles between the two extinct members of the DCS maaker and Antoine, 2012). However, low genetic diversity can clade (Merck’s and woolly rhinoceroses) and both representa- also be a consequence of particular life-history traits and/or tives of Rhinoceros, suggesting gene flow between these two long-term small population size (Westbury et al., 2018, 2019; pairs of lineages. We found no evidence of excess shared Xue et al., 2015). To investigate this, we calculated genome- derived alleles between the Sumatran rhinoceros and either Rhi- wide heterozygosity (GWH) for all eight rhinoceros species and noceros species (Figure S2D), despite their closer geographic compared these estimates with GWH in a range of other animal proximity (Figure 1). This suggests either no gene flow or similar species, including ruminants and, more broadly, mammals. We levels of gene flow between the Sumatran rhino and both Rhi- assessed whether GWH levels are lower in genomes recovered noceros species. These differential patterns of gene flow may from present-day animals (i.e., black, white, greater one-horned, explain the discrepancy of the phylogenetic placement of the and Sumatran rhinoceroses) compared to GWH in genomes Sumatran rhinoceros as sister to the genus Rhinoceros based recovered from specimens that pre-date the human-mediated on nuclear DNA versus sister to all other members of the subfam- declines during the last 100 years (i.e., the nearly 180-year-old ily Rhinocerotinae in some mitochondrial phylogenies (Willerslev Javan rhinoceros genome as well as the genomes from the three et al., 2009; Margaryan et al., 2020). extinct species). We estimated GWH based on transversions only to limit the Timing of divergence between species potential influence of DNA damage on estimates from the We used fossil data to calibrate our phylogeny and estimate ancient and historical genomes. However, for comparability lineage divergence times (Figure 2A). This resulted in a 65 with published results for other taxa, which incorporate all vari- mya estimate for the common ancestor of horses, tapirs, and able sites, we recalibrated our estimates based on the expected rhinoceroses and a 36 mya estimate for the common transition/transversion ratio (see Figure S3A). Our results ancestor of the extinct rhinoceros subfamily Elasmotheriinae showed that present-day rhinoceros genomes exhibit signifi- and the extant subfamily Rhinocerotinae. The three major cantly lower GWH compared to the historical Javan and extinct clades within the Rhinocerotinae subfamily diverged 16 mya genomes (one-way ANOVA, n = 8, F = 7.4, p = 0.04). On the other (Figure 2A), at the end of the early Miocene and around the hand, our comparison with a broad range of animals shows that time of the Miocene climatic optimum (17–14 mya), a period rhinoceroses in general display comparatively low levels of that was 3 C–4 C warmer than present (Lewis et al., 2008; GWH, especially relative to not only the combined dataset of Sosdian et al., 2020). Diversification occurred after the forma- all animals but also ruminants and other large herbivores (Fig- tion of the terrestrial connection between the Afro-Arabian ure 3A). The only mammalian family displaying lower average and Eurasian landmasses 20 mya (Van Couvering and Del- levels of GWH was the Felidae (Figure S3B), which is not unex- son, 2020). We hypothesize that this land bridge enabled pected, as carnivores/predators are generally less abundant dispersal events followed by vicariance, as is well documented than herbivores/prey (Owen-Smith, 2015). These findings are with the immigration into Africa from Eurasia of early rhinocer- robust to choice of reference genome used in our analyses otids, giraffids, suids, and viverrids and the emigration from (see Figure S3C). Africa to Eurasia of apes, deinotheres, and elephantoids, To better contextualize the observed levels of GWH, we char- among others (Van Couvering and Delson, 2020). acterized the inbreeding levels in our genomes through Cell 184, 1–12, September 16, 2021 5 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article distributions of runs of homozygosity (ROH). To evaluate the Previous studies have suggested that maintenance of a low robustness of our results, we first explored the effect of population size over extended periods of time allows for the excluding transitions on ROH inference (Figure S4A). This anal- purging of deleterious alleles while keeping low genome-wide ysis indicated that a reliance solely on transversions shifts the levels of genetic diversity (Westbury et al., 2018, 2019; Xue distribution of ROH segments to longer stretches but that this et al., 2015). Our finding that all rhinoceros species have had a is only a problem for species with low overall GWH. Since all small Ne during extended periods of their history could thus indi- the resequenced genomes (Javan, Siberian unicorn, Merck’s, cate a similar scenario. and woolly rhinoceroses) exhibit higher GWH than at least three To investigate how mutational load in Rhinocerotidae com- of the de novo assembled genomes (black, white, Sumatran, and pares to other species, we calculated ratios of loss-of-function greater one-horned rhinoceros), a reliance solely on transver- (LoF; generally highly deleterious) versus synonymous mutations sions for these individuals may somewhat artificially inflate the for seven of the rhinoceros species and 30 mammalian species length of ROH segments but should not bias our overall from diverse clades (Table S4). The results show that the levels interpretations. of mutational load in rhinoceroses fall within the range observed We did not detect any ROH segments in the Siberian unicorn, among other present-day mammals (Figure 5). However, the which may reflect its exceptional phylogenetic distance to the extinct and historic rhinoceros genomes (Siberian unicorn, white rhinoceros genome against which it was mapped, artifi- woolly, and Javan rhinoceroses) displayed a significantly higher cially inflating heterozygous sites (Figure S3C). In contrast, number of LoF mutations (Figure 5; one-way ANOVA, n = 7, ROH segments >2 Mb were detected in all other rhinoceros spe- F = 29.0, p = 0.003) compared to the present-day rhinoceros ge- cies, where all species except the woolly rhinoceros also ex- nomes (i.e., black, white, Sumatran, and greater one-horned hibited stretches as long as 5 Mb. However, we observed signif- rhinoceroses). We thus find no evidence for an accumulation of icantly higher inbreeding levels (one-way ANOVA, n = 7, F = 36.7, mutational load within the last decades for those species that p = 0.002), which are proportional to the overall length of ROH have gone through recent population declines. Although specu- segments, in the genomes from present-day rhinoceroses lative, we therefore hypothesize that extant rhinoceroses may compared to the genomes of the historical Javan and extinct have undergone some purging of mutational load in connection Merck’s and woolly rhinoceroses (Figure 3B). with their demographic declines in the last 100 years. However, Overall, the comparisons of GWH and inbreeding levels sug- this hypothesis requires further testing, for example through gest that recent population declines caused by heavy anthro- intra-specific comparisons of historical and modern genomes pogenic pressure in the 20th century (e.g., Player, 1973) re- that span these declines and, in some cases, recoveries (Sán- sulted in marked losses in genetic diversity, as well as chez-Barreiro et al., 2021; von Seth et al., 2021). Meanwhile, it increased inbreeding levels. However, the genomes from his- is worth noting that some of the species used in the analysis torical and extinct species, which were sampled either prior were mapped onto their closely related species with available to their recent population collapse or many thousands of years reference genomes, and the qualities of genome assemblies before their extinction, also exhibit low levels of GWH when and annotations vary between species as well, both of which compared to other animal species (Figure 3A). Moreover, the could influence the accuracy of gene effect estimation. observation of a moderate amount of long ROH segments in the genomes from the Javan, Merck’s, and woolly rhinoceroses DISCUSSION is consistent with background inbreeding in these species. We thus hypothesize that limited genetic diversity and moderate Our combined rhinoceros genome dataset represents a valuable inbreeding levels are intrinsic features of rhinoceros life history, resource for both the current and future study of the evolution where low population densities and limited dispersal result in and biology of these species, including characterization of the increased genetic drift as well as occasional mating between genetic basis of rhinoceros phenotypes (Table S5). For example, relatives. we uncovered frameshift mutations in IFT43 (intraflagellar trans- port 43) that could contribute to rhinoceroses’ poor eyesight. Demography and mutational load IFT43 is involved in the formation and maintenance of cilia, which To further assess the genomic background of the overall low are important for the development and function of the light-sen- GWH and moderate inbreeding levels in Rhinocerotidae, we sitive tissue at the back of the eye (the retina) (Arts et al., 2011). modeled changes in effective population size (Ne) throughout In the case of the Javan and greater one-horned rhinoceros, the Pleistocene using the pairwise sequentially Markovian coa- our genome sequences also provide a basis for further species- lescent (PSMC) model (Li and Durbin, 2011; Figure 4). Although specific conservation genetic analyses. By analyzing this previous studies have reconstructed the demographic histories dataset, we resolved a long-standing debate related to the evolu- for a subset of species, including black rhinoceros (Moodley tionary history of living and recently extinct rhinoceroses and pro- et al., 2020), white rhinoceros (Tunstall et al., 2018), Sumatran vided evidence that relatively low genomic heterozygosity and rhinoceros (Mays et al., 2018), and woolly rhinoceros (Lord moderate inbreeding levels may represent their long-term natural et al., 2020), our combined analysis allows exploration of both state. These findings suggest that low levels of diversity and high shared and unique responses through time. Overall, it is striking inbreeding observed in present-day rhinoceros genomes can that all eight species displayed either a general continual only partially be attributed to recent declines. This may be posi- decrease in Ne over the last two million years or a continuously tive news for conservation, since it implies that recent declines small Ne over extended time periods. may have had less impact on the genetic aspects of population 6 Cell 184, 1–12, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS Figure 4. Demographic trajectory of the eight rhinoceros species Each curve represents one species, with thin lines depicting bootstraps. The x axis corresponds to time before present in years on a log scale. LGM, last glacial maximum. We assumed two different substitution rates (substitutions/site/generation, m) and generation times (g) for our rhinoceros species and set the sample age of the three extinct species to 50,000 years (see Table S1). The y axis corresponds to the effective population size (Ne). Species are grouped by geographic distribution: Africa, black and white rhinoceroses; South Asia, Sumatran, Javan, and greater one-horned rhinoceroses; and northern Eurasia, Siberian unicorn, Merck’s, and woolly rhinoceroses. viability than previously thought. Nonetheless, extant rhinocer- that needs to be bridged before we can fully understand the oses undoubtedly face enormous challenges in the future, princi- evolutionary history of the rhinoceros family. pally due to anthropogenic and environmental effects. A major priority for rhinoceros conservation will be to halt illegal poaching STAR+METHODS and ensure that there is sufficient carrying capacity for population recovery. Our study highlights how genomics can complement Detailed methods are provided in the online version of this paper such actions by enabling monitoring of ongoing changes in ge- and include the following: netic variation, inbreeding, and mutational load. d KEY RESOURCES TABLE Limitations of the study d RESOURCE AVAILABILITY Given the historic and ancient nature of the specimens from four B Lead contact of the species studied, their DNA quality was not suitable for de B Materials availability novo assembly; thus, their genome sequences were recovered B Data and code availability through mapping against other species. This process can intro- d EXPERIMENTAL MODEL AND SUBJECT DETAILS duce biases in downstream analyses that can arise due to differ- B Source organisms ential mapping efficiencies influenced by phylogenetic distance d METHOD DETAILS to the reference genome, ancient DNA damage, and short read B DNA isolation and sequencing lengths. Therefore, although we took several steps to alleviate B Genome assembly and annotation such influences, we highlight that this should be kept in mind. B Genome alignment and phylogenetic inference Furthermore, it should be noted that the rhinoceros family was B Molecular dating once a speciose group, and only a small fraction of these were B Alignment tool selection studied here. Thus, given that we lack genome sequences B Gene flow from the majority of the clade, clearly a huge gap will remain B Heterozygosity estimation Cell 184, 1–12, September 16, 2021 7 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article Least concern, Data deficient, Near threatened Vulnerable, Endangered, Critically endangered LoF mutation rate - missense/silent Figure 5. Genome variant effect of a broad range of mammal species The black line represents a linear regression fit that showed a significant correlation between loss-of-function (LoF) mutation rate and rate of missense/silent (r = 0.72, p = 4.05e-7). Note that the Javan rhinoceros sample dates to 1838, so we set its conservation status as data deficient. The rhinoceros species are shown in turquoise lettering. LoF mutations are here represented by nonsense mutations, and the rate is used as a proxy for the accumulation of mutation load. Missense/ silent means the rate of number of missense mutations relative to that of silent mutations. For species with multiple individuals, this figure shows their mean values. Note that the result of Merck’s rhinoceros was not included due to its abnormal transversion/transition ratio (Figure S4B). B Demographic inference and Runs of Homozygosity Ekaterinburg) for providing the sample of Siberian unicorn. M.T.P.G. was sup- (ROH) estimation ported by European Research Council (ERC) Consolidator grant 681396 B Genetic load and identification for rhinoceros specific (Extinction Genomics). E.D.L. was supported by Independent Research Fund Denmark grant 8021-00218B. A.C. was supported by an Australian frameshift mutation Research Council Laureate Fellowship (FL140100260). T.M.B. is supported B Test of heterozygosity, PSMC and ROH estimation for by funding from the ERC under the European Union’s Horizon 2020 research the non-modern rhino samples and innovation program (grant agreement 864203), grant BFU2017-86471-P d QUANTIFICATION AND STATISTICAL ANALYSIS (MINECO/FEDER, UE), ‘‘Unidad de Excelencia Marı́a de Maeztu’’ funded by d ADDITIONAL RESOURCES the AEI (CEX2018-000792-M), Howard Hughes International Early Career, and Secretaria d’Universitats i Recerca and CERCA Programme del Departa- ment d’Economia i Coneixement de la Generalitat de Catalunya (GRC 2017 SUPPLEMENTAL INFORMATION SGR 880). L.D. was supported by the Swedish Research Council (2017- 04647) and Formas (2018-01640). We thank Dmitry Bogdanov and Roger Supplemental information can be found online at https://doi.org/10.1016/j.cell. Hall for giving us permission to use their rhinoceros artwork. 2021.07.032. AUTHOR CONTRIBUTIONS ACKNOWLEDGMENTS S.L., L.D., and M.T.P.G. conceived the project and designed the research. The authors acknowledge support from the Science for Life Laboratory, the M.-H.S.S., K.J.M., S.V., P.K., I.K., A.C., B.S., and G.Z. provided archaeological Garvan Institute of Medical Research, the Knut and Alice Wallenberg Founda- work, logistics, and/or ancient samples and data. F.S.-B., Y.M., M.B., T.v.d.V., tion, and the National Genomics Infrastructure funded by the Swedish O.R., C.S., and L.G.R.B.-v.S. coordinated logistics of and/or provided modern Research Council and Uppsala Multidisciplinary Center for Advanced Compu- samples and data. M.-H.S.S., N.D., K.J.M., P.D.H., J.D.K., J.v.S., H.H., C.G., tational Science for assistance with massively parallel sequencing and access G.M., and C.Y. conducted laboratory work. S.L., R.A.-O., D.A.D., M.V.W., and to the UPPMAX computational infrastructure. We thank the Natural History L.C. conducted analyses of data with considerable input from A.M., T.v.d.V., Museum at the University of Oslo for providing the Javan rhinoceros sample. S.G., P.D.H., T.M.-B., P.-O.A., L.D., and M.T.P.G. S.L., M.V.W., N.D., We thank the Museum of the Institute of Plant and Animal Ecology (UB RAS, K.J.M., P.D.H., D.A.D., F.S.-B., A.M., B.D.C., Y.M., K.R., A.L., T.M.-B., S.G., 8 Cell 184, 1–12, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS E.D.L., R.R.D., B.S., P.-O.A., L.D., and M.T.P.G. interpreted results and wrote ulation extinctions in a small mammal indicate Late Pleistocene ecosystem the paper with input from all other authors. instability. Proc. Natl. Acad. Sci. USA 109, 20532–20536. Briggs, A.W., Stenzel, U., Johnson, P.L.F., Green, R.E., Kelso, J., Prüfer, K., DECLARATION OF INTERESTS Meyer, M., Krause, J., Ronan, M.T., Lachmann, M., and Pääbo, S. (2007). Pat- terns of damage in genomic DNA sequences from a Neandertal. Proc. Natl. The authors declare no competing interests. Acad. Sci. USA 104, 14616–14621. Briggs, A.W., Stenzel, U., Meyer, M., Krause, J., Kircher, M., and Pääbo, S. Received: March 26, 2021 (2010). Removal of deaminated cytosines and detection of in vivo methylation Revised: June 16, 2021 in ancient DNA. Nucleic Acids Res. 38, e87. Accepted: July 23, 2021 Published: August 24, 2021 Brotherton, P., Endicott, P., Sánchez, J.J.M., Beaumont, M., Barnett, R., Aus- tin, J., and Cooper, A. (2007). Novel high-resolution characterization of ancient SUPPORTING CITATIONS DNA reveals C > U-type base modification events as the sole cause of post mortem miscoding lesions. Nucleic Acids Res. 35, 5717–5728. The following references appear in the supplemental information: Carbone Brown, F.H., and McDougall, I. (2011). Geochronology of the Turkana depres- et al. (2014); Chen et al. (2019); de Manuel et al. (2016); Dong et al. (2017); Du- sion of northern Kenya and southern Ethiopia. Evol. Anthropol. 20, 217–227. rand et al. (2011); Freedman et al. (2014); Groenen et al. (2012); Harr et al. (2016); Kardos et al. (2018); Karsten et al. (2011); Liu et al. (2014); Lynch Buckley, M. (2015). Ancient collagen reveals evolutionary history of the et al. (2015); Mallick et al. (2016); Nater et al. (2017); Palkopoulou et al. endemic South American ‘ungulates’. Proc. Biol. Sci. 282, 20142671, (2018); Prado-Martinez et al. (2013); (Reddy et al., 2015); Robinson et al. 20142671. (2016); Smith et al. (2017); Svardal et al. (2017); Teng et al. (2017); Veeramah Cappellini, E., Jensen, L.J., Szklarczyk, D., Ginolhac, A., da Fonseca, R.A., et al. (2015); Wang et al. (2019); Wilkie et al. (2013); Yim et al. (2014); Yoo Stafford, T.W., Holen, S.R., Collins, M.J., Orlando, L., Willerslev, E., et al. et al. (2017); (Zhang et al., 2014); Zoonomia Consortium (2020); (Cho (2012). Proteomic analysis of a pleistocene mammoth femur reveals more et al., 2013) than one hundred ancient bone proteins. J. Proteome Res. 11, 917–926. Cappellini, E., Welker, F., Pandolfi, L., Ramos-Madrigal, J., Samodova, D., REFERENCES Rüther, P.L., Fotakis, A.K., Lyon, D., Moreno-Mayar, J.V., Bukhsianidze, M., et al. (2019). Early Pleistocene enamel proteome from Dmanisi resolves Ste- Antoine, P.-O., Downing, K.F., Crochet, J., Duranthon, F., Flynn, L.J., Mari- phanorhinus phylogeny. Nature 574, 103–107. vaux, L., Metais, G., Rajpar, A.R., and Roohi, G. (2010). A revision of Acerathe- Carbone, L., Harris, R.A., Gnerre, S., Veeramah, K.R., Lorente-Galdos, B., rium blanfordi Lydekker, 1884 (Mammalia: Rhinocerotidae) from the Early Huddleston, J., Meyer, T.J., Herrero, J., Roos, C., Aken, B., et al. (2014). Miocene of Pakistan: postcranials as a key. Zool. J. Linn. Soc. 160, 139–194. Gibbon genome and the fast karyotype evolution of small apes. Nature 513, Antoine, P.-O., Métais, G., Orliac, M.J., Crochet, J.-Y., Flynn, L.J., Marivaux, 195–201. L., Rajpar, A.R., Roohi, G., and Welcomme, J.-L. (2013). Mammalian Neogene biostratigraphy of the Sulaiman Province, Pakistan. In Fossil mammals of Asia: Caroe, C., Gopalakrishnan, S., Vinner, L., Mak, S.S.T., Sinding, M.S., Sama- Neogene Biostratigraphy and Chronology, X.-m. Wang, L.J. Flynn, and M. For- niego, J.A., Wales, N., Sicheritzponten, T., and Gilbert, M.T.P. (2017). Sin- telius, eds. (Columbia University Press), pp. 400–422. gle-tube library preparation for degraded DNA. Methods Ecol. Evol. 9, 410–419. Antoine, P.-O., Reyes, M.C., Amano, N., Claude, J., Bautista, A.P., Vos, J.d., and Ingicco, T. (2021). A new clade of rhinoceroses from the Pleistocene of Cerdeño, E. (1998). Diversity and evolutionary trends of the Family Rhinocer- Asia sheds light on mainland mammal dispersals to the Philippines. Zool. J. otidae (Perissodactyla). Palaeogeogr. Palaeoclimatol. Palaeoecol. 141, 13–34. Linn. Soc. Published online May 22, 2021. https://doi.org/10.1093/zoolin- Chen, L., Qiu, Q., Jiang, Y., Wang, K., Lin, Z., Li, Z., Bibi, F., Yang, Y., Wang, J., nean/zlab009. Nie, W., et al. (2019). Large-scale ruminant genome sequencing provides in- Arts, H.H., Bongers, E.M.H.F., Mans, D.A., van Beersum, S.E., Oud, M.M., sights into their evolution and distinct traits. Science 364, 364. Bolat, E., Spruijt, L., Cornelissen, E.A.M., Schuurs-Hoeijmakers, J.H.M., de Cho, Y.S., Hu, L., Hou, H., Lee, H., Xu, J., Kwon, S., Oh, S., Kim, H.-M., Jho, S., Leeuw, N., et al. (2011). C14ORF179 encoding IFT43 is mutated in Sense- Kim, S., et al. (2013). The tiger genome and comparative analysis with lion and nbrenner syndrome. J. Med. Genet. 48, 390–395. snow leopard genomes. Nat. Commun. 4, 2433. Auwera, G.A.V.d., Carneiro, M.O., Hartl, C., Poplin, R., Angel, G.d., Lev- Cingolani, P., Platts, A., Wang, L., Coon, M., Nguyen, T., Wang, L., Land, S.J., y-Moonshine, A., Jordan, T., Shakir, K., Roazen, D., Thibault, J., et al. (2013). Lu, X., and Ruden, D.M. (2012). A program for annotating and predicting the From FastQ data to high-confidence variant calls: the Genome Analysis Toolkit effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of best practices pipeline. Curr. Protoc. Bioinformatics. 43, 11.10.11–11.10.33. Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin) 6, 80–92. Bai, B., Meng, J., Zhang, C., Gong, Y.-X., and Wang, Y.-Q. (2020). The origin of Dabney, J., Knapp, M., Glocke, I., Gansauge, M.T., Weihmann, A., Nickel, B., Rhinocerotoidea and phylogeny of Ceratomorpha (Mammalia, Perissodac- Valdiosera, C., Garcı́a, N., Pääbo, S., Arsuaga, J.L., and Meyer, M. (2013). tyla). Commun. Biol. 3, 509. Complete mitochondrial genome sequence of a Middle Pleistocene cave Becker, D., Antoine, P., and Maridet, O. (2013). A new genus of Rhinocerotidae bear reconstructed from ultrashort DNA fragments. Proc. Natl. Acad. Sci. (Mammalia, Perissodactyla) from the Oligocene of Europe. J. Syst. Palaeontol- USA 110, 15758–15763. ogy 11, 947–972. F. Darwin, ed. (1909). The foundations of the origin of species. Two essays Blanchette, M., Kent, W.J., Riemer, C., Elnitski, L., Smit, A.F., Roskin, K.M., written in 1842 and 1844 (Cambridge University Press). Baertsch, R., Rosenbloom, K., Clawson, H., Green, E.D., et al. (2004). Aligning multiple genomic sequences with the threaded blockset aligner. Genome Res. de Manuel, M., Kuhlwilm, M., Frandsen, P., Sousa, V.C., Desai, T., Prado-Mar- 14, 708–715. tinez, J., Hernandez-Rodriguez, J., Dupanloup, I., Lao, O., Hallast, P., et al. (2016). Chimpanzee genomic diversity reveals ancient admixture with bono- Böhme, M., Aiglstorfer, M., Antoine, P., Appel, E., Havlik, P., Metais, G., bos. Science 354, 477–481. Schneider, S., Setzer, F., Tappert, R., and Tran, D.N. (2013). Na Duong (northern Vietnam) – an exceptional window into Eocene ecosystems from Dinerstein, E. (2003). The return of the unicorns: the natural history and conser- Southeast Asia. Zitteliana 53, 121–167. vation of the greater one-horned rhinoceros (Columbia University Press). Brace, S., Palkopoulou, E., Dalén, L., Lister, A.M., Miller, R., Otte, M., Ger- Dinerstein, E., and McCracken, G.F. (1990). Endangered greater one-horned monpré, M., Blockley, S.P.E., Stewart, J.R., and Barnes, I. (2012). Serial pop- rhinoceros carry high levels of genetic variation. Conserv. Biol. 4, 417–422. Cell 184, 1–12, September 16, 2021 9 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article Dong, D., Lei, M., Hua, P., Pan, Y.-H., Mu, S., Zheng, G., Pang, E., Lin, K., and Heissig, K. (2012). The American genus Penetrigonias Tanner & Martin, 1976 Zhang, S. (2017). The genomes of two bat species with long constant fre- (Mammalia: Rhinocerotidae) as a stem group elasmothere and ancestor of Me- quency echolocation calls. Mol. Biol. Evol. 34, 20–34. noceras Troxell, 1921 (Zitteliana). dos Reis, M., Donoghue, P.C.J., and Yang, Z. (2014). Neither phylogenomic Holt, C., and Yandell, M. (2011). MAKER2: an annotation pipeline and genome- nor palaeontological data support a Palaeogene origin of placental mammals. database management tool for second-generation genome projects. BMC Biol. Lett. 10, 20131003. Bioinformatics 12, 491. Doyle, V.P., Young, R.E., Naylor, G.J.P., and Brown, J.M. (2015). Can we iden- Junier, T., and Zdobnov, E.M. (2010). The Newick utilities: high-throughput tify genes with Increased Phylogenetic Reliability. Syst. Biol. 64, 824–837. phylogenetic tree processing in the UNIX shell. Bioinformatics 26, 1669–1670. Durand, E.Y., Patterson, N., Reich, D., and Slatkin, M. (2011). Testing for Kahlke, R.-D., and Lacombat, F. (2008). The earliest immigration of woolly rhi- ancient admixture between closely related populations. Mol. Biol. Evol. 28, noceros (Coelodonta tologoijensis, Rhinocerotidae, Mammalia) into Europe 2239–2252. and its adaptive evolution in Palaearctic cold stage mammal faunas. Quat. Ellis, S., and Talukdar, B. (2019). Rhinoceros unicornis. In the IUCN red list of Sci. Rev. 27, 1951–1961. threatened species 2019: eT19496A18494149 (IUCN). https://dx.doi.org/10. Kalyaanamoorthy, S., Minh, B.Q., Wong, T.K.F., von Haeseler, A., and Jermiin, 2305/IUCN.UK.2019-3.RLTS.T19496A18494149.en. L.S. (2017). ModelFinder: fast model selection for accurate phylogenetic esti- Emslie, R. (2020). Ceratotherium simum. In the IUCN red list of threatened spe- mates. Nat. Methods 14, 587–589. cies 2020: eT4185A45813880 (IUCN). https://dx.doi.org/10.2305/IUCN.UK. Kardos, M., Åkesson, M., Fountain, T., Flagstad, Ø., Liberg, O., Olason, P., 2020-1.RLTS.T4185A45813880.en. Sand, H., Wabakken, P., Wikenros, C., and Ellegren, H. (2018). Genomic con- Ersmark, E., Orlando, L., Sandoval-Castellanos, E., Barnes, I., Barnett, R., Stu- sequences of intensive inbreeding in an isolated wolf population. Nat. Ecol. art, A.J., Lister, A.M., and Dalen, L. (2015). Population demography and ge- Evol. 2, 124–131. netic diversity in the Pleistocene cave lion. Open Quat. 1, 1–14. Karsten, M., Van Vuuren, B.J., Goodman, P.S., and Barnaud, A. (2011). The Freedman, A.H., Gronau, I., Schweizer, R.M., Ortega-Del Vecchyo, D., Han, E., history and management of black rhino in KwaZulu-Natal: A population genetic Silva, P.M., Galaverni, M., Fan, Z., Marx, P., Lorente-Galdos, B., et al. (2014). approach to assess the past and guide the future. Anim. Conserv. 14, 363–370. Genome sequencing highlights the dynamic early history of dogs. PLoS Genet. Katoh, K., and Standley, D.M. (2013). MAFFT multiple sequence alignment 10, e1004016. software version 7: improvements in performance and usability. Mol. Biol. Geraads, D. (2005). Pliocene Rhinocerotidae (Mammalia) from Hadar and Di- Evol. 30, 772–780. kika (Lower Awash, Ethiopia), and a revision of the origin of modern african Khan, M.K.B.M., and van Strien, N.J. (1997). Asian rhinos: status survey and rhinos. J. Vertebr. Paleontol. 25, 451–461. conservation action planVolume 32 (IUCN). Geraads, D. (2020). Perissodactyla (Rhinocerotidae and Equidae) from Kana- Kie1basa, S.M., Wan, R., Sato, K., Horton, P., and Frith, M.C. (2011). Adaptive poi. J. Hum. Evol. 140, 102373. seeds tame genomic sequence comparison. Genome Res. 21, 487–493. Gilbert, M.T.P., Tomsho, L.P., Rendulic, S., Packard, M., Drautz, D.I., Sher, A., Kirillova, I.V., Chernova, O.F., Der Made, J.V., Kukarskih, V.V., Shapiro, B., Der Tikhonov, A., Dalén, L., Kuznetsova, T., Kosintsev, P., et al. (2007). Whole- Plicht, J.V., Shidlovskiy, F.K., Heintzmann, P.D., Van Kolfschoten, T., and genome shotgun sequencing of mitochondria from ancient hair shafts. Sci- Zanina, O.G. (2017). Discovery of the skull of Stephanorhinus kirchbergensis ence 317, 1927–1930. (Jäger, 1839) above the Arctic Circle. Quat. Res. 88, 537–550. Gnerre, S., Maccallum, I., Przybylski, D., Ribeiro, F.J., Burton, J.N., Walker, Korf, I. (2004). Gene finding in novel genomes. BMC Bioinformatics 5, 59, 59. B.J., Sharpe, T., Hall, G., Shea, T.P., Sykes, S., et al. (2011). High-quality draft Korneliussen, T.S., Albrechtsen, A., and Nielsen, R. (2014). ANGSD: analysis of assemblies of mammalian genomes from massively parallel sequence data. next generation sequencing data. BMC Bioinformatics 15, 356. Proc. Natl. Acad. Sci. USA 108, 1513–1518. Kosintsev, P., Mitchell, K.J., Devièse, T., van der Plicht, J., Kuitems, M., Pet- Gopalakrishnan, S., Samaniego Castruita, J.A., Sinding, M.S., Kuderna, rova, E., Tikhonov, A., Higham, T., Comeskey, D., Turney, C., et al. (2019). Evo- L.F.K., Räikkönen, J., Petersen, B., Sicheritz-Ponten, T., Larson, G., Orlando, lution and extinction of the giant rhinoceros Elasmotherium sibiricum sheds L., Marques-Bonet, T., et al. (2017). The wolf reference genome sequence light on late Quaternary megafaunal extinctions. Nat. Ecol. Evol. 3, 31–38. (Canis lupus lupus) and its implications for Canis spp. population genomics. BMC Genomics 18, 495. Kosintsev, P.A., Zykov, S.V., Tiunov, M.P., Shpansky, A.V., Gasilin, V.V., Gim- ranov, D.O., and Devjashin, M.M. (2020). The first find of Merck’s rhinoceros Green, R.E., Krause, J., Briggs, A.W., Maricic, T., Stenzel, U., Kircher, M., Pat- (Mammalia, Perissodactyla, Rhinocerotidae, Stephanorhinus kirchbergensis terson, N., Li, H., Zhai, W., Fritz, M.H.-Y., et al. (2010). A draft sequence of the Jäger, 1839) remains in the Russian far east. Dokl. Biol. Sci. 491, 47–49. Neandertal genome. Science 328, 710–722. Kuhn, R.M., Haussler, D., and Kent, W.J. (2013). The UCSC genome browser Groenen, M.A.M., Archibald, A.L., Uenishi, H., Tuggle, C.K., Takeuchi, Y., and associated tools. Brief. Bioinform. 14, 144–161. Rothschild, M.F., Rogel-Gaillard, C., Park, C., Milan, D., Megens, H.-J., et al. (2012). Analyses of pig genomes provide insight into porcine demography Lewis, A.R., Marchant, D.R., Ashworth, A.C., Hedenas, L., Hemming, S.R., and evolution. Nature 491, 393–398. Johnson, J.V., Leng, M.J., Machlus, M.L., Newton, A.E., and Raine, J.I.  Babiker, H., (2008). Mid-Miocene cooling and the extinction of tundra in continental Harr, B., Karakoc, E., Neme, R., Teschke, M., Pfeifle, C., Pezer, Z., Antarctica. Proceedings of the National Academy of Sciences of the United Linnenbrink, M., Montero, I., Scavetta, R., et al. (2016). Genomic resources for States of America 105, 10676–10680. wild populations of the house mouse, Mus musculus and its close relative Mus spretus. Sci. Data 3, 160075. Li, H. (2011). A statistical framework for SNP calling, mutation discovery, asso- ciation mapping and population genetical parameter estimation from Havmøller, R.G., Payne, J., Ramono, W., Ellis, S., Yoganand, K., Long, B., sequencing data. Bioinformatics 27, 2987–2993. Dinerstein, E., Williams, A.C., Putra, R.H., and Gawi, J. (2016). Will current conservation responses save the critically endangered Sumatran rhinoceros Li, H. (2013). Aligning sequence reads, clone sequences and assembly contigs Dicerorhinus sumatrensis? Oryx 50, 355–359. with BWA-MEM. arXiv. Genomics 1303, 3997. Heintzman, P.D., Zazula, G.D., MacPhee, R., Scott, E., Cahill, J.A., McHorse, Li, H. (2014). Toward better understanding of artifacts in variant calling from B.K., Kapp, J.D., Stiller, M., Wooller, M.J., Orlando, L., et al. (2017). A new high-coverage samples. Bioinformatics 30, 2843–2851. genus of horse from Pleistocene North America. eLife 6, e29944. Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with Bur- Heissig, K. (1972). Paläontologische und geologische Untersuchungen im Ter- rows-Wheeler transform. Bioinformatics 25, 1754–1760. tiär von Pakistan. V. Rhinocerotidae (Mamm.) aus den unteren und mittleren Li, H., and Durbin, R. (2011). Inference of human population history from indi- Siwalik-Schichten. Abh. Bayerischen Akad. Wiss. 152, 1–112. vidual whole-genome sequences. Nature 475, 493–496. 10 Cell 184, 1–12, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS Liu, S., Lorenzen, E.D., Fumagalli, M., Li, B., Harris, K., Xiong, Z., Zhou, L., Kor- Orlando, L., Ginolhac, A., Raghavan, M., Vilstrup, J., Rasmussen, M., Magnus- neliussen, T.S., Somel, M., Babbitt, C., et al. (2014). Population genomics sen, K., Steinmann, K.E., Kapranov, P., Thompson, J.F., Zazula, G., et al. reveal recent speciation and rapid evolutionary adaptation in polar bears. (2011). True single-molecule DNA sequencing of a pleistocene horse bone. Cell 157, 785–794. Genome Res. 21, 1705–1719. Lord, E., Dussex, N., Kierczak, M., Dı́ez-Del-Molino, D., Ryder, O.A., Stanton, Orlando, L., Leonard, J.A., Thenot, A., Laudet, V., Guérin, C., and Hänni, C. D.W.G., Gilbert, M.T.P., Sánchez-Barreiro, F., Zhang, G., Sinding, M.S., et al. (2003). Ancient DNA analysis reveals woolly rhino evolutionary relationships. (2020). Pre-extinction demographic stability and genomic signatures of adap- Mol. Phylogenet. Evol. 28, 485–499. tation in the woolly rhinoceros. Curr. Biol. 30, 3871–3879.e7. Owen-Smith, N. (2015). Mechanisms of coexistence in diverse herbivore– Luo, R., Liu, B., Xie, Y., Li, Z., Huang, W., Yuan, J., He, G., Chen, Y., Pan, Q., carnivore assemblages: demographic, temporal and spatial heterogeneities Liu, Y., et al. (2015). Erratum: SOAPdenovo2: an empirically improved mem- affecting prey vulnerability. Oikos 124, 1417–1426. ory-efficient short-read de novo assembler. Gigascience 4, 30, 30. Pääbo, S. (1989). Ancient DNA: extraction, characterization, molecular clon- Lynch, V.J., Bedoya-Reina, O.C., Ratan, A., Sulak, M., Drautz-Moses, D.I., ing, and enzymatic amplification. Proc. Natl. Acad. Sci. USA 86, 1939–1943. Perry, G.H., Miller, W., and Schuster, S.C. (2015). Elephantid genomes reveal Palkopoulou, E., Mallick, S., Skoglund, P., Enk, J., Rohland, N., Li, H., Omrak, the molecular bases of woolly mammoth adaptations to the arctic. Cell Rep. A., Vartanyan, S., Poinar, H., Götherström, A., et al. (2015). Complete genomes 12, 217–228. reveal signatures of demographic and genetic declines in the woolly mammoth. Curr. Biol. 25, 1395–1400. Mak, S.S.T., Gopalakrishnan, S., Carøe, C., Geng, C., Liu, S., Sinding, M.S., Kuderna, L.F.K., Zhang, W., Fu, S., Vieira, F.G., et al. (2017). Comparative per- Palkopoulou, E., Lipson, M., Mallick, S., Nielsen, S., Rohland, N., Baleka, S., formance of the BGISEQ-500 vs Illumina HiSeq2500 sequencing platforms for Karpinski, E., Ivancevic, A.M., To, T.H., Kortschak, R.D., et al. (2018). A palaeogenomic sequencing. Gigascience 6, 1–13. comprehensive genomic history of extinct and living elephants. Proc. Natl. Acad. Sci. USA 115, E2566–E2574. Mallick, S., Li, H., Lipson, M., Mathieson, I., Gymrek, M., Racimo, F., Zhao, M., Chennagiri, N., Nordenfelt, S., Tandon, A., et al. (2016). The Simons Genome Parks, M., and Lambert, D. (2015). Impacts of low coverage depths and post- Diversity Project: 300 genomes from 142 diverse populations. Nature 538, mortem DNA damage on variant calling: a simulation study. BMC Genomics 201–206. 16, 19, 19. Pickford, M., Senut, B., and Hadoto, D. (1993). Geology and palaeobiology of Margaryan, A., Sinding, M.S., Liu, S., Vieira, F.G., Chan, Y.L., Nathan, S.K.S.S., the Albertine Rift valley, Uganda-Zaire. Volume I: geology; Géologie et paléo- Moodley, Y., Bruford, M.W., and Gilbert, M.T.P. (2020). Recent mitochondrial biologie de la vallée du Rift albertin, Ouganda, Zaı̈re. Volume I: Géologie (CI- lineage extinction in the critically endangered Javan rhinoceros. Zool. J. Linn. FEG Occasional Publications). Soc. 190, 372–383. Player, I. (1973). The White Rhino Saga, 1st (Stein and Day). Marshall, D.C., Simon, C., and Buckley, T.R. (2006). Accurate branch length estimation in partitioned Bayesian analyses requires accommodation of Prado-Martinez, J., Sudmant, P.H., Kidd, J.M., Li, H., Kelley, J.L., Lorente-Gal- among-partition rate variation and attention to branch length priors. Syst. dos, B., Veeramah, K.R., Woerner, A.E., O’Connor, T.D., Santpere, G., et al. Biol. 55, 993–1003. (2013). Great ape genetic diversity and population history. Nature 499, 471–475. Mays, H.L., Jr., Hung, C.M., Shaner, P.J.L., Denvir, J., Justice, M., Yang, S.F., Prüfer, K., Munch, K., Hellmann, I., Akagi, K., Miller, J.R., Walenz, B., Koren, S., Roth, T.L., Oehler, D.A., Fan, J., Rekulapally, S., and Primerano, D.A. (2018). Sutton, G., Kodira, C., Winer, R., et al. (2012). The bonobo genome compared Genomic analysis of demographic history and ecological niche modeling in with the chimpanzee and human genomes. Nature 486, 527–531. the endangered Sumatran rhinoceros Dicerorhinus sumatrensis. Curr. Biol. 28, 70–76.e4. Reddy, P.C., Sinha, I., Kelkar, A., Habib, F., Pradhan, S.J., Sukumar, R., and Galande, S. (2015). Comparative sequence analyses of genome and transcrip- Mendes, F.K., and Hahn, M.W. (2016). Gene tree discordance causes tome reveal novel transcripts and variants in the Asian elephant Elephas max- apparent substitution rate variation. Syst. Biol. 65, 711–721. imus. J. Biosci. 40, 891–907. Meyer, M., and Kircher, M. (2010). Illumina sequencing library preparation for Robinson, J.A., Ortega-Del Vecchyo, D., Fan, Z., Kim, B.Y., vonHoldt, B.M., highly multiplexed target capture and sequencing. Cold Spring Harb. Protoc. Marsden, C.D., Lohmueller, K.E., and Wayne, R.K. (2016). Genomic flatlining 2010, pdb.prot5448. in the endangered island fox. Curr. Biol. 26, 1183–1189. Mirarab, S., Reaz, R., Bayzid, M.S., Zimmermann, T., Swenson, M.S., and Rookmaaker, K., and Antoine, P.-O. (2012). New maps representing the histor- Warnow, T. (2014). ASTRAL: genome-scale coalescent-based species tree ical and recent distribution of the African species of rhinoceros: Diceros bicor- estimation. Bioinformatics 30, i541–i548. nis, Ceratotherium simum and Ceratotherium cottoni. Pachyderm 52, 91–96. Moodley, Y., Westbury, M.V., Russo, I.M., Gopalakrishnan, S., Rakotoarivelo, Rose, K.D., Holbrook, L.T., Rana, R.S., Kumar, K., Jones, K.E., Ahrens, H.E., A., Olsen, R.-A., Prost, S., Tunstall, T., Ryder, O.A., Dalén, L., and Bruford, Missiaen, P., Sahni, A., and Smith, T. (2014). Early Eocene fossils suggest M.W. (2020). Interspecific gene flow and the evolution of specialisation in black that the mammalian order Perissodactyla originated in India. Nat. Commun. and white rhinoceros. Mol. Biol. Evol. 37, 3105–3117. 5, 5570. Muizon, C.d., and Ladevèze, S. (2020). Cranial anatomy of Andinodelphys co- Sánchez-Barreiro, F., Gopalakrishnan, S., Ramos-Madrigal, J., Westbury, chabambensis, a stem metatherian from the early Palaeocene of Bolivia. Geo- M.V., de Manuel, M., Margaryan, A., Ciucani, M.M., Vieira, F.G., Patramanis, diversitas 42, 597–739. Y., Kalthoff, D.C., et al. (2021). Historical population declines prompted signif- Nater, A., Mattle-Greminger, M.P., Nurcahyo, A., Nowak, M.G., de Manuel, M., icant genomic erosion in the northern and southern white rhinoceros. Mol. Desai, T., Groves, C., Pybus, M., Sonay, T.B., Roos, C., et al. (2017). Morpho- Ecol. Published online June 27, 2021. https://doi.org/10.1111/mec.16043. metric, behavioral, and genomic evidence for a new orangutan species. Curr. Sayyari, E., Whitfield, J.B., and Mirarab, S. (2018). DiscoVista: Interpretable vi- Biol. 27, 3576–3577. sualizations of gene tree discordance. Mol. Phylogenet. Evol. 122, 110–115. Nguyen, L.T., Schmidt, H.A., von Haeseler, A., and Minh, B.Q. (2015). IQ- Schubert, M., Lindgreen, S., and Orlando, L. (2016). AdapterRemoval v2: rapid TREE: a fast and effective stochastic algorithm for estimating maximum-likeli- adapter trimming, identification, and read merging. BMC Res. Notes 9, 88. hood phylogenies. Mol. Biol. Evol. 32, 268–274. Seppey, M., Manni, M., and Zdobnov, E.M. (2019). BUSCO: Assessing O’Leary, M.A., Bloch, J.I., Flynn, J.J., Gaudin, T.J., Giallombardo, A., Giannini, Genome Assembly and Annotation Completeness. Methods Mol. Biol. 1962, N.P., Goldberg, S.L., Kraatz, B.P., Luo, Z.X., Meng, J., et al. (2013). 227–245. The placental mammal ancestor and the post-K-Pg radiation of placentals. Shpansky, A.V., and Boeskorov, G.G. (2018). Northernmost Record of the Science 339, 662–667. Merck’s Rhinoceros Stephanorhinus kirchbergensis (Jäger) and Taxonomic Cell 184, 1–12, September 16, 2021 11 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article Status of Coelodonta jacuticus Russanov (Mammalia, Rhinocerotidae). Pale- Wang, G.-D., Zhang, M., Wang, X., Yang, M.A., Cao, P., Liu, F., Lu, H., Feng, ontol. J. 52, 445–462. X., Skoglund, P., Wang, L., et al. (2019). Genomic approaches reveal an Simpson, J.T., Wong, K., Jackman, S.D., Schein, J.E., Jones, S.J.M., and endemic subpopulation of gray wolves in southern China. iScience 20, Birol, I. (2009). ABySS: a parallel assembler for short read sequence data. 110–118. Genome Res. 19, 1117–1123. Welker, F., Collins, M.J., Thomas, J.A., Wadsley, M., Brace, S., Cappellini, E., Slater, G.S.C., and Birney, E. (2005). Automated generation of heuristics for Turvey, S.T., Reguero, M., Gelfo, J.N., Kramarz, A., et al. (2015). Ancient pro- biological sequence comparison. BMC Bioinformatics 6, 31. teins resolve the evolutionary history of Darwin’s South American ungulates. Smith, C.D., Edgar, R.C., Yandell, M.D., Smith, D.R., Celniker, S.E., Myers, Nature 522, 81–84. E.W., and Karpen, G.H. (2007). Improved repeat identification and masking Welker, F., Smith, G.M., Hutson, J.M., Kindler, L., Garcia-Moreno, A., Villa- in Dipterans. Gene 389, 1–9. luenga, A., Turner, E., and Gaudzinski-Windheuser, S. (2017). Middle Pleisto- Smith, S.D., Kawash, J.K., Karaiskos, S., Biluck, I., and Grigoriev, A. (2017). cene protein sequences from the rhinoceros genus Stephanorhinus and the Evolutionary adaptation revealed by comparative genome analysis of woolly phylogeny of extant and extinct Middle/Late Pleistocene Rhinocerotidae. mammoths and elephants. DNA Res. 24, 359–369. PeerJ 5, e3033. Sosdian, S.M., Babila, T.L., Greenop, R., Foster, G.L., and Lear, C.H. (2020). Westbury, M., Baleka, S., Barlow, A., Hartmann, S., Paijmans, J.L.A., Kramarz, Ocean Carbon Storage across the middle Miocene: a new interpretation for A., Forasiepi, A.M., Bond, M., Gelfo, J.N., Reguero, M.A., et al. (2017). A mito- the Monterey Event. Nat. Commun. 11, 134. genomic timetree for Darwin’s enigmatic South American mammal Macrau- Stamatakis, A. (2014). RAxML version 8: a tool for phylogenetic analysis and chenia patachonica. Nat. Commun. 8, 15951. post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. Westbury, M.V., Hartmann, S., Barlow, A., Wiesel, I., Leo, V., Welch, R., Parker, Stanke, M., Diekhans, M., Baertsch, R., and Haussler, D. (2008). Using native D.M., Sicks, F., Ludwig, A., Dalén, L., and Hofreiter, M. (2018). Extended and and syntenically mapped cDNA alignments to improve de novo gene finding. continuous decline in effective population size results in low genomic diversity Bioinformatics 24, 637–644. in the world’s rarest hyena species, the brown hyena. Mol. Biol. Evol. 35, Steiner, C.C., and Ryder, O.A. (2011). Molecular phylogeny and evolution of 1225–1237. the Perissodactyla. Zool. J. Linn. Soc. 163, 1289–1303. Westbury, M.V., Petersen, B., Garde, E., Heide-Jørgensen, M.P., and Loren- Suyama, M., Torrents, D., and Bork, P. (2006). PAL2NAL: robust conversion of zen, E.D. (2019). Narwhal genome reveals long-term low genetic diversity protein sequence alignments into the corresponding codon alignments. Nu- despite current large abundance size. iScience 15, 592–599. cleic Acids Res. 34, W609-12. Wilkie, G.S., Davison, A.J., Watson, M., Kerr, K., Sanderson, S., Bouts, T., Svardal, H., Jasinska, A.J., Apetrei, C., Coppola, G., Huang, Y., Schmitt, C.A., Steinbach, F., and Dastjerdi, A. (2013). Complete genome sequences of Jacquelin, B., Ramensky, V., Müller-Trutwin, M., Antonio, M., et al. (2017). elephant endotheliotropic herpesviruses 1A and 1B determined directly from Ancient hybridization and strong adaptation to viruses across African vervet fatal cases. J. Virol. 87, 6700–6712. monkey populations. Nat. Genet. 49, 1705–1713. Teng, H., Zhang, Y., Shi, C., Mao, F., Cai, W., Lu, L., Zhao, F., Sun, Z., and Willerslev, E., Gilbert, M.T.P., Binladen, J., Ho, S.Y., Campos, P.F., Ratan, A., Zhang, J. (2017). Population genomics reveals speciation and introgression Tomsho, L.P., da Fonseca, R.R., Sher, A., Kuznetsova, T.V., et al. (2009). Anal- between brown Norway rats and their sibling species. Mol. Biol. Evol. 34, ysis of complete mitochondrial genomes from extinct and extant rhinoceroses 2214–2228. reveals lack of phylogenetic resolution. BMC Evol. Biol. 9, 95. Thorne, J.L., Kishino, H., and Painter, I.S. (1998). Estimating the rate of evolu- Xue, Y., Prado-Martinez, J., Sudmant, P.H., Narasimhan, V., Ayub, Q., Szpak, tion of the rate of molecular evolution. Mol. Biol. Evol. 15, 1647–1657. M., Frandsen, P., Chen, Y., Yngvadottir, B., Cooper, D.N., et al. (2015). Moun- Tougard, C., Delefosse, T., Hänni, C., and Montgelard, C. (2001). Phylogenetic tain gorilla genomes reveal the impact of long-term population decline and relationships of the five extant Rhinoceros species (Rhinocerotidae, Perisso- inbreeding. Science 348, 242–245. dactyla) based on mitochondrial cytochrome b and 12S rRNA genes. Mol. Phy- Yang, Z. (2007). PAML 4: phylogenetic analysis by maximum likelihood. Mol. logenet. Evol. 19, 34–44. Biol. Evol. 24, 1586–1591. Tunstall, T., Kock, R., Vahala, J., Diekhans, M., Fiddes, I., Armstrong, J., Paten, Yang, D.Y., Eng, B., Waye, J.S., Dudar, J.C., and Saunders, S.R. (1998). Tech- B., Ryder, O.A., and Steiner, C.C. (2018). Evaluating recovery potential of the nical note: improved DNA extraction from ancient bones using silica-based northern white rhinoceros from cryopreserved somatic cells. Genome Res. 28, spin columns. Am. J. Phys. Anthropol. 105, 539–543. 780–788. Van Couvering, J.A., and Delson, E. (2020). African land mammal ages. Yates, A.D., Achuthan, P., Akanni, W., Allen, J., Allen, J., Alvarez-Jarreta, J., J. Vertebr. Paleontol. 40, e1803340. Amode, M.R., Armean, I.M., Azov, A.G., Bennett, R., et al. (2020). Ensembl 2020. Nucleic Acids Res. 48 (D1), D682–D688. Vandenberghe, N., Hilgen, F.J., and Speijer, R.P. (2012). The Paleogene period. In The Geologic Time Scale 2012, F.M. Gradstein, J.G. Ogg, M.D. Yim, H.-S., Cho, Y.S., Guang, X., Kang, S.G., Jeong, J.-Y., Cha, S.-S., Oh, H.- Schmitz, and G.M. Ogg, eds. (Elsevier Science), pp. 855–921. M., Lee, J.-H., Yang, E.C., Kwon, K.K., et al. (2014). Minke whale genome and Veeramah, K.R., Woerner, A.E., Johnstone, L., Gut, I., Gut, M., Marques- aquatic adaptation in cetaceans. Nat. Genet. 46, 88–92. Bonet, T., Carbone, L., Wall, J.D., and Hammer, M.F. (2015). Examining phylo- Yoo, D., Kim, K., Kim, H., Cho, S., Kim, J.N., Lim, D., Choi, S.-G., Choi, B.-H., genetic relationships among gibbon genera using whole genome sequence and Kim, H. (2017). The genetic origin of short tail in endangered Korean dog, data using an approximate bayesian computation approach. Genetics 200, DongGyeongi. Sci. Rep. 7, 10048. 295–308. Zhang, W., Fan, Z., Han, E., Hou, R., Zhang, L., Galaverni, M., Huang, J., Liu, von Seth, J., Dussex, N., Dı́ez-Del-Molino, D., van der Valk, T., Kutschera, V.E., H., Silva, P., Li, P., et al. (2014). Hypoxia adaptations in the grey wolf (Canis Kierczak, M., Steiner, C.C., Liu, S., Gilbert, M.T.P., Sinding, M.S., et al. (2021). lupus chanco) from Qinghai-Tibet Plateau. PLoS Genet. 10, e1004466. Genomic insights into the conservation status of the world’s last remaining Su- matran rhinoceros populations. Nat. Commun. 12, 2393. Zheng, Y., and Janke, A. (2018). Gene flow analysis method, the D-statistic, is robust in a wide parameter space. BMC Bioinformatics 19, 10. Wang, K., Lenstra, J.A., Liu, L., Hu, Q., Ma, T., Qiu, Q., and Liu, J. (2018). Incomplete lineage sorting rather than hybridization explains the inconsistent Zoonomia Consortium (2020). A comparative genomics multitool for scientific phylogeny of the wisent. Commun. Biol. 1, 169. discovery and conservation. Nature 587, 240–245. 12 Cell 184, 1–12, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Biological samples Javan rhinoceros This study NHM, UiO 734; NCBI:txid102233 Merck’s rhinoceros This study F-4160; NCBI:txid2003782 Siberian Unicorn This study IPAE 915/2804; NCBI:txid2491732 Woolly rhinoceros This study ND036; NCBI:txid222863 Greater one-horned rhinoceros This study KB14498/SB137; NCBI:txid9809 Black rhinoceros Moodley et al., 2020 GCA_013634535.1; NCBI:txid9805 Sumatran rhinoceros Lord et al., 2020 GCA_014189135.1; NCBI:txid89632 White rhinoceros SAMN00778988 GCF_000283155.1; NCBI:txid73337 Deposited data Raw sequence data and assemblies This study NCBI Project Number PRJNA687817 Custom scripts for genome alignment This study https://github.com/liushanlin/rhinoceros- analysis comparative-genome Software and algorithms AdapterRemoval v2.2.2 Schubert et al., 2016 https://github.com/MikkelSchubert/ adapterremoval ALLPATHS-LG v.52485 Gnerre et al., 2011 http://software.broadinstitute.org/allpaths-lg/ blog/?page_id=12 ABySS v.2.0.2 Simpson et al., 2009 https://github.com/bcgsc/abyss SOAPdenovo2 Luo et al., 2015 https://github.com/aquaskyline/SOAPdenovo2 BUSCO v.5.0.0 Seppey et al., 2019 https://busco.ezlab.org/ RepeatMasker Open-4.0.2013-2015 Arian Smit & Robert Hubley http://www.repeatmasker.org/ RepeatRunner Smith et al., 2007 https://www.yandell-lab.org/software/ repeatrunner.html Augustus v3.2.3 Stanke et al., 2008 http://bioinf.uni-greifswald.de/augustus/ SNAP Korf, 2004 https://github.com/KorfLab/SNAP Maker v2.31.10 Holt and Yandell, 2011 https://www.yandell-lab.org/software/maker.html BWA v0.7.17-r1188 Li and Durbin, 2009 https://github.com/lh3/bwa LAST v1061 Kie1basa et al., 2011 https://gitlab.com/mcfrith/last MULTIZ v11.2 Blanchette et al., 2004 http://www.bx.psu.edu/miller_lab/ ANGSD v0.924) Korneliussen et al., 2014 http://www.popgen.dk/angsd/index.php/ANGSD MAFFT v7.310 Katoh and Standley, 2013 https://mafft.cbrc.jp/alignment/software/ linuxportable.html UCSC tools Kuhn et al., 2013 https://github.com/ucscGenomeBrowser/kent RaxML v8.2.10 Stamatakis, 2014 https://github.com/stamatak/standard-RAxML ASTRAL v5.6.3 Mirarab et al., 2014 https://github.com/smirarab/ASTRAL DiscoVista v1.0 Sayyari et al., 2018 https://github.com/esayyari/DiscoVista NW_utiles v1.5.0 Junier and Zdobnov, 2010 https://github.com/tjunier/newick_utils/wiki EXONERATE v2.2.0 Slater and Birney, 2005 https://www.ebi.ac.uk/about/ vertebrate-genomics/software/exonerate PAL2NAL v14 Suyama et al., 2006 https://github.com/drostlab/orthologr/tree/ master/inst/pal2nal/pal2nal.v14 MCMCTree v4.8 Yang, 2007 http://abacus.gene.ucl.ac.uk/software/paml. html#download (Continued on next page) Cell 184, 1–12.e1–e7, September 16, 2021 e1 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article Continued REAGENT or RESOURCE SOURCE IDENTIFIER PSMC Li and Durbin, 2011 https://github.com/lh3/psmc SnpEff Cingolani et al., 2012 https://github.com/pcingola/SnpEff Bcftools v1.8 Li, 2011 https://samtools.github.io/bcftools/ GATK v3.8 Auwera et al., 2013 https://gatk.broadinstitute.org/hc/ IQ-Tree v1.6 Nguyen et al., 2015 https://github.com/Cibiv/IQ-TREE RESOURCE AVAILABILITY Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact: Shanlin Liu (shanlin.liu@cau.edu.cn). Materials availability This study did not generate any new unique reagents. Data and code availability Raw sequencing data and genome assemblies can be accessed at NCBI databases under project number Bio- Project:PRJNA687817. The custom scripts for genome alignment analysis for the ancient and historical rhinoceros have been depos- ited in https://github.com/liushanlin/rhinoceros-comparative-genome. EXPERIMENTAL MODEL AND SUBJECT DETAILS Source organisms The specimen of the greater one-horned [Indian] rhinoceros (Rhinoceros unicornis) derives from cell culture held at the San Diego Zoo (ID = KB14498, SB137) and originally derives from a captive born female individual, that in turn was derived from captive born parents (Dam-100288 and Sire-100289). The specimen of black rhinoceros (Diceros bicornis) was provided by the Zululand Rhino Reserve, South Africa (ID = 46373), and is detailed further in the original paper reporting its genome (Moodley et al., 2020). The Sumatran rhi- noceros (Dicerorhinus sumatrensis) sample derives from blood sampled from an individual named Kertam by a staff at the Sabah Wildlife Department in Borneo, Malaysia. Further details on this specimen are detailed in the original paper that released its genome (Lord et al., 2020). The Javan rhinoceros (Rhinoceros sondaicus) genome derives from three sub-samples of dried soft tissue (two) and bone (one) taken from a skull collected in 1838 in Java, that is currently kept in the collections of the Natural History Museum at the University of Oslo (Natural History Museum, Oso; accession Museum id: 734). The woolly rhinoceros (Coelodonta antiquitatis) specimen (ID = ND036; a femur) was collected along the Rakvachan River (N69 17,80’ E167 38,53’), on the Kyttyk Peninsula, Chu- kotka, Russia. The sample has been radiocarbon dated twice, yielding radiocarbon age estimates of 46200 ± 2300 14C years BP (OxA-36569) and 51980 ± 4900 14C years BP (MAG-2095). The Siberian unicorn (Elasmotherium sibiricum) genome was generated from a radius subsample of a specimen (IPAE 915/2804) that originated from Tobolsk, Russia (58 N 68 E) and has a radiocarbon age of > 49,200 14C years BP (OxA-34900). The genome of the Merck’s rhinoceros (Stephanorhinus kirchbergensis) derives from a first molar (M1) tooth root that was subsampled from a complete cranium (National Alliance of Shidlovskiy ‘‘Ice Age,’’ Ice Age Museum, Moscow; accession F-4160) recovered from the Chondon River valley in Yakutia, Russia (N70 120 E137 ) and dated to between 48,000 and 70,000 years BP (Kirillova et al., 2017; Cappellini et al., 2019). Lastly the white rhinoceros (Ceratotherium simum simum) genome was taken from the Broad Institute assembly that is publicly available at https://www.ncbi.nlm.nih.gov/assembly/ GCF_000283155.1/. METHOD DETAILS DNA isolation and sequencing DNA of the Javan rhinoceros was extracted and prepared for sequencing in the ancient DNA laboratories at the GLOBE Institute, University of Copenhagen, following standard clean lab procedures (Cappellini et al., 2012; Orlando et al., 2011). DNA was extracted from dried soft tissue using a SDS - DTT - proteinase k buffer (Gilbert et al., 2007) and from bone using a EDTA - urea - proteinase k buffer (Ersmark et al., 2015). Libraries were constructed following the BEST protocol (Caroe et al., 2017) with the modifications as in Mak et al. (2017), and sequenced on the BGISeq500 platform at BGI Shenzhen, with a sequencing strategy of 50 PE. e2 Cell 184, 1–12.e1–e7, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS DNA was extracted from the woolly rhinoceros specimen at the ancient DNA laboratory at the Swedish Museum of Natural History from 50 mg of bone powder collected using a Dremel drill. Then, we extracted DNA using a modified version of protocol C in Yang et al. (1998) as described in Brace et al. (2012). Sequencing libraries were prepared following the BEST library build protocol (Caroe et al., 2017), where the adaptor oligos were custom-designed for the BGISeq 500 Sequencing Platform (Mak et al., 2017). Libraries were amplified in a 50 mL reaction containing 5 U AmpliTaq Gold polymerase (Applied Biosystems, Foster City, CA), 13 AmpliTaqGold buffer, 2.5 mM MgCl2, 0.4 mg/mL bovine serum albumin (BSA), 0.2 mM each dNTP, 0.2 mM BGI forward primer (Mak et al., 2017), 0.2 mM BGI reverse index-primer (Mak et al., 2017) and 10 mL of library DNA template. Amplified libraries were sequenced on a BGISeq 500 platform at BGI Shenzhen with a SE100 sequencing strategy. We extracted DNA from the Merck’s rhinoceros M1 tooth root, as previously reported by Kirillova et al. (2017) and generated data as described in Cappellini et al. (2019). Briefly, in addition to the DNA extract obtained in Kirillova et al., a second DNA extraction was performed following the previously reported methods (Dabney et al., 2013) in the University of California Santa Cruz Paleogenomics laboratory. We constructed four libraries (single-indexed and double-stranded) for each DNA extract following the methods reported by Meyer and Kircher (2010). The total of eight indexed libraries were pooled and sent to SciLifeLab (Stockholm, Sweden) for sequencing on two lanes of the Illumina HiSeq-X platform (150 PE). In addition, a subsample was processed in the ancient DNA lab- oratories at the GLOBE Institute, University of Copenhagen, following the same protocols as described with the same library con- struction method for the Javan rhinoceros, and sequenced on three lanes of PE100 data on a BGISeq 500 platform at BGI Shenzhen. A subsample was taken from a radius (IPAE 915/2804) of the Siberian unicorn and DNA was extracted at the Australian Centre for Ancient DNA (University of Adelaide) as previously reported by Kosintsev et al. (2019). A new Illumina sequencing library was created for this study following the protocol described by Meyer and Kircher (2010) with an additional step added to excise deaminated cy- tosines with the USER enzyme mix (New England Biolabs) as described by Briggs et al., (2010). The library was split into 8 separate PCR reactions to minimize PCR bias and maintain library complexity. Each PCR of 25 mL contained 1 3 HiFi buffer, 2.5 mM MgSO4, 1 mM dNTPs, 0.5 mM each primer (containing a unique combination of 7-mer i5 and i7 indexes), 0.1 U Platinum Taq Hi-Fi polymerase and 3 mL DNA. The cycling conditions were 94 C for 6 min; 7 cycles of 94 C for 30 s, 60 C for 30 s, and 68 C for 40 s; followed by 68 C for 10 min. Following PCR, replicates were pooled and purified using 1.1x volume AxyPrep magnetic beads, eluted in 30 mL EB buffer, and quantified using a TapeStation (Agilent Technologies). The library was sequenced on three lanes of an Illumina HiSeq X Ten at the Garvan Institute of Medical Research (Sydney, Australia). 4 additional PCR reactions were conducted in Stockholm, Sweden. Each PCR of 25 mL contained 1x AccuPrime reaction mix, 0.3 mM IS4 amplification primer, 0.3 mM P7 indexing primer, 7 U AccuPrime Pfx (Thermo Scientific) polymerase. The cycling conditions were 95 C for 2 min, 14 cycles at 95 C for 15 s, 60 C for 30 s and 68 C for 1 min. PCR replicates were pooled and purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA), eluted in 36 mL EB Buffer, and quantified using a high-sensitivity DNA chip on a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). The library was subsequently sequenced on one Illumina HiSeqX lane with a 2 3 151 bp setup in the High Output mode at SciLifeLab (Stockholm, Sweden). For the historic/ancient rhinoceros samples, although multiple sequencing libraries were constructed and some were sequenced with different read length strategies, comparison tests of the different datasets found no bias differences (detailed in the Additional resources section) and thus we merged them for the subsequent analyses. High molecular weight DNA was extracted from a cell culture specimen from the San Diego Zoo (ID = KB14498, SB137) for the greater one-horned rhinoceros using a Kingfisher duo prime extraction robot. Paired-end Truseq PCR-free libraries of insert size of 180bp and 670bp were constructed and sequenced on the Illumina HiSeq X platform at SciLifeLab (Stockholm, Sweden) generating approximately 400 million paired-end (2x150bp) each. Additionally, three mate-pair libraries (3, 5 and 20 kb) from two specimens, one from San Diego Zoo (KB17733) and one from Rotterdam Zoo (IR_104724), were constructed and sequenced for approximately 100 million paired-end reads (2x150bp) each on the Illumina HiSeq X at SciLifeLab (Stockholm, Sweden). Genome assembly and annotation We merged the paired-end reads for all the ancient and historic rhino samples (Siberian unicorn, Javan, Merck’s and woolly rhinoc- eros) using a modified version of AdapterRemoval (Schubert et al., 2016) which masked the conflict bases that have identical sequencing qualities to Ns and removed collapsed reads < 30 bp. Then, we mapped the short reads from those rhinoceros species onto their corresponding genome references (Table S2) and obtained their genomes (consensus sequences) using the doFasta mod- ule (-doFasta 2) in ANGSD (version 0.924) (Korneliussen et al., 2014) with a minimum requirement of mapping quality and base quality of 20. A minimum depth of 5 was set for the Javan and woolly rhinoceros, and the minimum depth was set as 3 for the Siberian unicorn and Merck’s rhinoceros due to their relatively lower coverages (303 for the former pair, versus 103 for the latter pair). For the greater one-horned rhinoceros, we obtained its genome using a mate-pair assembly. Since assembly quality can vary for different datasets, we used three different assemblers and evaluated their performance. The following assemblers for short read sequence data were used: ALLPATHS-LG v.52485 (HAPLOIDIFY = True) (Gnerre et al., 2011), ABySS v.2.0.2 (-k = 61) (Simpson et al., 2009) and SOAPdenovo2 v. 2.04 (-K 61) (Luo et al., 2015). Out of the three assemblers, ALLPATHS-LG was selected for down- stream analyses as it produced the most gene-complete and the contiguous assembly (with a scaffold N50 of 27.7 Mbp. Gene completeness was measured with BUSCO v.5.0.0 (Seppey et al., 2019) using the ‘‘mammalia_odb10’’ ortholog dataset, which showed a low degree of missing, fragmented and duplicated genes: ‘‘C:96.2% [S:95.6%, D:0.6%], F:1.0%, M:2.8%, n:9226.’’ Cell 184, 1–12.e1–e7, September 16, 2021 e3 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article The genome annotation of white rhinoceros is its original release with its genome assembly. Since genome annotations of both the black and Sumatran rhinoceros were not released with their genome assemblies, we annotated the genomes of the black, Sumatran and greater one-horned rhinoceroses as follows. First, repeats were masked using RepeatMasker (Smit, AFA, Hubley, R & Green, P. RepeatMasker Open-4.0.2013-2015) with all (‘‘model_org=all’’) species being included in the RepBase. Then, we masked the trans- posable element proteins using RepeatRunner (Smith et al., 2007) and the repeat protein library te_proteins.fasta (downloaded from http://weatherby.genetics.utah.edu/data/maker_tutorial.tgz). After that, we applied Ensemble (version 64) (Yates et al., 2020) for homolog prediction with its amino acid datasets including Sus scrofa, Tursiops truncatus, Bos taurus, Drosophila melanogaster, Vicugna pacos and Homo sapiens. De novo gene prediction was achieved with: a) Augustus (version 3.2.3) (Stanke et al., 2008), where ‘‘human’’ was used as the gene prediction species model for Augustus; b) and SNAP (Korf, 2004) using the provided mammal model (mam54). Finally, Maker (version 2.31.10) (Holt and Yandell, 2011) was used to weigh and merge different evidence, and obtain protein-coding-gene sets for the black, Sumatran and greater one-horned rhinoceros. Genome alignment and phylogenetic inference The four de novo assembled rhinoceros genomes and the genome of tapir (Tapirus terrestris, GCA_004025025.1) were aligned to the horse reference genome (Equus caballus, GCF_000002305.2) using LAST with parameters of -m 100 -E 0.05 (Kie1basa et al., 2011), and each pairwise alignment was chained and netted to form high quality blocks. After that, the non-syntenic regions were filtered out from each alignment. Finally, the pairwise local alignments were combined to generate a final multi-species whole genome alignment (MWGA) using MULTIZ without reference fixed and with a radius of 50 in dynamic programming (Blanchette et al., 2004). After generating the MWGA for all the de novo assemblies, a sliding window-based phylogenetic analysis was conducted along the horse genome with a window size of 100 kb. A set of UCSC genome browser tools (Kuhn et al., 2013), including mafsInRegion, maf- Filter, maf2fasta, and several in-house PERL scripts were applied to divide the whole genome alignment into sub-regions. Then, for each sub-region, we included the genome sequences of the non-de novo sequenced rhinoceros samples, via extracting the corre- sponding sequences based on the region coordinate information on their reference genomes, and conducted multiple sequence alignment using MAFFT (Katoh and Standley, 2013). We also removed the alignments with effective length < 1,000 bp and effective ratio < 0.5, of which the effective sites represent nucleotide sites that do not include missing information (Ns) for any taxa in the mul- tiple sequence alignment. After that, the GTR+CAT model of RaxML (Stamatakis, 2014) was used to build Maximum Likelihood (ML) phylogenetic trees for each window. Finally, the species tree was generated using the multi-species coalescent model-based soft- ware ASTRAL III (Mirarab et al., 2014), by combining all the regional trees. Due to the poor quality of the endogenous DNA recovered from the historical and ancient samples, their genome sequences could only be reconstructed through mapping against a reference genome. To test whether the choice of reference genome would play any role in shaping either the consensus sequence recovered, or phylogenetic inference based on the data, we used the 100 PE reads generated from the woolly rhinoceros and the Siberian unicorn to examine reference bias via aligning their short reads onto different references including the horse (Equus caballus), white, greater one-horned, and Sumatran rhinoceros, respectively. Following genome alignment, 82 genome regions of length > 20,000 bp were selected to infer phylogenetic trees (Figure S1A). Those regions were randomly scattered across the genome, thus having little chance of being found in the same recombination blocks. We calculated the frequency of three topologies around each focal internal branch of the species tree (Figure 2B) using DiscoVista (Sayyari et al., 2018) and NW_utils (Junier and Zdobnov, 2010), for the 100 kb windows based gene trees. A multi-species coalescent simulation was also applied to determine the expected gene tree distributions on the basis of the dated species tree using the Phy- base package following Wang et al. (2018). Then, we inspected the congruency between the frequency of the three topologies in- ferred from empirical genomic data and that generated from the simulation for the Rhinocerotinae lineage. To validate the robustness of our species tree, we also reconstructed the phylogenetic relationship across each chromosome inde- pendently (reference by horse’s chromosomes) for all the rhinoceros species that have de novo genome assemblies and two out- groups of tapir and horse using the aforementioned 100 kb sliding windows based method. Then, for each chromosome we inferred a species tree and calculated the tree topology frequency using DiscoVista (Sayyari et al., 2018). As a genomic region with length of 100 kb may contain multiple recombination breakpoints, we sampled a short alignment with length of 5,000 bp within each 100 kb sliding window to infer gene trees using RaxML (Stamatakis, 2014) with a substitution model setting of GTR+CAT and 100 bootstrap replicates. Then we filtered out gene trees with nodes of bootstrap support < 85 to guarantee congruent signals throughout each subregion. Finally, the frequency of three topologies around the phylogenetic discordance branch obversed based on 100 kb sliding windows were calculated using DiscoVista (Sayyari et al., 2018). Molecular dating Orthologs of the species with de novo genome assemblies were extracted from the former synteny checked MWGA: 1) gene location and the corresponding exon regions in the horse genome annotation file were used to extract CDSs from the MWGA alignment asking for a coverage R 80%; 2) the amino acid sequences obtained from the horse genome serving as query were used to find correspond- ing homologs for each species using protein2genome model in EXONERATE (Slater and Birney, 2005); 3) exons that were not shared by all species were removed to improve the gene alignment accuracy; 4) amino acid sequences were aligned using MAFFT (Katoh and Standley, 2013) and in turn used to guide the CDS alignment using PAL2NAL (Suyama et al., 2006). For the samples without de novo assemblies, we obtained their orthologs according to the gene and exon location information from their reference genomes. In e4 Cell 184, 1–12.e1–e7, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS order to further diminish the potential influence of reference bias, we obtained two gene sequences for each ortholog by mapping their reads onto two different reference genomes (Table S2), and merged the two gene sequences while masking conflicts as ‘‘N’’ (X for amino acids). We inferred the evolutionary timescale of rhinoceros lineages using a set of 3,513 orthologs from the set identified in the previous section (8,820,642 nucleotides). The set of orthologs was selected to minimize the possible bias in molecular rates arising from an- alyses of finite number of sites or excessive discrepancies in molecular rates across lineages in particular loci (Marshall et al., 2006). To select orthologs, we first performed a fast tree search for each ortholog, using only NNI tree rearrangements, under a GTR+R4 substitution model (Kalyaanamoorthy et al., 2017) implemented in IQ-TREE v1.6 (Nguyen et al., 2015). Orthologs were considered as adequate for molecular dating analyses if they contained low degrees of rate variation among lineages (Doyle et al., 2015), and if their gene trees were not too greatly distinct from the species tree (Mendes and Hahn 2016). Specifically, the orthologs retained had gene trees with coefficients of variation in root-to-tip length (non-clocklikeness) < 0.1, and Robinson-Foulds distances to our inferred species tree % 2. Our analyses included time-calibrations on four splits of the species tree, each based on several lines of evidence from the fossil record. The nodes were calibrated using uniform distributions with soft maximum bounds, placing a 0.025 probability on older ages. Placement of maximum bound for a clade assumes that less inclusive clades cannot be older than the oldest fossils of more inclusive clades. (i) We calibrated the split between Elasmotheriinae and Rhinocerotinae as occurring between 35 Ma and 44 Ma. This minimum bound is supported by the earliest record of Rhinocerotinae, with Epiaceratherium naduongense, at between 35-39 Ma (Böhme et al., 2013), which is unambiguously a species nested within Rhinocerotinae. Additional evidence for this minimum bound is the earliest species assigned to Elasmotheriinae, Penetrigonas dakotensis, which appears at ca. 38 Ma (Heissig, 2012). Furthermore, morpho-anatomical cladistic analysis has placed the 38 My-old fossil Subhyracodon occidentalis as a representative of Elasmother- iinae (Becker et al., 2013). The maximum bound is supported by fossil and phylogenetic evidence placing with high confidence the minimum age of Rhinocerotidae at 44 Ma (Bai et al., 2020). (ii) The split between Rhinoceros and Dicerorhinus was calibrated as occur- ring between 13 Ma and 23 Ma. The minimum age is based on remains from the middle Siwaliks of Pakistan dating from ca. 13 Ma and unambiguously assigned to Dicerorhinus, under the name of D. aff. sumatrensis (Heissig, 1972). Other remains from the same region and age confirm the occurrence of Dicerorhinus as early as the middle Miocene (13.7-11.65 Ma; P.-O.A., unpublished data). The maximum bound is informed by the first attested occurrence of Rhinocerotina in the fossil record (one-horned rhinocerotine Gainda- therium cf. browni; Bugti Hills, Pakistan: 22.6 Ma; Antoine et al., 2013). (iii) The split between Rhinoceros unicornis and Rhinoceros sondaicus was calibrated as occurring between 1.9 Ma and 5.3 Ma. The minimum bound is supported by the earliest occurrence of both species in the fossil record (Antoine et al., 2021). A morphology-based phylogenetic analysis of Rhinocerotina, with a compre- hensive sample within Rhinoceros, retrieves the following topology: (R. sondaicus,(R. sinensis,(R. unicornis,(R. kendengindicus, R. platyrhinus)))), with the first occurrence of Rhinoceros platyrhinus estimated at around the Pliocene–Pleistocene transition (2.58 Ma; Antoine et al., 2021). The maximum age constraint coincides with the earliest estimated age of Rhinoceros from the same morphology-based phylogenetic analyses. (iv) The split between Ceratotherium simum and Diceros bicornis was calibrated to occur between 5.3 Ma and 7.3 Ma. This range is consistent with records of Diceros bicornis recognized in upper Miocene deposits (> 5.3 Ma) at Lothagam (Kenya, 6.54-5.2 Ma; Brown and McDougall, 2011) and Albertine (Uganda, 7.25-5.3 Ma; Pickford et al., 1993). Molecular data tend to identify the timing of the early radiation of mammals as being older than fossil evidence (e.g., O’Leary et al., 2013; dos Reis et al., 2014). For this reason, we performed analyses that included a fifth node-calibration on the age of the crown of all extant Perissodactyla. This calibration was a soft maximum bound at 66 Ma with a 0.01 prior probability of an older age. The age of this bound is based on the absence of unambiguous crown placental mammals before this time. All analyses were repeated after excluding this calibration for comparison (Figure S5A). Additional molecular dating analyses were performed by including data from the tapir genome (Figure S5A), therefore including representatives of all living families of Perissodactyla (Equidae, Tapiridae, and Rhinocerotidae). The addition of these data was fol- lowed by filtering loci for biased rate estimates as described above, which led to a dataset of similar size (3,163 orthologs with 7,325,631 nucleotides). The species tree topology in these augmented analyses included the tapir as the sister to Rhinocerotidae, while the horse remained as the outgroup. We included a sixth calibration in addition to those described above, this time on the timing of the split between Rhinocerotidae and Tapiridae. The calibration had a hard minimum bound and soft maximum bound (0.025 prior probability on older ages), ranging between 54 Ma and 64 Ma. The minimum age is based on the earliest appearance of Tapiroidea in the fossil record shortly after the Paleocene-Eocene transition, around 54 Ma, with Vastanolophus holbrooki and Cambaylophus vast- anensis in Asia (Rose et al., 2014) and 53.5 Ma with Heptodon in North America (Vandenberghe et al., 2012). The maximum age is based on the earliest record of Notoungulata at 65 Ma (Tiupampa, Bolivia; Muizon and Ladevèze, 2020). The extinct South Amer- ican orders Notoungulata and Litopterna form a clade, sister to Perissodactyla within Panperissodactyla (Buckley, 2015; Welker et al., 2015; Westbury et al., 2017), therefore suggesting a similar, or earlier, age for stem Perissodactyla. Bayesian dating analyses were performed using a GTR+G substitution model and an uncorrelated -gamma relaxed clock model as implemented in MCMCtree, part of PAML v4.8 (Yang 2007). We further addressed heterogeneity in molecular evolutionary processes by partitioning the molecular clock and substitution models into each of the three codon positions (three partition subsets). We improved the efficiency of the analysis using approximate Bayesian computation (Thorne et al., 1998). The posterior distribution was estimated using Markov chain Monte Carlo (MCMC) samples. After a burn-in phase of 105 MCMC steps, samples were drawn Cell 184, 1–12.e1–e7, September 16, 2021 e5 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll OPEN ACCESS Article every 103 MCMC steps over a total of 107 steps. We verified convergence to the stationary distribution by comparing parameter es- timates from four independent runs. Effective sample sizes were verified to be above 200 for all estimated parameters. Alignment tool selection The BWA-ALN algorithm (Li and Durbin, 2009) is commonly used to align such short reads generated in aDNA studies to reference genomes, owing to its stable performance for ultra-short length (typically shorter than 100 bp). However, simulation data showed read mapping success rate of the relatively longer reads (60 bp and 80 bp) dropped about three fold for high divergence reads (3% divergence level) using BWA-ALN (Parks and Lambert, 2015). We therefore compared the performance of BWA-ALN against BWA-MEM (Li, 2013), on the sequencing data obtained for the four historical/ancient samples (Figure S5B). For BWA-ALN, we disabled read seeding (-l 1024) to enhance error tolerance, while we marked shorter split hits as secondary alignments and removed them from the outputs for the BWA-MEM method. Gene flow To explore for potential gene flow between the different rhinoceros species, we computed D-statistics that assesses genetic affinities between taxa based on patterns of shared derived alleles (Green et al., 2010). For a given tree of topology (((H1, H2), H3), O), where O represents an outgroup, under the null hypothesis of no gene flow, shared derived alleles between H1 and H3, or H2 and H3, can only derive from incomplete lineage sorting (ILS) and are expected to be symmetric between the two pairs. Therefore, ancient gene flow is assumed between a pair in which imbalanced aggregation of shared derived alleles are detected. However, we should be aware of the limitations of D-statistics analysis between highly divergent lineages - certain assumptions could be violated in cases such as repeated or independent mutations (Prüfer et al., 2012) and differences in lineage specific mutation rates (Zheng and Janke, 2018). We independently mapped all raw reads from each rhinoceros species and the horse to three different reference genomes (greater one-horned rhinoceros, white rhinoceros and Sumatran rhinoceros) and then calculated D-statistics using the doAbbababa function module in ANGSD (Korneliussen et al., 2014) specifying the horse as the outgroup ancestral allele and the following filtering parameters: 1) use only sites where the reads of the outgroup all have the same base; 2) minimum base quality and mapping quality of 20; 3) use only the transversion sites, all individuals have a minimum depth of 3, and maximum depth of 70; 4) a block size of 5 Mb to estimate standard errors using jackknife procedure. We further assessed whether less efficient mapping of the horse outgroup to the rhinoceros reference genomes may be driving some of our conflicting D-statistics results that differed based on the reference genome choice. We did this by calculating the pair- wise distance of all rhinoceros species to the outgroup horse three times using all individuals mapped to the three different reference genomes. We calculated pairwise distances in 1 Mb sliding windows across all scaffolds > 1 Mb in size using a consensus base call approach (-doIBS 2) in ANGSD. We additionally applied the following filters; only include sites that have coverage in all individuals, a minimum base quality and mapping quality of 20, and only consider transversion differences. Heterozygosity estimation We estimated genome-wide heterozygosity for all the rhinoceros species based on a Site Frequency Spectrum (SFS), of which one diploid individual will generate three allelic states: homozygous ancestral allele state (AA), heterozygotic state (AB) and homozygous derived allele state (BB). Therefore, the whole genome level heterozygosity rate can be calculated as AB/(AA+AB+BB). We applied the doSaf function (-doSaf 1) module in ANGSD (version 0.924) (Korneliussen et al., 2014) to calculate genome-wide heterozygosity for each rhinoceros species with parameters of: unique mapping, a minimum mapping quality of 20, a minimum base quality of 20, a minimum and maximum depth value of 1/3 and 2 times of the average genome depths, respectively, for each species. For the ancient and historic samples (Javan, Merck’s, woolly rhinoceros and Siberian unicorn), transition sites were removed (-noTrans 1) to eliminate the potential influences of DNA damage derived from cytosine deamination. Demographic inference and Runs of Homozygosity (ROH) estimation We applied the Pairwise Sequentially Markovian Coalescent (PSMC) method (Li and Durbin, 2011) to infer the history of population size changes in the ancestors of all the eight rhino species. For the species for which genomes were de novo assembled, we aligned 303 genome coverage of the shotgun reads to recover the diploid genome heterozygosity information with depth filters of R 1/3 and % 2 3 of the average depths. Then, PSMC was used to estimate the distribution of the time to the most recent common ancestor (TMRCA) between the two alleles across all chromosomes using the density information of heterozygous sites across each diploid genome. Changes in effective population size (Ne) were inferred assuming a substitution rate (m) of 2.2 3 108 substitutions/site/gen- eration and a generation time (g) of 25 years for the African rhinoceroses (Figure 3A; Moodley et al., 2020) and a m of 2.34 3 108 and g of 12 years for the other rhinoceros species (Figures 3B and 3C; Dinerstein and McCracken, 1990; Lord et al., 2020; Mays et al., 2018). The consistency of the demographic results was tested by performing 100 bootstrap replicates as shown in Figure 4. ROH segments were recorded by summarizing genome regions of which the TMRCA dated to a recent time period (best fit K value of % 2) following the method in Palkopoulou et al. (2015). e6 Cell 184, 1–12.e1–e7, September 16, 2021 Please cite this article in press as: Liu et al., Ancient and modern genomes unravel the evolutionary history of the rhinoceros family, Cell (2021), https://doi.org/10.1016/j.cell.2021.07.032 ll Article OPEN ACCESS Genetic load and identification for rhinoceros specific frameshift mutation First, gene annotation files of the de novo sequenced genomes (white, black, Sumatran and greater one-horned rhinoceroses) were imported into SnpEFF (Cingolani et al., 2012) to create a gene function database. Then, we obtained the heterozygous sites for each rhinoceros species using BCFTOOLS (Li, 2011) requiring: minimum and maximum depth of 1/3 and 2x of the average depth; mapping and base quality R 20; variance quality R 30; read supporting number of minor allele account for R 20% of the total depth. After that, functional classes of the variances were estimated using SnpEFF (Cingolani et al., 2012) with default settings. We trimmed 10 bp from the 50 and 30 ends of each read for the Merck’s rhinoceros and Siberian unicorn to alleviate the impact of DNA damages, because both the samples showed abnormal transition/transversion ratios - as high as 5.367 without end trimming (Figure S4B). However, as the transition/transversion ratio of the Merck’s rhinoceros stayed abnormally high after trimming, which is consistent with previous re- ports of its DNA damage pattern (Figure S4 in Cappellini et al., 2019), we excluded it from further comparisons To compare the rate of missense and loss-of-function mutations in rhinoceros to that of other mammals we obtained published re- sequencing data for 402 mammalian genomes from 30 species and mapped these to the phylogenetically closest available reference genome for each species (detailed in Table S4) using BWA-MEM v0.7.17 (Li, 2013). We then obtained and filtered variant calls for each individual using GATK HaplotypeCaller v3.8 following the ‘‘short variant discovery best-practices guidelines including ‘‘hard filtering’’ (Auwera et al., 2013). Additionally, we only retained within-species bi-allelic sites, and removed all indels and sites found at a frequency that was either less than one third, or greater than three times, the genome-wide autosomal coverage (Li, 2014). To investigate the putative genetic background of unique rhinoceros biological adaptations, we explored frameshift mutations shared across the rhinoceros lineages. We inspected the frameshift mutations in the alignment results generated by the EXONERATE analysis mentioned above. First, all the frameshift mutations with their location information (Gene IDs and Exon IDs) were extracted out for rhinoceros species with de novo assemblies being available (black, white, Sumatran and greater one-horned rhinoceroses). Then, we filtered out the frameshift mutations that share the same locations and exist in all four rhinoceros species, and recorded their gene information for further examination. Test of heterozygosity, PSMC and ROH estimation for the non-modern rhino samples Since four of the rhinoceros genomes were generated from historic and ancient samples, data could only be generated from them by short read re-sequencing. Furthermore, given they represent historical or ancient samples, they are expected to contain damaged DNA (principally cytosine deamination (Pääbo, 1989)) manifested as elevated C > T/G > A transitions at the 50 and 30 ends of sequence reads (Briggs et al., 2007; Brotherton et al., 2007). Therefore, we tested the viability of restricting a number of the analyses performed to only the transversions, including heterozygosity estimation, historical demographical inference, and ROH estimation (Figures S3 and S4). Also, the influences of the other factors such as reference selection and sequencing depth were evaluated. For the ROH estimation, we applied the ROH inference method to the four rhinoceros species with de novo genome assemblies, using the data- sets both with and without transitions. For the PSMC estimation, we mapped short reads representing various data volumes (different genome coverages, 113, 173 and 353) from the black rhinoceros onto the white rhinoceros genome to test whether it is possible to obtain accurate demographic history results using short read data from species that have a closely-related reference genome available. QUANTIFICATION AND STATISTICAL ANALYSIS Analyses related to genome assembly, assembly quality evaluation, short read alignment and evolutionary relationship inferences can be found in the Method details section. ADDITIONAL RESOURCES The results of some additional tests could be useful in understanding the sequence dataset for the historic/ancient rhinoceros sam- ples and the comparative genome analyses, and can be found at: https://github.com/liushanlin/rhinoceros-comparative-genome/ blob/main/additional_resources.md Cell 184, 1–12.e1–e7, September 16, 2021 e7 ll Article OPEN ACCESS Supplemental figures Figure S1. Phylogenetic relationship inferences, related to Figure 2 (A): Cladogram generated using different rhinoceros genomes as mapping references. Red branches represent different consensus sequences for the woolly rhi- noceros generated when using different reference genomes for the mapping (second species named in each tip label), while blue represent the Siberian unicorn (Elasmotherium sibiricum). The cladogram shows a clear effect of reference bias when short reads are mapped on to a very distantly-related genome in order to obtain the consensus sequences. For example, while the position of the woolly rhinoceros is correct when mapped to either the white, Sumatran or greater one-horned rhinoceros, it is placed erroneously sister to all other rhinoceroses when mapped against the horse genome. Note that relationships were inferred under the multi- species coalescent as implemented in the software ASTRAL, which does not provide branch lengths in units of molecular substitutions nor terminal lengths, and thus we only show a tree topology. Branch labels indicate local posterior probabilities. (B): Tree topologies obtained by mapping the Siberian unicorn against different reference genomes (I: white rhinoceros; II: greater one-horned rhinoceros) and using only the genomes that were de novo assembled (III). Branch labels indicate local posterior probabilities estimated in ASTRAL. (C): The distribution of gene trees across the genome using the horse’s chromosome ID as reference. The upper panel shows chromosomes that placed majority support on the species tree shown in Figure 2, while the bottom panel displays chromosomes with a greater frequency of an alternative gene tree. (D): Comparison of the frequency of the three bipartitions between 100 kb sliding windows and 5,000 bp subregions. Branches 4 and 5 correspond to the tree in panel C in Figure 2. ll OPEN ACCESS Article Figure S2. D-statistic analyses, related to Figure 2 GOR: greater one-horned rhinoceros; SUR: Sumatran rhinoceros; WHR: white rhinoceros; SIU: Siberian unicorn; BLR: black rhinoceros; JAR: Javan rhinoceros; MER: Merck’s rhinoceros; WOR: woolly rhinoceros. Each dot represents a test using the topology and species in the subtitles. For example, (((GOR, SUR), BLR), Horse) represents one of the topologies in Panel A-III and is shown as one of the dots with the x axis of BLR. For a preset topology (((H1,H2),H3),H4), a negative D- statistics value means a closer relationship between H1 and H3 compared to that of H2 and H3, while a positive value means that H2 is closer to H3 than H1 is. Z score reflects the significance of the test with a value > 3 being considered as statistically significant. (A): D-statistics score inference bias generated from reference selection. Species names on the x axis represent different H3 settings in the topology presets. It shows biases in the D-statistics scores that may have arisen as a consequence of the choice of reference genome against which the data was mapped prior to the (legend continued on next page) ll Article OPEN ACCESS analysis. For example, when the species used as the mapping reference is placed at the H1 position, in every case, we observed an excess of shared derived alleles between the H2 and H3 species. We suggest that the deep divergence between H1 and H4 means that only more conserved regions, and regions more similar to the mapping reference will map. This will artificially make the alleles in the reference individual look more similar to the outgroup ancestral allele, resulting in more ABBA patterns and ultimately the biased D-statistics results that we see. (B) showed that regardless of the reference genome selected (i.e GOR, SUM, or WHR), a decreased pairwise distance between the horse and the reference genome (and the sister species of the reference to a lesser extent), relative to the other non-reference rhinoceros species. Taken together, using either the H1 or H2 individual as the mapping reference for this analysis will result in biased and likely incorrect inferences. Therefore, to avoid such obvious biases caused by the lower mapping efficiency of the outgroup horse genome, we only use D-statistics to infer gene flow when the reference genome is either in the H3 position, the sister species to the H3, or absent from the test entirely; (B): The whole genome pairwise genetic distance between horse and the different rhinoceros species, when initially mapped against different rhinoceros reference genomes (indicated as red color, GOR (upper), SUR (middle) and WHR (bottom)). y axis represents the pairwise distance to horse; (C): Excess of genetic affinities between Rhinoceros species (GOR and JAR) and the African rhinoceros (WHR and BLR). Species names on the x axis represent different H2 settings in the topology preset; (D): Excess of genetic affinity between Rhinoceros and the two extinct DCS clade members (MER and WOR), but not with the Sumatran rhinoceros using two mapping references (GOR and WHR). Species names on the x axis represent different H2 settings in the topology preset. ll OPEN ACCESS Article Figure S3. Heterozygosity estimations, related to Figure 3 GOR: greater one-horned rhinoceros; SUR: Sumatran rhinoceros; WHR: white rhinoceros; BLR: black rhinoceros; SIU: Siberian unicorn; JAR: Javan rhinoceros; MER: Merck’s rhinoceros; WOR: woolly rhinoceros. (A): Difference between heterozygosity ratio estimated taking into account all the sites and only the transversion sites. Numbers represent the median values for each group. For each box, the central line represents median value, and box limits represent first and third quartile, and whisker extends from the hinge to their largest values, but < 1.5*IQR (Interquartile range). (B): Whole genome heterozygosity levels in different taxonomic groups. Note that all the humans and rhinoceros were excluded in the current figure. The yellow horizontal line represents a low heterozygosity of 0.15% similar to that of the rhinoceroses. In the bottom panel we only include families in Artiodactyla, Carnivora, (legend continued on next page) ll Article OPEN ACCESS Primates and Rodentia that have multiple species of low heterozygosity level (upper panel), and families that have > 2 species in our dataset (Table S3 for all the details). (C): Factors influencing heterozygosity ratio estimation. Numbers represent the median values for each group. Influence of genome coverage (Depth of 303, 153 and 103 that indicated by different shades of blue) was tested by mapping reads from the black rhinoceros onto the reference genome of white rhinoceros. We noticed that different reference genomes contributed little to the variation of heterozygosity estimation, although when compared to the heterozygosity rate estimated using its own genome as a reference, using other closely related species tends to augment the estimation to a small degree and the inflation worsen when mapped to a higher divergent reference genome. ll OPEN ACCESS Article Figure S4. Tests for ROH, PSMC, and gene effect estimations, related to Figures 3, 4 and 5 (A): ROH estimation comparison between the usage of all sites versus only transversions. This comparison shows that using transversions only leads to the detection of fewer ROH segments in total but tends to enlarge the ROH segments to longer stretches, and the differences are more obvious for the two rhinoceros species with relatively lower whole genome heterozygosity levels - the greater one-horned and white rhinoceroses; (B): Transition/Transversion (Ts/Tv) ratio of different rhinoceros species. Trim-Ten (dark khaki) represent the Ts/Tv ratios calculated based on reads that were trimmed 10 bp from both the 50 and 30 end, while noTrimming (light blue) represent the rates obtained without read trimming; (legend continued on next page) ll Article OPEN ACCESS (C): Demographic history estimated for black rhino. BR refers to black rhinoceros as reference genome, while WR refers to white rhinoceros as reference genome. FN means the false negative rate set for the PSMC curve plotting. The 353, 173 and 113 represent the genome coverages for each test. It shows that adequate coverage (~303) short reads, or inadequate coverage with proper FN settings can generate reliable demographic estimation; (D): Demographic trajectories of different rhinoceros species. Results were estimated using both all heterozygous sites with a substitution rate of 1.95x108 per generation, and transversion heterozygous sites with a substitution rate of 0.65x108 per generation. We observed that the transversion only heterozygosity sites can also deliver comparable demographic results. ll OPEN ACCESS Article Figure S5. Analyses of molecular dating and alignment tool selection, related to Figure 2 and STAR Methods (A): Molecular dating using different datasets and time-calibrations. The x-axes show the times of divergence in units of Million years ago (Mya). GOR: greater one-horned rhinoceros; SUR: Sumatran rhinoceros; WHR: white rhinoceros; SIU: Siberian unicorn; BLR: black rhinoceros; JAR: Javan rhinoceros; MER: Merck’s rhinoceros; WOR: woolly rhinoceros. (legend continued on next page) ll Article OPEN ACCESS (B): Comparison between the bwaALN and bwaMEM algorithms based on aligning one million randomly selected reads against their closely related reference genomes. The taxa used as reference for Siberian unicorn, Javan, woolly and Merck’s rhinoceros are white, greater one-horn and Sumatran rhinoceroses, respectively. ALN and MEM on the x axis represent the number of reads that mapped successfully (Mapping Quality value R 10) using the corresponding tools but failed with the other, and ‘‘random’’ represents the reads that had mapping quality values < 10 using both tools. ‘‘same’’ and ‘‘different’’ represent the number of reads that were mapped onto the same and different locations, respectively, using the two methods.