1. Evolutionary Biology
  2. Genetics and Genomics

Radiocarbon and genomic evidence for the survival of Equus Sussemionus until the late Holocene

  1. Dawei Cai  Is a corresponding author
  2. Siqi Zhu
  3. Mian Gong
  4. Naifan Zhang
  5. Jia Wen
  6. Qiyao Liang
  7. Weilu Sun
  8. Xinyue Shao
  9. Yaqi Guo
  10. Yudong Cai
  11. Zhuqing Zheng
  12. Wei Zhang
  13. Songmei Hu
  14. Xiaoyang Wang
  15. He Tian
  16. Youqian Li
  17. Wei Liu
  18. Miaomiao Yang
  19. Jian Yang
  20. Duo Wu
  21. Ludovic Orlando  Is a corresponding author
  22. Yu Jiang  Is a corresponding author
  1. Bioarchaeology Laboratory, Jilin University, China
  2. Key Laboratory of Animal Genetics, Breeding and Reproduction of Shaanxi Province, College of Animal Science and Technology, Northwest A&F University, China
  3. Heilongjiang Provincial Institute of Cultural Relics and Archaeology, China
  4. Shaanxi Provincial Institute of Archaeology, China
  5. Ningxia Institute of Cultural Relics and Archaeology, China
  6. College of Earth and Environmental Sciences; MOE Key Laboratory of Western China's Environmental Systems, Lanzhou University, China
  7. Centre d’Anthropobiologie et de Génomique de Toulouse (CAGT), CNRS UMR 5288, Université de Toulouse, Université Paul Sabatier, Toulouse, France, France
https://doi.org/10.7554/eLife.73346
Share this article
Cite this article
  1. Dawei Cai
  2. Siqi Zhu
  3. Mian Gong
  4. Naifan Zhang
  5. Jia Wen
  6. Qiyao Liang
  7. Weilu Sun
  8. Xinyue Shao
  9. Yaqi Guo
  10. Yudong Cai
  11. Zhuqing Zheng
  12. Wei Zhang
  13. Songmei Hu
  14. Xiaoyang Wang
  15. He Tian
  16. Youqian Li
  17. Wei Liu
  18. Miaomiao Yang
  19. Jian Yang
  20. Duo Wu
  21. Ludovic Orlando
  22. Yu Jiang
(2022)
Radiocarbon and genomic evidence for the survival of Equus Sussemionus until the late Holocene
eLife 11:e73346.
https://doi.org/10.7554/eLife.73346
  2. Download BibTeX
  3. Download .RIS

Abstract

The exceptionally rich fossil record available for the equid family has provided textbook examples of macroevolutionary changes. Horses, asses, and zebras represent three extant subgenera of Equus lineage, while the Sussemionus subgenus is another remarkable Equus lineage ranging from North America to Ethiopia in the Pleistocene. We sequenced 26 archaeological specimens from Northern China in the Holocene that could be assigned morphologically and genetically to Equus ovodovi, a species representative of Sussemionus. We present the first high-quality complete genome of the Sussemionus lineage, which was sequenced to 13.4× depth of coverage. Radiocarbon dating demonstrates that this lineage survived until ~3500 years ago, despite continued demographic collapse during the Last Glacial Maximum and the great human expansion in East Asia. We also confirmed the Equus phylogenetic tree and found that Sussemionus diverged from the ancestor of non-caballine equids ~2.3–2.7 million years ago and possibly remained affected by secondary gene flow post-divergence. We found that the small genetic diversity, rather than enhanced inbreeding, limited the species’ chances of survival. Our work adds to the growing literature illustrating how ancient DNA can inform on extinction dynamics and the long-term resilience of species surviving in cryptic population pockets.

Editor's evaluation

This article represents multiple milestones in our understanding of the evolution and extinction of Pleistocene equids, including revising the timing of extinction and clarifying the evolutionary history of Equus (Sussemionus) ovodovi. The discovery of the late persistence of non-caballine equid taxa in Northern China until deep into the late Holocene is particularly important. This finding will be of broad interest to the paleontology, paleoecology, archaeology, and paleogenomic communities and should stimulate important future research into equid extinction processes.

https://doi.org/10.7554/eLife.73346.sa0

Introduction

Today, all of the seven extant species forming the horse family belong to one single genus, Equus. It emerged in North America some 4.0–4.5 million years ago (Orlando et al., 2013), and first spread into Eurasia ~2.6 million years ago (Mya), via the Beringia land bridge (Lindsay et al., 1980). This first vicariance and expansion out of America led to the emergence of the ancestors of zebras, hemiones, and donkeys, a group collectively known as non-caballine (or stenonine) equids. Another expansion through Beringia occurred ~0.8–1.0 Mya (Vershinina et al., 2021), which allowed caballine equids (i.e., those most closely related to the horse) to enter into the Old World, where they persisted until the modern era and were domesticated ~5500–4200 years ago (Gaunitz et al., 2018; Librado et al., 2021; Outram et al., 2009).

In the recent years, ancient DNA (aDNA) data have revealed that the genetic diversity of non-caballine Equus was considerably larger in the past than it is today (Librado and Orlando, 2021; Orlando, 2020). This was further confirmed as the first mitochondrial DNA (mtDNA) data of Equus (Sussemionus) were collected (hereafter referred to as Sussemiones) (Eisenmann, 2010). This lineage radiated across North America, Africa, and Siberia, and developed multiple adaptations to a whole range of arid and humid environments (Eisenmann, 2010). Sussemiones were first believed to have become extinct during the Middle Pleistocene as the last known specimen showing typical morpho-anatomical characters dated back to ~500 kya (thousand years ago) in southeastern Siberia, Russia (Vasiliev, 2013). However, DNA results obtained on multiple osseous remains within the radiocarbon range and showing morphological traits reminiscent of the Eurasian Sussemiones species indicated that the lineage in fact survived until the Late Pleistocene (Druzhkova et al., 2017; Orlando et al., 2009; Vilstrup et al., 2013; Yuan et al., 2019). Early publications indicated survival dates 40–50 kya in southeastern Siberia, Russia (Proskuryakova cave) (Orlando et al., 2009; Vilstrup et al., 2013), ~32 kya at the Denisova cave (Druzhkova et al., 2017), and ~12.6 kya at northeastern China (Yuan et al., 2019).

Despite an abundant fossil material, only a limited number of Sussemiones specimens have been investigated for ancient mitochondrial DNA (aDNA). These studies showed that Sussemiones formed a non-caballine lineage that may have diverged first from the lineage ancestral to zebras, hemiones, and donkeys. However, the exact placement of Sussemiones could not be fully resolved and has remained contentious (Heintzman et al., 2017; Orlando et al., 2009; Vilstrup et al., 2013). In this study, we have carried out archaeological excavations in three Holocene sites in China (Honghe, Heilongjiang Province; Muzhuzhuliang, Shaanxi Province; Shatangbeiyuan, Ningxia Province) (Figure 1 and Supplementary file 1a) and uncovered equine samples showing morphological features that may be characteristic of Sussemiones.

Figure 1 with 2 supplements see all
Download asset Open asset
Geographic distribution of E. (Sussemionus) ovodovi specimens investigated at the DNA level and pre-molar and molar morphology.

(A) E. (Sussemionus) ovodovi geographic range. The three red circles indicate the archaeological sites analyzed in this study. The site (Honghe) that delivered the complete genome sequence at 13.4-fold average depth of coverage (HH06D) is highlighted with a square. The black circles indicate sites that provided complete mitochondrial genome sequences in previous studies (Druzhkova et al., 2017; Orlando et al., 2009; Vilstrup et al., 2013; Yuan et al., 2019). The temporal range covered by the different samples analyzed is given in years before present (YBP) and follows the name of each site. Numbers between parentheses indicate the number of samples for which DNA sequence data could be generated. (B) Facies masticatoria dentis of P2, M3, p2, and m3 for the E. (Susseminous) ovodovi samples of the Honghe site analyzed here (a), E. Sussemionus (Eisenmann, 2010) (b), and E. caballus (Laboratory specimen) (c). 1, 4 protocones; 2, 5 metacones; 3 caballine notch. Teeth from the right side are shown, except for E. Sussemionus. The erupted teeth of the samples of the Honghe site appear to be smaller than those of the E. Sussemionus specimen.

Equine assemblages dating to prior to the late Shang dynasty (ca. 3300 years ago) have documented the presence of wild horses in Northern China during the Late Pleistocene (Yuan and Flad, 2006). The taxonomic status and/or stratigraphic placement of the rare material attributed to Neolithic and early Shang contexts remained, however, contentious, leaving the possibility that Sussemiones or other equid taxa (co-)existed in China at the time, especially at the sites investigated in this study. At the Honghe site (47.20°N, 123.62°E), excavation fieldwork of nearly 20,000 m2 has uncovered a late Neolithic settlement site dated to ~3400–4400 years ago, which belonged to a unique, rich fishing and hunting culture characteristic of northeastern China (Figure 1—figure supplement 1). The scale of the moated settlement indicated that there was already social management and relatively high productivity and building technology (Zhang et al., 2020). The Muzhuzhuliang site (38.83°N, 110.50°E) belonged to the ‘Longshan culture.’ It was dated to ~3800–4300 years ago and represents the most complete moated settlement hitherto excavated in the late Neolithic age of Northern China, showing a mixed subsistence economy involving agriculture, animal husbandry, and hunting (Wang et al., 2015). Finally, the Shatangbeiyuan site (35.63°N, 105.11°E) belonged to the early cultural relics of Neolithic ‘Qijia culture,’ which was dated to ~3900–4200 years ago. While millet represented the main crop produced at that time, findings including stone and bone arrowheads have also supported the presence of hunting (Fan et al., 2017). No obvious signs of domestication, including paleopathologies related to horseback riding, bridling, or chariotry (Bendrey, 2007; Taylor and Tuvshinjargal, 2018), were found amongst the equine specimens investigated at the three sites. In contrast, slash marks could be identified on some of the bones (HH13H, HH26H, and MZ104H), together with indications of bone marrow extraction (Figure 1—figure supplement 2). These findings suggest that these specimens were hunted.

In this study, we have sequenced the complete nuclear genome of Sussemiones specimens. This allowed us to not only solve the phylogenetic placement of Sussemiones within the Equus evolutionary tree, but also to time their divergence relative to other non-caballine equids, as well as to reconstruct their demographic trajectory until their extinction during the mid-Holocene.

Results

Archaeological samples and sequencing data

All the equine specimens investigated in this study showed morphological and genetic signatures (short fragments of the mitochondrial hypervariable region) distinct from those of extant horses and donkeys (Figure 1—figure supplement 2). The morphological differences were especially marked in the second and third molars, which appeared to be smaller than in modern horses, and reminiscent of the third molars paracones and metacones observed in Sussemiones specimens (Figure 1B). Combined, these samples were radiocarbon dated to 3456–4460 calibrated years before the present (cal BP), including a mid-second millennium BCE date for the most recent sample, HH13H (3270 ± 30 uncal. BP, i.e., 3456–3616 cal BP) (Supplementary file 1b). They could, thus, represent some of the latest surviving Sussemiones individuals prior to their extinction.

We next aimed at genetically characterizing the taxonomic status of these specimens using high-throughput DNA sequencing technologies. We extracted ancient DNA from a total of 26 specimens and sequenced the whole nuclear genome at ~0.002–13.4 times coverage, including four samples from Honghe that provided 13.4×, 3.9×, 1.1×, and 1.0× nuclear genome (Supplementary file 1a). Comparison of the X chromosome and autosomal coverage revealed the presence of 15 male and 11 female individuals (Supplementary file 1c).

Taxonomic status

To assess whether the sequenced specimens belonged to the same taxonomic group or comprised different species, we carried out a principal component analysis (PCA), including all the equine species sequenced at the genome level (depth of coverage ≥1×) (Figure 2A, Figure 2—figure supplements 1 and 2). For this, we downloaded 11 previously-published equine genomes representing all extant species of equids and the extinct quagga zebra (Huang et al., 2015; Jónsson et al., 2014; Kalbfleisch et al., 2018; Orlando et al., 2013; Renaud et al., 2018; Supplementary file 1d). All the Chinese specimens analyzed in this study were found to cluster together along the first two PCA components, in a group that was distinct from all other equine species (Figure 2A, Figure 2—figure supplement 1) but closer to non-caballine equine species than to the horse (Figure 2A). This suggested that they were all members of a unique taxonomic group, most related to non-caballine equids.

Figure 2 with 9 supplements see all
Download asset Open asset
Genetic affinities within the genus Equus.

The Honghe (HH), Muzhuzhuliang (MZ), and Shatangbeiyuan (BY) specimens are shown in red, while Asian asses, African asses, zebras, and horses are shown in purple, blue, green, and black, respectively. (A) Principal component analysis (PCA) based on genotype likelihoods, including horses and all other extant non-caballine lineages (16,293,825 bp, excluding transitions). Only specimens whose genomes were sequenced at least to 1.0× average depth of coverage are included. (B) Maximum likelihood tree based on six mitochondrial partitions (representing a total of 16,591 bp). Those E. ovodovi sequences that were previously published are shown in red. The tree was rooted using Hippidion saldiasi and Haringtonhippus francisci as outgroups. Node supports were estimated from 1000 bootstrap pseudo-replicates and are displayed only if greater than 50%. The black line indicates the mitochondrial clades A and B. (C) Maximum likelihood tree based on sequences of 19,650 protein-coding genes, considering specimens sequenced at least at a 3.0× average depth of coverage (representing 32,756,854 bp).

Maximum likelihood (ML) and Bayesian phylogenetic analyses including the nearly complete 17 mitochondrial genomes reported in this study (Supplementary file 1a, depth of coverage above 1×) confirmed their clustering with non-caballine equids, within a single monophyletic group that also included five previously characterized Sussemiones specimens (Figure 2B, Figure 2—figure supplements 3 and 4, Supplementary file 1e). This grouping was supported with maximal (100%) bootstrap values. This, and the PCA clustering, indicated that the different excavation sites investigated in this study in fact all provided specimens that belonged to the E. (Sussemionus) ovodovi species.

We also used complete mitogenomes to assess the diversity of maternal lineages present in the E. (Sussemionus) ovodovi lineage. Phylogenetic analyses showed two major clades, of which only one (clade B) was previously characterized. The other clade (A) consisted of two and three individuals from Muzhuzhuliang and Honghe, respectively (Figure 2B).

To further assess phylogenetic affinities, we used the two genomes characterized to at least 3× average depth of coverage (HH04D and HH06D) to place Sussemiones within the equine phylogenetic tree. To achieve this, we used ML phylogenetic reconstruction and an alignment of the coding sequences of the protein-coding genes (Figure 2C, Figure 2—figure supplement 5). This showed that the Chinese ancient specimens branched off before the radiation leading to modern asses and zebras (Figure 2C). Similar tree topologies were recovered using whole-genome SNPs by TreeMix (Pickrell and Pritchard, 2012; Figure 2—figure supplements 6 and 7, Supplementary file 1f). Combined with the analysis of the occlusal surface of the molars, in particular the absence of the caballine notch, the shape of metacones and protocones, and the reduced tooth size (Figure 1B), our analyses consistently supported the material analyzed as small specimens of the extinct Equus (Sussemionus) ovodovi. We, thus, concluded that this lineage survived in China during the Holocene, and until 3477–3637 cal BP, which is ~9000 years after the latest known specimen to date (Druzhkova et al., 2017; Orlando et al., 2009; Vilstrup et al., 2013; Yuan et al., 2019).

Interspecies admixture and demographic modeling

Bifurcating trees fail to capture possible admixture events between lineages. Yet, previous research has unveiled pervasive admixture within equids, even amongst extant equids showing different chromosomal numbers (Jónsson et al., 2014). We thus next assessed whether the genomic data showed evidence for gene flow between Sussemiones and other non-caballine equids. To achieve this, we first applied D-statistics (Soraggi et al., 2018) to the genome sequence underlying 26 individual genomes and detected that E. ovodovi shared an excess of derived polymorphisms with asses than relative to zebras (Figure 3—figure supplements 1 and 2). This suggested that at least one admixture event could have taken place between Sussemiones and the ancestor of asses after their divergence from zebras.

We next leveraged the comparative genome panel and the ancient E. (Sussemionus) ovodovi genome characterized to high depth of coverage (HH06D) to reconstruct the equine demographic history using G-PhoCS (Gronau et al., 2011). More specifically, we first selected members of each equine lineage representing a total number of 10 genomes, conditioned analyses on 15,324 ‘neutral’ loci, and assumed that the genus Equus emerged some 4.0–4.5 Mya, following previous estimates (Orlando et al., 2013). G-PhoCS analyses confirmed previous work indicating that the zebras and asses linages diverged ~2.0 Mya and that the deepest divergence within zebras and asses took place prior to ~1.5 Mya (Jónsson et al., 2014; Figure 3). It revealed that the Sussemiones lineage diverged from the ancestors of extant non-caballine equids ~2.3–2.7 Mya, in line with the fossil record (Eisenmann, 2010). Allowing for migrations provided support for gene flow between Sussemiones and the ancestor of asses and zebras (Figure 3). However, weak to no migrations were detected between Sussemiones and extant equids (Supplementary file 1g). Importantly, the admixture between Sussemiones and the ancestor of extant asses seems to have been stronger than that between Sussemiones and the ancestor of extant zebras, in line with the results of D-statistics. G-PhoCS also supported the presence of significant unidirectional gene flow prior to ~2.3–2.7 Mya, from the horse branch into the ancestral branch to all non-caballine equids, including Sussemiones (probability of gene flow 2.2–8.8%, Supplementary file 1h). This is consistent with previous HMMCoal analyses applied to whole-genome sequences of all extant equine species, which indicated significant gene flow between the deepest branches of the Equus phylogenetic tree until 3.4 Mya, mostly from a caballine lineage into the ancestor of all non-caballine equids (Jónsson et al., 2014).

Figure 3 with 3 supplements see all
Download asset Open asset
Demographic model for extinct and extant equine lineages as inferred by G-PhoCS (Gronau et al., 2011).

Node bars represent 95% confidence intervals. The width of each branch is scaled with respect to effective population sizes (Ne). Independent Ne values were estimated for each individual branch of the tree, assuming constant effective sizes through time. Migration bands and probabilities of migration (transformed from total migration rates) are indicated with solid arrows. The red triangle indicates the earliest Sussemionus evidence found in the fossil record. (Images: E. caballus by Infomastern, E. a. somalicus by cuatrok77, E. kiang by Dunnock_D, E. a. africanus by Jay Galvin, E. hemionus by Cloudtail the Snow Leopard, E. z. hartmannae by calestyo, E. b. quagga by Internet Archive Book Images, E. b. boehmi by GRIDArendal, and E. grevyi by 5of7.)

Dynamic demographic profiles, heterozygosity, and inbreeding levels

We next leveraged the high-coverage Sussemiones genome characterized here to further explore the demographic dynamics until extinction. When modeled as constant through time, population sizes in G-PhoCS indicated that most lineages, including Sussemiones, consisted of small populations, excepting the Burchell’s zebra (Supplementary file 1i). Pairwise sequential Markovian coalescent (PSMC) analyses, however, provided evidence for population size variation through time. First, the PSMC demographic trajectory of Sussemiones was found to diverge from that of other non-caballine equids (specifically, E. hemionus) after ~2.0 Mya, confirming the divergence date estimate retrieved by G-PhoCS (Supplementary file 1i). Second, the Sussemiones demographic trajectory was found to have constantly increased during the last million year but to have remained relatively low for a long period of time, until it reached a peak between 74 and 84 kya. It was, then, followed by an ~45-fold collapse until 13 kya (Figure 4). The lineage maintained extremely reduced population sizes through the Last Glacial Maximal (LGM, 19–26 kya) (Clark et al., 2009) and the Holocene, until it ultimately became extinct.

Figure 4 with 2 supplements see all
Download asset Open asset
Pairwise sequential Markovian coalescent (PSMC) profiles (100 bootstrap pseudo-replicates) of four Eurasian equine species (E. ovodovi HH06D, E. caballus TWI [Kalbfleisch et al., 2018], E. hemionus ONA, and E. kiang KIA) (Jónsson et al., 2014).

The y-axis represents the effective population size (×10,000), and the x-axis is scaled in millions of years before present. Faded lines show bootstrap values.

Importantly, the sample sequenced to sufficient coverage (HH06D) showed minimal heterozygosity and moderate inbreeding levels identified by the fraction of the segments within runs of homozygosity (ROH) (Figure 5). Strikingly, this is true in spite of the increased DNA damage error rates of this genome (Figure 2—figure supplement 9), which likely inflate our estimates. The limited population sizes and resulting genetic diversity, rather than particularly enhanced inbreeding, may, thus, have limited the chances of survival of the species and have ultimately led to extinction.

Figure 5 with 2 supplements see all
Download asset Open asset
Heterozygosity and inbreeding levels of different equine lineages.

(A) Individual heterozygosity outside runs of homozygosity (ROH). (B) Fraction of the genome in ROH. Estimates were obtained excluding transitions and are shown together with their 95% confidence intervals. The colors mirror those from Figure 2.

Discussion

Phylogenetic placement of Equus (Sussemionus) ovodovi

In this study, we have characterized the first nuclear genomes of the now-extinct equine lineage, E. (Sussemionus) ovodovi, the last surviving member of the subgenus Sussemionus. We demonstrated that this lineage survived in China well into the Holocene with the most recent specimens analyzed dating to ~3456–3616 cal BP. This is almost 9000 years after the latest specimens previously documented in the fossil record (Druzhkova et al., 2017; Vilstrup et al., 2013; Yuan et al., 2019). Our work, thus, shows that Sussemionus represents the last currently known Equus subgenus to become extinct. Our work also adds to the list of recently identified members of the horse family that were still alive at the time horses and donkeys were first domesticated, ~5500 years ago (Fages et al., 2019; Gaunitz et al., 2018; Rossel et al., 2008). In contrast to those divergent members that were identified in Siberia (Equus lenensis) and Iberia (IBE), which both belonged to the horse species (Fages et al., 2019; Schubert et al., 2014a), Sussemiones members were most closely related to non-caballine equids. This is in agreement with previous studies (Der Sarkissian et al., 2015; Druzhkova et al., 2017; Heintzman et al., 2017; Orlando et al., 2009; Vilstrup et al., 2013; Yuan et al., 2019), which could, however, not fully resolve the exact phylogenetic placement of this species within non-caballines as topological tests based on mitochondrial genomes received low confidence support (Der Sarkissian et al., 2015; Druzhkova et al., 2017; Heintzman et al., 2017; Orlando et al., 2009; Vilstrup et al., 2013; Yuan et al., 2019). Our study solved this question by reporting the first whole-genome phylogeny of Sussemiones, which confirmed with maximal bootstrap support this species as a basal lineage of non-caballine equids.

Suitable habitat and geographic distribution

Previous zooarchaeological and environmental research indicated an ecological range for Sussemiones overlapping with the grasslands located east of the Altay Mountains and west of the Yenisei River during the Late Pleistocene (Khenzykhenova et al., 2016; Malikov, 2016; Plasteeva et al., 2015; Shchetnikov et al., 2015; Slon et al., 2018). Recent research also reported species occurrence in northeastern China (~12,600–40,200 YBP), where similar climatic and ecological conditions were found at the time (Yuan et al., 2019). It could, thus, be speculated that Sussemiones were adapted to an environment with moderately dry climatic conditions and steppe landscapes (Yuan et al., 2019). However, our study identified Sussemiones specimens in three late Holocene sites from China characterized by mild and humid environmental conditions. Additionally, two distinct mitochondrial haplogroups from 22 individuals have been defined from the six known sites, suggesting possible population structure across various geographic areas and adaptation to local environments. It also suggests that the species could adapt to a wider variety of habitats than previously hypothesized and rejects the contention that the species became extinct as it could not survive in warmer climatic conditions (Yuan et al., 2019).

Interestingly, the Sussemiones specimens identified in this study were excavated from sites in northeastern China located at almost the same latitude as those Sussemiones localities known so far from Russia, but also at lower latitudes (Figure 1A). This implies that the geographic range of E. ovodovi was larger than previously expected and included at least Northern China and Southern Siberia. Although in the absence of identified fossils from Mongolia, given that there is a lack of mitochondrial phylogeographic structure, we could speculate that the two regions were in contact at least maternally. Further work is necessary to establish whether or not the species survived in other pockets both within and outside China.

Demographic history with ancestral interspecific admixture

Our analyses reveal that the divergence between Sussemiones and the most recent common ancestor of all extant non-caballine equids took place ~2.3–2.7 Mya, prior to the divergence of zebras and asses. Post-divergence admixture events with the lineage ancestral to asses and zebras, on the one hand, and the lineage ancestral to all extant zebras, were also identified (Figure 3 and Supplementary file 1h). Our results, thus, reveal non-caballine ancestral lineages occupying partly sympatric distributions that were, consequently, different than those of their descendants, in which zebras are restricted in Africa and Asian asses in Asia. Whether the admixture events identified here directly involved the Sussemiones lineage or one (or more) ghost lineage(s) closely related to Sussemiones requires further research.

Limited genetic diversity before extinction

The demographic profile of Sussemiones shows that after the peak of population size culminating ~74 kya, Sussemiones went through a slow and continuous decline until 13 kya (Figure 4). This time period encompasses several major climate changes (especially the LGM, ~19–26 kya) (Clark et al., 2009) and the great human expansion to Eurasia (~35–45 kya) (Henn et al., 2012). The effective size of Sussemiones populations that survived in Northern China until at least ~3500 years ago, remained extremely small, as indicated by their extremely reduced heterozygosity levels compared to other extant and extinct equine species. As the inbreeding levels were not particularly high compared to some members of endangered equine species (Figure 5), the reduced genetic diversity available in the lineage may have compromised the long-term survival of the lineage, in a process partly reminiscent of what was previously described for the woolly mammoth (Palkopoulou et al., 2015). And considering the rapid expansion of domestic horses across Eurasia from about 2000 BC (Librado et al., 2021), this lineage was ultimately replaced under growing anthropogenic stress.

In conclusion, our study clarifies the phylogenetic placement, speciation timing, and evolutionary history of the now-extinct Sussemionus equine subgenus. This group did not remain in reproductive isolation from other equine lineages, but contributed to the genetic makeup of the ancestors of present-day asses, while receiving genetic material from the ancestors of African zebras. This supports geographic distributions at least partly overlapping at the time, thus, not identical to those observed today. The species demographic trajectory experienced a steady decline from ~74 kya and during a period witnessing both important climatic changes and the Great human expansion across Asia (Henn et al., 2012). It survived with minimal genetic diversity the Pleistocene-Holocene transition, and for at least eight millennia before it became extinct, which provides insights into the survival potential of large animals since the Holocene. Given the persistence of Sussemiones throughout the third and second millennia BCE, archaeologists must be exceedingly careful while assigning Asian zooarchaeological material to equine taxa until this period.

Materials and methods

Genome sequencing

Request a detailed protocol

Minimum number of individuals (MNI) was determined by assigning the frequency of hip bone and was calculated from the acetabular bone to avoid double counting. MNI was estimated to 31 individuals at Honghe, 4 at Muzhuzhuliang, and 4 at Shatangbeiyuan. DNA preservation conditions were compatible with the recovery of ancient DNA sequences from only 20 of the 31 Honghe samples, 3 of the 4 Muzhuzhuliang samples, and 3 of the 4 Shatangbeiyuan samples (Supplementary file 1a).

All pre-PCR procedures were conducted in a dedicated ancient DNA laboratory at Jilin University (JLU) that is physically separated from the post-PCR laboratory. To remove potential contaminant DNA, working areas and benches were frequently cleaned with bleach and UV exposure. Lab experiments were carried out wearing full-body suits, facemasks, and gloves. To detect contamination, mock controls were included at each experimental step, including DNA extraction, DNA library preparation, and PCR setup.

Prior to DNA extraction, the outer surface of the sample was cleaned with a brush. The cleaned sample was subsequently cut into smaller pieces and soaked in 10% bleach for 20 min, rinsed with ethanol and distilled water, and then UV-irradiated for 30 min on each side. Finally, powder was obtained using a dental drill (Traus 204, Korea). Ancient DNA was extracted from the sample powder by using a modified silica spin column method (Yang et al., 1998), in the dedicated ancient DNA facilities from JLU. For each specimen, a total of 200 mg powder was added with 3.9 ml EDTA (0.465 mol/L) and placed in the refrigerator at 4°C for 12 hr for decalcification, and then 0.1 mL proteinase K (0.4 mg/mL) were added and incubated overnight in a rotating hybridization oven at 50℃ (220 rpm). After centrifugation, the supernatant was transferred into an Amicon Ultra-4 centrifugal filter device (Merck Millipore Ltd, 10,000 Nominal Molecular Weight Limit), reduced to less than 100 µL, and purified with QIAquick PCR Purification Kit (QIAGEN), according to the manual instructions.

Before preparation of DNA libraries, we first PCR-targeted short fragments of the mitochondrial hypervariable region to select those samples positive for the presence of equine DNA (which was further confirmed through Sanger sequencing). For this, we used the oligonucleotide primers L15473 5′-CTTCCCCTAAACGACAACAA-3′ and reverse primer H15692 5′-TTTGACTTGGATGGGGTATG-3′; and forward primer L15571 5′-AATGGCCTATGTACGTCGTG-3′ and reverse primer H15772 5′-GGGAGGGTTGCTGATTTC-3′ from Juan et al., 2007, and the amplification conditions therein.

Double-stranded single-indexed libraries were prepared using NEBNext Ultra Ⅱ DNA Library Prep Kit for Illumina (NEB #E7645S) and NEBNext Multiplex Oligos for Illumina Index Primers Set 1 and 2 (NEB #E7335S, #E7500S), following the manufacturer’s instructions with minor modifications. Specifically, the extracted DNA (50 µL) were end-repaired and A-tailed by adding 7 μL of NEBNext Ultra II End Prep Reaction Buffer and 3 μl of NEBNext Ultra II End Prep Enzyme Mix, and incubated for 40 min at 20°C and then 30 min at 65°C. The adaptor was ligated to the dA-tailed DNA fragments by adding 30 µL of NEBNext Ultra II Ligation Master Mix, 1 µL of NEBNext Ligation Enhancer and 2.5 µL of NEBNext Adaptor for Illumina (dilution 1:10), and incubated for 20 min at 20°C. The adaptor was then linearized by adding 3 µL of USER Enzyme and performing an incubation for 15 min at 37°C. The adaptor-ligated DNA were cleaned without size selection using the MinElute PCR Purification Kit (QIAGEN, Germany), following the instructions provided by the manufacturer. PCR enrichment was performed by using 30 µL of NEBNext Ultra II Q5 Master Mix, 1 µL of Index Primer, 1 µL of Universal PCR Primer, and 18 µL of adaptor-ligated DNA. PCR cycling conditions comprised an initial denaturation at 98°C for 30 s, 14–16 cycles of 98°C for 10 s, 65°C for 75 s, and a final extension at 65°C for 5 min. PCR-amplified DNA libraries were purified using Agencourt AMPure XP Beads, following the manufacturer’s instructions, and Illumina sequencing was performed on HiSeq X Ten platform using 150 bp paired-end reads. Overall, we sequenced a total of 28 DNA libraries and generated 2,727,843,803 read pairs (https://www.ebi.ac.uk/ena/browser/view/PRJEB44527?show=reads).

Radiocarbon dating

Request a detailed protocol

Radiocarbon dating of the samples was performed at the Beta Analytic Radiocarbon Dating Laboratory, Miami, FL. Bone or tooth pieces about 2 g were sampled in the bone and sent for subsequent dating of collagen (not ultrafiltered). Calibration was carried out using OxCalOnline (https://c14.arch.ox.ac.uk/oxcal.html) and the IntCal20 calibration curve. Calibrated dates are provided in Supplementary file 1b.

Data processing

Request a detailed protocol

Sequencing reads were processed and aligned against the horse (EquCab3.0 Kalbfleisch et al., 2018) and donkey (Renaud et al., 2018) reference genomes using the PALEOMIX pipeline (Schubert et al., 2014b) with default parameters, except that we followed the recommendations from Schubert et al., 2012 and disabled seeding. Briefly, paired-end (PE) reads longer than 25 nucleotides were trimmed with AdapterRemoval v2.2 (Schubert et al., 2016) and aligned against the reference genomes using BWA (Li and Durbin, 2009), retaining alignments with mapping qualities superior to 25. PCR duplicates were then removed using Picard (http://broadinstitute.github.io/picard/) (Broad Institute, 2019). Finally, all ancient and modern reads were locally realigned around indels using GATK (McKenna et al., 2010).

Postmortem DNA damage and average sequencing error rates were determined with mapDamage2.0 (Jónsson et al., 2013; Figure 2—figure supplement 8) and ANGSD (Korneliussen et al., 2014; Figure 2—figure supplement 9), respectively. Further rescaling and trimming procedures were implemented following Gaunitz et al., 2018 to limit the impact of remnant nucleotide misincorporations in subsequent analyses. For each of the DNA libraries examined, the base composition of the position preceding read starts on the horse reference genome showed an excess of Guanine and, to a lesser extent, of Adenine residues (Figure 2—figure supplement 8). This is in line with depurination driving postmortem DNA fragmentation (Briggs et al., 2007). Additionally, error rate estimates for each nucleotide substitution class indicated the predominance of C→T and G→A misincorporations (Figure 2—figure supplement 9). Such misincorporation rates were particularly inflated towards read ends, but not read starts (Figure 2—figure supplement 8). This is in line with the DNA nucleotide misincorporation profiles expected for the type of DNA library constructed (Seguin-Orlando et al., 2015), which was caused by the Q5 polymerase being unable to read through 5' uracils, thereby excluding the typical 5' excess of C-to-T. MapDamage profiles were, thus, consistent with Cytosine deamination at 5'-overhanging ends as the most prominent postmortem DNA degradation reactions (Jónsson et al., 2013).

GATK HaplotypeCaller was used to obtain individual gvcf files with “--minPruning 1 --minDanglingBranchLength 1” to increase sensitivity. Then, we employed GATK GenotypeGVCFs for genotyping with the option “--includeNonVariantSites” in order to retain non-variant loci. The vcf files were further filtered in TreeMix and G-PhoCS analysis.

Principal component analysis (PCA)

Request a detailed protocol

The genotype likelihood framework implemented in ANGSD helped mitigate various error rates in ancient and modern genomes. Using EquCab3 (Kalbfleisch et al., 2018) as the reference genome, ANGSD was run using the following options: “-only_proper_pairs 1 -uniqueOnly 1 -remove_bads 1 -minQ 20 -minMapQ 25C 50 -baq 1 -skipTriallelic 1GL 2 -SNP_pval 1e-6 -rmTrans 1”. This provided a dataset consisting of a total of 16,293,825 transversions when the horse was included, and 10,094,431 transversions when the horse was excluded (i.e., when analyses were restricted to non-caballine genomes only). In these analyses, only specimens sequenced to an average depth of coverage ≥1× were retained. PCA was carried out using the PCAngsd package (Meisner and Albrechtsen, 2018; Figure 2A). To assess the impact of potential reference bias, all analyses were repeated after mapping the sequence data against the donkey reference (Figure 2—figure supplement 2).

Phylogenetic inference

Mitochondrial phylogeny

Request a detailed protocol

Cleaned reads were mapped against the horse mitochondrial genome (GenBank accession no. NC_001640), following the same procedure as when mapping against the nuclear genome. Samples showing an average depth of coverage <1× were disregarded, leaving a total of 17 individuals for further analyses. After removing duplicates, consensus mitochondrial sequences were generated using ANGSD (-doFasta 2 -doCounts 1 -setMinDepth 3 -uniqueOnly 1 -remove_bads 1 -minQ 25 -minMapQ 25). Multiple alignment was performed together with the comparative mtDNA sequences downloaded from GenBank (Supplementary file 1e) using MUSCLE v3.8.31 (Edgar, 2004), with default parameters. The alignments were then split into six partitions (first, secoond, and third codon positions, rRNA, tRNA, and control region) by Partition Finder v2.1.1 (Lanfear et al., 2012).

Two ML trees based on all six partitions and excluding the control region (positions 15,469–16,660 of the horse reference mitochondrial genome) were both reconstructed using RAxML-NG v.0.9.0 (Kozlov et al., 2019) with GTR+GAMMA substitution model. A total of 1000 bootstrap pseudo-replicates were carried out to assess node robustness (Figure 2—figure supplement 3). BEAST 2.6.6 (Bouckaert et al., 2019) was used to perform Bayesian phylogenetic reconstruction and to estimate split times. The six partitions described above were used, for which the best substitution model was determined using modelgenerator (version 0.85, Keane et al., 2006) and a Bayesian information criterion. We calibrated the tree using tip dates (see Supplementary file 1j) and an age of 4–4.5 Mya for the root of crown group E. caballus (normal prior, mean 4.25 Mya, SD: 0.15 Mya) (Orlando et al., 2013). We applied together with the birth-death model and a relaxed molecular clock (log normal) for 1000 million generations (sampling frequency = 1 every 1000), while forced monophyly for all main lineages, including donkeys, hemiones, horses, ovodovi, and zebras. Convergence was assessed visually using Tracer v1.6 (with all individual ESS >200), and posterior date estimates were retrieved using 25% as burn-in. The final consensus tree was produced by TreeAnnotator 2.6.6 (Drummond and Rambaut, 2007) as the maximum clade credibility tree from 100,000 randomly sampled trees obtained using LogCombiner v2.6.6 (Bouckaert et al., 2019) (burn-in = 20%). The final tree was plotted using ITOL (Letunic and Bork, 2016; Figure 2—figure supplement 4).

Autosomal phylogeny

Request a detailed protocol

As for autosomes, we reconstructed an ML phylogenetic tree as implemented in the PALEOMIX phylo_pipeline, which is dedicated to phylogenomic reconstructions (Schubert et al., 2014b). This analysis was based on the coding sequence (CDS) of protein-coding genes annotated in EquCab3.0, partitioning data according to first, second, and third codon positions. ML phylogenetic inference was performed using ExaML v3.0.21 (Kozlov et al., 2015) and RAxML v8.2.12 (Stamatakis, 2014) under the GAMMA substitution model with 100 bootstrap pseudo-replicates (Figure 2C, Figure 2—figure supplement 5A). We also repeated the same procedure after mapping against the donkey reference genome, which returned the same topology (Figure 2—figure supplement 5B).

Additionally, we extracted biallelic single-nucleotide polymorphisms (SNPs) from the dataset generated in the section ‘Variant calling’ using bcftools v1.9 (Li et al., 2009). Both variant datasets obtained following mapping against the horse and donkey reference genomes were used in this analysis to rule out reference bias. We applied filters composed of minimum Phred-scaled quality score quality (QUAL) = 20, sites for all individuals below 2 or twice the mean coverage, and allowed up to three individuals with missing data per site. After disregarding transitions, a total of 18,803,101 (mapping against horse genome) and 19,459,070 (mapping against donkey genome) transversions were finally used as input for TreeMix (Pickrell and Pritchard, 2012) with parameters “-k 500 -root TWI”, and considering an increasing number of migrations edges (0 ≤ m ≤ 3; Figure 2—figure supplements 6 and 7, Supplementary file 1f).

Admixture analyses with D-statistics

Request a detailed protocol

D-statistics were calculated to investigate potential introgression between E. ovodovi and other non-caballines (Figure 3—figure supplement 1) using the doAbbababa2 program in ANGSD (Soraggi et al., 2018). Individuals were grouped according to their respective species. D-statistics were computed in the form (((H1, H2), H3), Outgroup) considering only the autosomal sites from bam files mapping against the horse reference with the following options: “-minQ 20 -minMapQ 25 -remove_bads 1 -only_proper_pairs 0 -uniqueOnly 1 -baq 1C 50”. The horse reference genome was used as the Outgroup. H1 and H2 denoted any non-caballine genomes except E. ovodovi, while H3 denoted the E. ovodovi. Confidence intervals were estimated applying a jackknife procedure and 5 Mb windows. Z-scores with absolute values higher than 3 were considered to be statistically significant. To rule out possible reference bias, we also rerun the same analysis using sequence alignments against the donkey reference genome (Figure 3—figure supplement 2).

G-PhoCS demographic model

Data preparation and filtering

Request a detailed protocol

In order to model the equine evolutionary history, we selected a total of 10 individuals representing each individual lineage and used their high-coverage genomes as input for G-PhoCS (Gronau et al., 2011). Genotypes were called by GATK and candidate ‘neutral’ loci were identified by applying the following filters:

  1. The simple repeats track available for the reference genome was obtained from Ensembl v99 release; corresponding regions were masked.

  2. All exons of protein-coding genes were discarded together with their 10 kb flanking regions; this was done based on the GTF format annotation file of the reference genome available from Ensembl v99 Genome Browser.

  3. We identified conserved noncoding elements (CNEs) using phastCons scores (based on the 20-way Conservation track provided for the mammal clade according to the genomic coordinates of the human reference) downloaded from the Table Browser of UCSC. All CNEs and their 100 bp flanking regions were masked using liftOver to convert human genome coordinates into EquCab3.0 horse genome coordinates.

  4. Exons of noncoding RNA genes together with their 1 kb flanking regions were removed, based on the annotations available for the reference genome.

  5. Gaps in the reference genome were disregarded.

Besides the various filters described above, regions/sites likely to be enriched for misaligned bases, and to have high false-negative rates during read alignment or variant detection, were masked as missing data. More specifically, different individuals could be treated differently depending on the genotyping results obtained as described in the section ‘Variant calling,’ depending on the presence of (1) indels, (2) triallelic sites, (3) positions with depth of coverage twice the mean depth recorded for each individual, and (4) transition sites.

We selected 1 kb loci located with minimum inter-locus distance of 30 kb from the intervals that pass all the criteria described above. Then, consensus sequences were generated for each individual from the vcf file generated in the section ‘Variant calling’ using bcftools ‘consensus’ command, with IUPAC codes indicating heterozygous genotypes (--iupac-codes) and ‘N’ representing masked sites (--mask and --missing ‘N’).

Finally, we excluded contiguous intervals if the total amount of missing bases was greater than 50% of the region length, resulting in a final collection of 15,324 loci using the horse reference genome (autosomes only). Neighbor-joining trees were constructed to confirm the topology before the inferring the population divergence (Figure 3—figure supplement 3).

MCMC setup

Request a detailed protocol

We used default global settings (Gronau et al., 2011), including a gamma prior distribution (α = 1, β = 10,000) for all mutation-scaled population sizes (θ) and a gamma prior distribution (α = 0.002, β = 0.00001) for all mutation-scaled migration rate (m). The initial parameter value of mutation-scaled divergence times (τ) was first set individually for each population. Then, we ran ~100,000–200,000 iteration tests and manually evaluated the convergence by checking the achieve acceptance ratios (i.e., accept if around 30–70%) or using Tracer v1.6 (http://tree.bio.ed.ac.uk/software/tracer/). For each test, we updated the input of the initial τ and all fine-tuned parameters based on previous results to get the appropriate value. The final results in Figure 3 are based on 500,000 MCMC iterations, considering the first 10% as burn-in.

Parameter calibration

Request a detailed protocol

We assumed an average generation time (g) of 8 years. The coalescent time of the Equus (4.0–4.5 Mya) (Orlando et al., 2013) was used to bound the mutation rate μ (per site per year). Effective population sizes (Ne) and divergence times (T) were estimated by scaling θ and τ parameter using g and μ (Supplementary file 1g), and the following formula: Ne = θ/(4 μg) and T = τ/μ (Gronau et al., 2011).

Inferring gene flow

Request a detailed protocol

Total migration rates (M) were estimated by a mutation-scaled parameter (m) given by M = mτm, where τm is the mutation-scaled time span of the migration band. The total migration rate gives an accumulated migration rate over a long period of time, which can be superior to 100%. We then converted such rates, M, into a probability of migration using the formula P = 1-e-M (where p is the probability of gene flow), according to the method presented in vonHoldt et al., 2016.

The migration model implemented in G-PhoCS makes it possible to detect gene flow between any two lineages by introducing migration bands manually to the demographic model. However, it remains difficult to detect weak migration events. Additionally, scenarios including a large number of migration bands can lead to spurious results. To address this, we first inferred a demographic model with no migration bands, and then introduced several migration bands corresponding to five independent scenarios (Supplementary file 1g). A significant migration band was considered supported if both the 95% Bayesian credible interval of total migration rate (M) did not include 0% values and if the mean value of M was estimated to be greater than 0.03%.

Settings for the migration bands between extant caballines are based on previous research (Jónsson et al., 2014). The significant migration band from horse to the non-caballine ancestor was identified (Supplementary file 1g), in line with previous work (Jónsson et al., 2014). However, no other non-negligible (M > 3%) migration bands were found in our analyses (Supplementary file 1g).

We then tried to estimate the migration events between E. ovodovi and other branches. We added all possible migration bands between E. ovodovi and extant non-caballine branches into the demographic model except the migration bands between E. ovodovi and the ancestor of extant non-caballines as the model is often underpowered to infer migration between sister populations. All of the migration bands were separated into four demographic models. Only three migration bands were shown to be significant (Supplementary file 1g).

Finally, the total four migration bands were combined into one demographic model (Supplementary file 1h) and compared the estimates to the one including no migration (Supplementary file 1i).

We caution that the analyses carried out using TreeMix and G-PhoCS returned partly discordant results. This may be due to TreeMix modeling pulses of admixture in contrast to G-PhoCS, in which situations of continuous gene flow can be accommodated. Additionally, gene flow affecting the two deepest tree branches can be directly accommodated by reducing their divergence. Therefore, the deep admixture inferred by G-PhoCS from the caballine branch into the ancestral branch of Sussemiones and other non-caballine equids cannot be expected to be identified through an individual migration edge with TreeMix as this could simply be modeled through a more limited divergence between both underlying lineages. The same holds true for the asymmetric gene flow inferred by G-PhoCS between Sussemiones and the branch ancestral to all extant asses; TreeMix is likely to only identify the resulting unidirectional contribution of these admixtures, which mainly sources to the branch ancestral to extant asses into Sussemiones; since G-PhoCS also infers additional admixture from the branch ancestral to extant zebras into Sussemiones, we can expect TreeMix to accommodate both sources of gene flow through a reduced divergence between Sussemiones and the branch ancestral to stenonines. Finally, it may reflect limitations pertaining to the two underlying data sets utilized, consisting, on the one hand, to the whole-genome SNP panel for TreeMix, further filtered for 15,324 candidate ‘neutral’ loci in G-PhoCS.

Demographic trajectories with PSMC

PSMC analyses

Request a detailed protocol

In order to reconstruct the past demographic dynamics of the E. ovodovi lineage, we applied the PSMC algorithm (version 0.6.5-r67) (Li and Durbin, 2011) to the sample HH06D (12.0×, mapping against horse reference), as well as three other Eurasian equine species (E. caballus TWI, E. hemionus ONA, and E. kiang KIA).

We first obtained the diploid consensus sequences after mapping against the horse genome for the autosomes of each specimens using bcftools ‘mpileup’ command and the ‘vcf2fq’ command from vcfutils.pl with the following filters: mapping quality ≥ 25; adjust mapping quality = 50; minimum depth of coverage = 8; maximum depth of coverage ≤ 99.5% quantile of the coverage distribution; minimum RMS mapping quality = 10; filtering window size of indels = 5.

After filtering the bases with Phred quality scores strictly lower than 35, we ran PSMC with the following command: ‘psmc -N25 -t15 -r5 -p “4+25*2+4 + 6”’. Calibration was carried out using a generation time of 8 years and mutation rate of 7.242 × 10–9 per generation per site, following previous work (Jónsson et al., 2014). However, as for the misincorporation pattern and high error rate of HH06D (Figure 2—figure supplements 8 and 9), we also performed analyses without transitions using mutation rates of 2.3728 × 10–9 that was obtained assuming that the most recent common ancestor of living equine species emerged 4 Mya (Orlando et al., 2013).

We found a great expansion of HH06D in the past 50,000 years when retaining transitions but not when conditioning on transversions (Figure 4—figure supplement 1). The former is thus likely spurious and at least partly driven by severe postmortem DNA damage signatures in the sequence data. We therefore only used the latter inference when considering the ancient HH06D specimen.

False-negative rate correction

Request a detailed protocol

The HH06D genome (12.0×) was corrected assuming a uniform false-negative rate (uFNR) following Orlando et al., 2013 as the average depth of coverage is lower than the recommended 20×. To identify the correction value of uFNR for HH06D, we randomly downsampled reads of the SOM genome (21.0×), using DownsampleSam function of Picard Tools to downscale sequence data to the same average depth of coverage as that obtained for HH06D. This indicated that a value of 0.22 was the most suitable uFNR value for rescaling the HH06D PSMC profile (Figure 4—figure supplement 2A). The KIA and the ONA genomes, which also showed limited coverage, were also rescaled following the same procedure (Figure 4—figure supplement 2B and C). Finally, PSMC confidence intervals were assessed from 100 bootstrap pseudo-replicates (Figure 4).

Heterozygosity inference and inbreeding

Request a detailed protocol

Global heterozygosity rates and inbreeding levels were inferred for high-coverage individuals (>10×) using ROHan (Renaud et al., 2019) with default parameters, except that transitions were excluded (--tvonly) (Figure 5—figure supplement 1). To limit the impact of remnant misincorporations, we used the attached estimateDamage.pl script to estimate damage for all ancient samples prior to heterozygosity computation. Inbreeding was co-estimated together with genome-wide heterozygosity levels from the total ROH length (Figure 5—figure supplement 2).

Data availability

Sequencing data have been deposited in the European Nucleotide Archive under the accession number PRJEB44527.

The following data sets were generated
    1. Zhu S
    2. Gong M
    3. Zhang N
    4. Wen J
    5. Liang Q
    6. Sun W
    7. Shao X
    8. Guo Y
    9. Cai Y
    10. Zheng Z
    11. Zhang W
    12. Hu S
    13. Wang X
    14. Tian H
    15. Li Y
    16. Liu W
    17. Yang M
    18. Yang J
    19. Wu D
    20. Orlando L
    21. Jiang Y
    22. Cai D
    (2021) European Nucleotide Archive
    ID PRJEB44527. Our sequence data provided 26 mitochondrial genomes and 3 complete nuclear genomes for Equus (Sussemionus) ovodovi.
    https://www.ebi.ac.uk/ena/browser/view/PRJEB44527?show=reads
The following previously published data sets were used
    1. Renaud G
    2. Petersen B
    3. Seguin-Orlando A
    4. Bertelsen M F
    5. Waller A
    6. Newton R
    7. Paillot R
    8. Bryant N
    9. Vaudin M
    10. Librado P
    11. Orlando L
    (2018) European Nucleotide Archive
    ID PRJEB24845. This study aims at improving the genome reference of the domestic donkey using the Chicago/HiRiSe technology.
    https://www.ebi.ac.uk/ena/browser/view/PRJEB24845?show=reads
    1. Jonsson H
    (2014) European Nucleotide Archive
    ID PRJEB7446. Speciation with gene flow in equids despite extensive chromosomal plasticity.
    https://www.ebi.ac.uk/ena/browser/view/PRJEB7446?show=reads
    1. Ginolhac A
    (2013) NCBI BioSample
    ID SAMN02179859. General Sample for Equus asinus asinus, Willy.
    https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN02179859
    1. Dugarjaviin M
    (2014) NCBI BioSample
    ID SAMN03010637. Model organism or animal sample from Equus hemionus.
    https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN03010637
    1. Achilli A
    (2012) NCBI GenBank
    ID 347361635. Mitochondrial genomes from modern horses reveal the major haplogroups that underwent domestication.
    https://www.ncbi.nlm.nih.gov/nuccore/347361635/
    1. Vilstrup JT
    (2013) NCBI GenBank
    ID JX312719. Mitochondrial phylogenomics of modern and ancient equids.
    https://www.ncbi.nlm.nih.gov/nuccore/JX312719
    1. Vilstrup JT
    (2013) NCBI GenBank
    ID JX312721. Mitochondrial phylogenomics of modern and ancient equids.
    https://www.ncbi.nlm.nih.gov/nuccore/JX312721
    1. Vilstrup JT
    (2013) NCBI GenBank
    ID JX312725. Mitochondrial phylogenomics of modern and ancient equids.
    https://www.ncbi.nlm.nih.gov/nuccore/JX312725
    1. Vilstrup JT
    (2013) NCBI GenBank
    ID JX312730. Mitochondrial phylogenomics of modern and ancient equids.
    https://www.ncbi.nlm.nih.gov/nuccore/JX312730
    1. Vilstrup JT
    (2013) NCBI GenBank
    ID JX312732. Mitochondrial phylogenomics of modern and ancient equids.
    https://www.ncbi.nlm.nih.gov/nuccore/JX312732
    1. Vilstrup JT
    (2013) NCBI GenBank
    ID JX312734. Mitochondrial phylogenomics of modern and ancient equids.
    https://www.ncbi.nlm.nih.gov/nuccore/JX312734
    1. Der Sarkissian C
    (2015) NCBI GenBank
    ID KM881671. Mitochondrial genomes reveal the extinct Hippidion as an outgroup to all living equids.
    https://www.ncbi.nlm.nih.gov/nuccore/KM881671
    1. Libradoa P
    (2015) NCBI GenBank
    ID KT368725. Tracking the origins of Yakutian horses and the genetic basis for their fast adaptation to subarctic environments.
    https://www.ncbi.nlm.nih.gov/nuccore/KT368725
    1. Orlando L
    (2016) NCBI GenBank
    ID KT757740. Recalibrating Equus evolution using the genome sequence of an early Middle Pleistocene horse.
    https://www.ncbi.nlm.nih.gov/nuccore/KT757740
    1. Orlando L
    (2016) NCBI GenBank
    ID KT757741. Recalibrating Equus evolution using the genome sequence of an early Middle Pleistocene horse.
    https://www.ncbi.nlm.nih.gov/nuccore/KT757741
    1. Druzhkova AS
    2. Makunin AI
    3. Vorobieva NV
    4. Vasiliev SK
    5. Ovodov ND
    6. Shunkov MV
    7. Trifonov VA
    8. Graphodatsky AS
    (2017) NCBI GenBank
    ID KY114520. Complete mitochondrial genome of an extinct Equus (Sussemionus) ovodovi specimen from Denisova cave (Altai, Russia).
    https://www.ncbi.nlm.nih.gov/nuccore/KY114520
    1. Xu X
    2. Gullberg A
    3. Arnason U
    (2016) NCBI GenBank
    ID X97337. The complete mitochondrial DNA (mtDNA) of the donkey and mtDNA comparisons among four closely related mammalian species-pairs.
    https://www.ncbi.nlm.nih.gov/nuccore/X97337

References

    1. Fages A
    2. Hanghøj K
    3. Khan N
    4. Gaunitz C
    5. Seguin-Orlando A
    6. Leonardi M
    7. McCrory Constantz C
    8. Gamba C
    9. Al-Rasheid KAS
    10. Albizuri S
    11. Alfarhan AH
    12. Allentoft M
    13. Alquraishi S
    14. Anthony D
    15. Baimukhanov N
    16. Barrett JH
    17. Bayarsaikhan J
    18. Benecke N
    19. Bernáldez-Sánchez E
    20. Berrocal-Rangel L
    21. Biglari F
    22. Boessenkool S
    23. Boldgiv B
    24. Brem G
    25. Brown D
    26. Burger J
    27. Crubézy E
    28. Daugnora L
    29. Davoudi H
    30. de Barros Damgaard P
    31. de Los Ángeles de Chorro Y de Villa-Ceballos M
    32. Deschler-Erb S
    33. Detry C
    34. Dill N
    35. do Mar Oom M
    36. Dohr A
    37. Ellingvåg S
    38. Erdenebaatar D
    39. Fathi H
    40. Felkel S
    41. Fernández-Rodríguez C
    42. García-Viñas E
    43. Germonpré M
    44. Granado JD
    45. Hallsson JH
    46. Hemmer H
    47. Hofreiter M
    48. Kasparov A
    49. Khasanov M
    50. Khazaeli R
    51. Kosintsev P
    52. Kristiansen K
    53. Kubatbek T
    54. Kuderna L
    55. Kuznetsov P
    56. Laleh H
    57. Leonard JA
    58. Lhuillier J
    59. Liesau von Lettow-Vorbeck C
    60. Logvin A
    61. Lõugas L
    62. Ludwig A
    63. Luis C
    64. Arruda AM
    65. Marques-Bonet T
    66. Matoso Silva R
    67. Merz V
    68. Mijiddorj E
    69. Miller BK
    70. Monchalov O
    71. Mohaseb FA
    72. Morales A
    73. Nieto-Espinet A
    74. Nistelberger H
    75. Onar V
    76. Pálsdóttir AH
    77. Pitulko V
    78. Pitskhelauri K
    79. Pruvost M
    80. Rajic Sikanjic P
    81. Rapan Papeša A
    82. Roslyakova N
    83. Sardari A
    84. Sauer E
    85. Schafberg R
    86. Scheu A
    87. Schibler J
    88. Schlumbaum A
    89. Serrand N
    90. Serres-Armero A
    91. Shapiro B
    92. Sheikhi Seno S
    93. Shevnina I
    94. Shidrang S
    95. Southon J
    96. Star B
    97. Sykes N
    98. Taheri K
    99. Taylor W
    100. Teegen W-R
    101. Trbojević Vukičević T
    102. Trixl S
    103. Tumen D
    104. Undrakhbold S
    105. Usmanova E
    106. Vahdati A
    107. Valenzuela-Lamas S
    108. Viegas C
    109. Wallner B
    110. Weinstock J
    111. Zaibert V
    112. Clavel B
    113. Lepetz S
    114. Mashkour M
    115. Helgason A
    116. Stefánsson K
    117. Barrey E
    118. Willerslev E
    119. Outram AK
    120. Librado P
    121. Orlando L
    (2019) Tracking Five Millennia of Horse Management with Extensive Ancient Genome Time Series
    Cell 177:1419–1435.
    https://doi.org/10.1016/j.cell.2019.03.049
    1. Fan J
    2. Wang X
    3. Yang J
    4. Liu S
    5. Zhu Z
    6. Lv J
    7. Chen G
    8. Jia C
    (2017)
    The Excavation of the Beiyuan Site at Shatang Town in Longde County, Ningxia in 2013
    Relics and Museolgy 6:3–12.
    1. Librado P
    2. Khan N
    3. Fages A
    4. Kusliy MA
    5. Suchan T
    6. Tonasso-Calvière L
    7. Schiavinato S
    8. Alioglu D
    9. Fromentier A
    10. Perdereau A
    11. Aury J-M
    12. Gaunitz C
    13. Chauvey L
    14. Seguin-Orlando A
    15. Der Sarkissian C
    16. Southon J
    17. Shapiro B
    18. Tishkin AA
    19. Kovalev AA
    20. Alquraishi S
    21. Alfarhan AH
    22. Al-Rasheid KAS
    23. Seregély T
    24. Klassen L
    25. Iversen R
    26. Bignon-Lau O
    27. Bodu P
    28. Olive M
    29. Castel J-C
    30. Boudadi-Maligne M
    31. Alvarez N
    32. Germonpré M
    33. Moskal-Del Hoyo M
    34. Wilczyński J
    35. Pospuła S
    36. Lasota-Kuś A
    37. Tunia K
    38. Nowak M
    39. Rannamäe E
    40. Saarma U
    41. Boeskorov G
    42. Lōugas L
    43. Kyselý R
    44. Peške L
    45. Bălășescu A
    46. Dumitrașcu V
    47. Dobrescu R
    48. Gerber D
    49. Kiss V
    50. Szécsényi-Nagy A
    51. Mende BG
    52. Gallina Z
    53. Somogyi K
    54. Kulcsár G
    55. Gál E
    56. Bendrey R
    57. Allentoft ME
    58. Sirbu G
    59. Dergachev V
    60. Shephard H
    61. Tomadini N
    62. Grouard S
    63. Kasparov A
    64. Basilyan AE
    65. Anisimov MA
    66. Nikolskiy PA
    67. Pavlova EY
    68. Pitulko V
    69. Brem G
    70. Wallner B
    71. Schwall C
    72. Keller M
    73. Kitagawa K
    74. Bessudnov AN
    75. Bessudnov A
    76. Taylor W
    77. Magail J
    78. Gantulga J-O
    79. Bayarsaikhan J
    80. Erdenebaatar D
    81. Tabaldiev K
    82. Mijiddorj E
    83. Boldgiv B
    84. Tsagaan T
    85. Pruvost M
    86. Olsen S
    87. Makarewicz CA
    88. Valenzuela Lamas S
    89. Albizuri Canadell S
    90. Nieto Espinet A
    91. Iborra MP
    92. Lira Garrido J
    93. Rodríguez González E
    94. Celestino S
    95. Olària C
    96. Arsuaga JL
    97. Kotova N
    98. Pryor A
    99. Crabtree P
    100. Zhumatayev R
    101. Toleubaev A
    102. Morgunova NL
    103. Kuznetsova T
    104. Lordkipanize D
    105. Marzullo M
    106. Prato O
    107. Bagnasco Gianni G
    108. Tecchiati U
    109. Clavel B
    110. Lepetz S
    111. Davoudi H
    112. Mashkour M
    113. Berezina NY
    114. Stockhammer PW
    115. Krause J
    116. Haak W
    117. Morales-Muñiz A
    118. Benecke N
    119. Hofreiter M
    120. Ludwig A
    121. Graphodatsky AS
    122. Peters J
    123. Kiryushin KY
    124. Iderkhangai T-O
    125. Bokovenko NA
    126. Vasiliev SK
    127. Seregin NN
    128. Chugunov KV
    129. Plasteeva NA
    130. Baryshnikov GF
    131. Petrova E
    132. Sablin M
    133. Ananyevskaya E
    134. Logvin A
    135. Shevnina I
    136. Logvin V
    137. Kalieva S
    138. Loman V
    139. Kukushkin I
    140. Merz I
    141. Merz V
    142. Sakenov S
    143. Varfolomeyev V
    144. Usmanova E
    145. Zaibert V
    146. Arbuckle B
    147. Belinskiy AB
    148. Kalmykov A
    149. Reinhold S
    150. Hansen S
    151. Yudin AI
    152. Vybornov AA
    153. Epimakhov A
    154. Berezina NS
    155. Roslyakova N
    156. Kosintsev PA
    157. Kuznetsov PF
    158. Anthony D
    159. Kroonen GJ
    160. Kristiansen K
    161. Wincker P
    162. Outram A
    163. Orlando L
    (2021) The origins and spread of domestic horses from the Western Eurasian steppes
    Nature 598:634–640.
    https://doi.org/10.1038/s41586-021-04018-9
  57. Book
    1. Taylor W
    2. Tuvshinjargal T
    (2018)
    Horseback riding, asymmetry, and anthropogenic changes to the equine skull: evidence for mounted riding in Mongolia’s late Bronze Age
    In: Bartosiewicz László, Gál Erika, editors. Care or Neglect?: Evidence of Animal Disease in Archaeology ; Proceedings of the 6th Meeting of the Animal Palaeopathology Working Group of the International Council for Archaeozoology (ICAZ). Oxbow Books. pp. 134–154.
    1. Wang Z
    2. Guo X
    3. Kang N
    4. Liu X
    5. Hu K
    6. Chen J
    (2015)
    Preliminary Report on the Excavation of Muzhuzhuliang Site in Shenmu, Shaanxi
    Archaeology and Cultural Relics 5:2015.
  64. Book
    1. Yuan J
    2. Flad R
    (2006)
    Research on Early Horse Domestication
    China: Oxbow Books.

Article and author information

Author details

  1. Dawei Cai

    Bioarchaeology Laboratory, Jilin University, Changchun, China
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing - original draft, Writing - review and editing
    Contributed equally with
    Siqi Zhu and Mian Gong
    For correspondence
    caidw@jlu.edu.cn
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6650-0217

  2. Siqi Zhu

    Bioarchaeology Laboratory, Jilin University, Changchun, China
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing - original draft, Writing - review and editing
    Contributed equally with
    Dawei Cai and Mian Gong
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6307-1189

  3. Mian Gong

    Key Laboratory of Animal Genetics, Breeding and Reproduction of Shaanxi Province, College of Animal Science and Technology, Northwest A&F University, Yangling, China
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing - original draft, Writing - review and editing
    Contributed equally with
    Dawei Cai and Siqi Zhu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1785-0621

  4. Naifan Zhang

    Bioarchaeology Laboratory, Jilin University, Changchun, China
    Contribution
    Investigation, Validation
    Competing interests
    No competing interests declared

  5. Jia Wen

    Key Laboratory of Animal Genetics, Breeding and Reproduction of Shaanxi Province, College of Animal Science and Technology, Northwest A&F University, Yangling, China
    Contribution
    Data curation, Formal analysis, Investigation, Software, Validation, Visualization
    Competing interests
    No competing interests declared

  6. Qiyao Liang

    Bioarchaeology Laboratory, Jilin University, Changchun, China
    Contribution
    Investigation, Validation, Visualization
    Competing interests
    No competing interests declared

  7. Weilu Sun

    Bioarchaeology Laboratory, Jilin University, Changchun, China
    Contribution
    Investigation, Validation
    Competing interests
    No competing interests declared

  8. Xinyue Shao

    Bioarchaeology Laboratory, Jilin University, Changchun, China
    Contribution
    Investigation, Validation
    Competing interests
    No competing interests declared

  9. Yaqi Guo

    Bioarchaeology Laboratory, Jilin University, Changchun, China
    Contribution
    Investigation, Validation
    Competing interests
    No competing interests declared

  10. Yudong Cai

    Key Laboratory of Animal Genetics, Breeding and Reproduction of Shaanxi Province, College of Animal Science and Technology, Northwest A&F University, Yangling, China
    Contribution
    Data curation, Investigation, Methodology, Software
    Competing interests
    No competing interests declared

  11. Zhuqing Zheng

    Key Laboratory of Animal Genetics, Breeding and Reproduction of Shaanxi Province, College of Animal Science and Technology, Northwest A&F University, Yangling, China
    Contribution
    Data curation, Investigation, Methodology, Software
    Competing interests
    No competing interests declared

  12. Wei Zhang

    Heilongjiang Provincial Institute of Cultural Relics and Archaeology, Harbin, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  13. Songmei Hu

    Shaanxi Provincial Institute of Archaeology, Xi’an, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  14. Xiaoyang Wang

    Ningxia Institute of Cultural Relics and Archaeology, Yinchuan, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  15. He Tian

    Heilongjiang Provincial Institute of Cultural Relics and Archaeology, Harbin, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  16. Youqian Li

    Heilongjiang Provincial Institute of Cultural Relics and Archaeology, Harbin, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  17. Wei Liu

    Heilongjiang Provincial Institute of Cultural Relics and Archaeology, Harbin, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  18. Miaomiao Yang

    Shaanxi Provincial Institute of Archaeology, Xi’an, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  19. Jian Yang

    Ningxia Institute of Cultural Relics and Archaeology, Yinchuan, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  20. Duo Wu

    College of Earth and Environmental Sciences; MOE Key Laboratory of Western China's Environmental Systems, Lanzhou University, Lanzhou, China
    Contribution
    Resources
    Competing interests
    No competing interests declared

  21. Ludovic Orlando

    Centre d’Anthropobiologie et de Génomique de Toulouse (CAGT), CNRS UMR 5288, Université de Toulouse, Université Paul Sabatier, Toulouse, France, Toulouse, France
    Contribution
    Conceptualization, Methodology, Supervision, Writing - review and editing
    For correspondence
    ludovic.orlando@univ-tlse3.fr
    Competing interests
    No competing interests declared

  22. Yu Jiang

    Key Laboratory of Animal Genetics, Breeding and Reproduction of Shaanxi Province, College of Animal Science and Technology, Northwest A&F University, Yangling, China
    Contribution
    Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing - original draft, Writing - review and editing
    For correspondence
    yu.jiang@nwafu.edu.cn
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4821-3585

Funding

Major Program of National Fund of Philosophy and Social Science of China (17ZDA221)

  • Dawei Cai

H2020 European Research Council (681605)

  • Ludovic Orlando

National Natural Science Foundation of China (31822052)

  • Yu Jiang

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank High-Performance Computing (HPC) of Northwest A&F University (NWAFU) for providing computing resources.

Copyright

© 2022, Cai, Zhu, Gong et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,610
    views
  • 410
    downloads
  • 10
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Dawei Cai
  2. Siqi Zhu
  3. Mian Gong
  4. Naifan Zhang
  5. Jia Wen
  6. Qiyao Liang
  7. Weilu Sun
  8. Xinyue Shao
  9. Yaqi Guo
  10. Yudong Cai
  11. Zhuqing Zheng
  12. Wei Zhang
  13. Songmei Hu
  14. Xiaoyang Wang
  15. He Tian
  16. Youqian Li
  17. Wei Liu
  18. Miaomiao Yang
  19. Jian Yang
  20. Duo Wu
  21. Ludovic Orlando
  22. Yu Jiang
(2022)
Radiocarbon and genomic evidence for the survival of Equus Sussemionus until the late Holocene
eLife 11:e73346.
https://doi.org/10.7554/eLife.73346

Share this article

https://doi.org/10.7554/eLife.73346

Categories and tags

Further reading

    1. Evolutionary Biology

    Fecal Microbiomes: Baboons, bacteria, and biological clocks address an age-old question

    Amanda D Melin
    Insight

    Studying the fecal microbiota of wild baboons helps provide new insight into the factors that influence biological aging.

    1. Cancer Biology
    2. Evolutionary Biology

    High-density sampling reveals volume growth in human tumours

    Arman Angaji, Michel Owusu ... Johannes Berg
    Research Article

    In growing cell populations such as tumours, mutations can serve as markers that allow tracking the past evolution from current samples. The genomic analyses of bulk samples and samples from multiple regions have shed light on the evolutionary forces acting on tumours. However, little is known empirically on the spatio-temporal dynamics of tumour evolution. Here, we leverage published data from resected hepatocellular carcinomas, each with several hundred samples taken in two and three dimensions. Using spatial metrics of evolution, we find that tumour cells grow predominantly uniformly within the tumour volume instead of at the surface. We determine how mutations and cells are dispersed throughout the tumour and how cell death contributes to the overall tumour growth. Our methods shed light on the early evolution of tumours in vivo and can be applied to high-resolution data in the emerging field of spatial biology.

    1. Evolutionary Biology

    Heterozygote advantage can explain the extraordinary diversity of immune genes

    Mattias Siljestam, Claus Rueffler
    Research Article

    The majority of highly polymorphic genes are related to immune functions and with over 100 alleles within a population, genes of the major histocompatibility complex (MHC) are the most polymorphic loci in vertebrates. How such extraordinary polymorphism arose and is maintained is controversial. One possibility is heterozygote advantage (HA), which can in principle maintain any number of alleles, but biologically explicit models based on this mechanism have so far failed to reliably predict the coexistence of significantly more than ten alleles. We here present an eco-evolutionary model showing that evolution can result in the emergence and maintenance of more than 100 alleles under HA if the following two assumptions are fulfilled: first, pathogens are lethal in the absence of an appropriate immune defence; second, the effect of pathogens depends on host condition, with hosts in poorer condition being affected more strongly. Thus, our results show that HA can be a more potent force in explaining the extraordinary polymorphism found at MHC loci than currently recognized.