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.

Decision letter

  1. George H Perry
    Senior and Reviewing Editor; Pennsylvania State University, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Ancient genomes redate the extinction of Sussemionus, a subgenus of Equus, to late Holocene" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and George Perry as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential Revisions include providing additional detail on existing datasets, the adjustment and re-processing of some analyses, and the addition of a section on archaeological context and expanded discussion on regional zooarchaeological implications of your findings. These points and more excellent feedback are detailed in the below reviews. In addition, I agree with reviewer #2s suggestion to revise your title for the next version of your paper.

Overall: Well done on this paper. We look forward to receiving and reviewing your revision. I will share it with the two reviewers to confirm that their comments have been suitably addressed.

Reviewer #2 (Recommendations for the authors):

I suggest the author's reconsider their manuscript title as it is the radiocarbon dates, rather than the ancient genomes, that re-date the extinction.

Given the specimens under study are isolated skeletal elements from a handful of stratigraphic horizons, it is possible that some of the specimens belonged to the same individual. There should be a description in the Methods about the minimum number of individuals at each locality. The mitogenomic data would prove valuable in this regard.

The authors use relative X-chromosome and autosomal coverage to distinguish male and female individuals (L104-106). However, these assignments should be justified with data. A new supplementary table with mean coverage of the autosomes and the X chromosome together with the ratio between the two would suffice.

There are issues with the BEAST-derived Bayesian phylogeny (Figure 2 —figure supplement 6). First, the ages of the tips do not appear to have been constrained with the known ages of the ancient specimens (or they are out by an order of magnitude). For example, JX312734 looks to be ~400,000 years old, when its age (based on Table S4) is 40,000 years old. Note also that ancient E. caballus individuals have been constrained as modern. Second, the constant population and strict clock models used (L538) are not suitable for the interspecies analysis used here. The authors should consider the Birth-death serially-sampled and relaxed clock models. Third, there is no description of how the molecular clock was calibrated (fossil, previous genomic estimate, and/or fixed mutation rate). I refer the authors to two publications for reference (dois: 10.1111/mec.15977, 10.7554/eLife.29944).

The authors present two D-statistics analyses based on alignment to either the horse or donkey reference genomes, which give very different results (Figure 3 —figure supplement 1 and 2). Presumably the larger D-stats when an African ass is included in the donkey reference analysis is an artifact of reference genome bias, and so the horse genome (outgroup) results should be considered more reliable. The authors should add a statement about this disagreement and an explanation in the Results, as this will not be clear to non-specialists, especially given the statement on L579-581 that the donkey reference genome analysis was included to 'eliminate the bias of the reference genome'.

It is stated that Z-scores of {greater than or equal to}3 were considered statistically significant on L577-579. However, the Z-score data are not presented and so it is not possible to determine which of the D-statistics in Figure 3 —figure supplement 1 and 2 are significant or not (L793, L801).

The G-PhoCS analysis suggests major introgression between the early non-caballine equid lineages. However, none of these events are recovered in the TreeMix analyses even with up to three migration edges considered. The apparently conflicting signal between these two analyses needs to be explained.

It will not be clear to non-specialists how the total migration rate for a single direction can be >100% (Tables S6 and S7). Please add a statement in the Methods as to why this occurs.

More details are needed in the G-PhoCs methods to enable reproducibility. Specifically, how were non-coding RNA genes identified and removed (L599-600), and what thresholds and methods were used to detect enriched misaligned bases and high false negative rates (L602-604)? In this regard, it may be helpful if the authors made their code available for these methods.

There is no discussion of the discordance between the mitochondrial and nuclear genome trees, which the G-PHoCS analysis seems to shed some light upon. I invite the authors to comment on this.

I almost missed that the authors have already made their raw sequence data publicly available (included in eLife manuscript information but not in the manuscript). To ensure readers can easily find the raw data, I suggest that the authors give a link to the European Nucleotide Archive BioProject code on L456.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Ancient DNA research redates the extinction of Sussemionus, a subgenus of Equus, to late Holocene" for further consideration by eLife. Your revised article has been evaluated by George Perry (Senior Editor) and the two reviewers of the previous version of your paper.

The manuscript has been improved but there are some remaining issues that need to be addressed. Detailed points of required revision are noted below, but in addition I will note that in my view your manuscript revisions in response to excellent and important points made by the reviewers in their original review are too superficial in multiple instances, e.g. in cases where expanded discussion or incorporation of particular concepts into your interpretation were requested, but the points were addressed with the addition of a short sentence only rather than taking the opportunity to maximally improve the manuscript, which is what we expect.

Thus, in your response, please detail the further revisions you made to the previous set of review comments, with the above in mind, in addition to further point-by-point responses to the specific comments below. This will be the final opportunity to revise your manuscript.

1. The issue with the title is not resolved.

2. Please provide expanded information on the Β Analytic 14C methods and results, if at all possible. (please contact the company for more specific information on how the samples were processed, for the sake of methodological completeness and data reproducibility).

3. The level of detail in the new archaeological background paragraph should be further improved. The revisions do accomplish the important goal of pointing the reader to the relevant background literature/citations, but acknowledgment and/or summary of the state of knowledge of the archaeological record of equids in the study region is incomplete. At the least this should include reference to Yuan and Flad's 2006 summary (Research on Early Horse Domestication in China. In Equids in Time and Space, ed. by Marjan Mashkour, pp. 124-131.)

4. The new sentence on 102-103 should be removed, and the sentence ending this paragraph on lines 104-105 needs to explain much more specifically what "no traces of domestication" means (e.g. no paleopathological problems?) and what "indicates they were hunted for food" (e.g. butchery patterns indicating meat removal? arrowheads imbedded in bone?).

5. It is unclear where the information about minimum number of individuals is coming from. The authors state that there were at least "31 individuals in the Honghe samples" (L460), yet there are only 20 samples from this locality in Supplementary File 1a. Further, the authors need to expand on the statement "the same process is repeated for the other two sites to ensure the specimens are unique" by including comparable counts.

6. We thank the authors for applying some of the suggested changes to the BEAST analysis. However, the tip dates for known sample ages have still not been constrained. Contrary to the rebuttal letter, there should not be any "deviations from the known ages" as these parameters should be fixed.

7. The explanation for the discordance between the G-PhoCS and TreeMix analyses needs to be stated in the manuscript.

8. G-PhoCS burn-in: although the MCMC run settings may use a burn-in of 0, the burn-in needs to be removed during post-processing (as is stated in the G-PhoCS user manual). For example, Vershinina et al. 2021 (doi: 10.1111/mec.15977) used a burn-in of 10%. This needs to be applied.

9. Figure 4 figure supplement 2: This is still confusing to the reader. The panel keys should be updated to reflect that all PSMC plots are based on E. a. somalicus.

10. Supplementary File 1c: correct 'gender' to 'sex'

11. Supplementary File 1e: correct 'yBC'. The age given for Haringtonhippus francisci is uncalibrated.

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

Author response

Reviewer #2 (Recommendations for the authors):

I suggest the author's reconsider their manuscript title as it is the radiocarbon dates, rather than the ancient genomes, that re-date the extinction.

Thank you for your suggestion, and we have revised the title accordingly: “Ancient DNA research redates the extinction of Sussemionus, a subgenus of Equus, to late Holocene”.

Given the specimens under study are isolated skeletal elements from a handful of stratigraphic horizons, it is possible that some of the specimens belonged to the same individual. There should be a description in the Methods about the minimum number of individuals at each locality. The mitogenomic data would prove valuable in this regard.

We thank the reviewer for pointing this out. The description about the minimum number of individuals was added in the Methods (lines 457-461): “Considering the preservation status and quantity, the minimum number of individuals was determined by assigning the frequency of hip bone and was calculated from the acetabular bone to avoid double-counting. Based on counts of skeletal elements, there is a minimum of 31 individuals in the Honghe samples. The same process is repeated for the other two sites to ensure the specimens are unique”. Meanwhile, the mitochondrial genome of each individual was visually checked to ensure they are unique according to your suggestion.

The authors use relative X-chromosome and autosomal coverage to distinguish male and female individuals (L104-106). However, these assignments should be justified with data. A new supplementary table with mean coverage of the autosomes and the X chromosome together with the ratio between the two would suffice.

We apologize for the lack of information here. We have now added Supplementary File 1c with mean coverage of the autosomes and the X chromosome together with the ratio between them.

There are issues with the BEAST-derived Bayesian phylogeny (Figure 2 —figure supplement 6). First, the ages of the tips do not appear to have been constrained with the known ages of the ancient specimens (or they are out by an order of magnitude). For example, JX312734 looks to be ~400,000 years old, when its age (based on Table S4) is 40,000 years old. Note also that ancient E. caballus individuals have been constrained as modern. Second, the constant population and strict clock models used (L538) are not suitable for the interspecies analysis used here. The authors should consider the Birth-death serially-sampled and relaxed clock models. Third, there is no description of how the molecular clock was calibrated (fossil, previous genomic estimate, and/or fixed mutation rate). I refer the authors to two publications for reference (dois: 10.1111/mec.15977, 10.7554/eLife.29944).

Thank you for this suggestion. First, we have set tip dates according to the ages of species in Supplementary File 1e. But considering an enormous range of time scales, few ancient specimens may showed some deviations from the known ages in Bayesian phylogeny.

Second, we have reconstructed Bayesian phylogeny using Birth-death model and relaxed molecular clock in Figure 2—figure supplement 4 according to your suggestion.

Third, we have added the description in lines 602-604: “we calibrated the tree using an age of 4–4.5 Mya for the root of crown group E. caballus (normal prior, mean 4.25 Mya, stdev: 0.15 Mya)” (see L. Orlando et al. (2013), https://www.nature.com/articles/nature12323).

The authors present two D-statistics analyses based on alignment to either the horse or donkey reference genomes, which give very different results (Figure 3 —figure supplement 1 and 2). Presumably the larger D-stats when an African ass is included in the donkey reference analysis is an artifact of reference genome bias, and so the horse genome (outgroup) results should be considered more reliable. The authors should add a statement about this disagreement and an explanation in the Results, as this will not be clear to non-specialists, especially given the statement on L579-581 that the donkey reference genome analysis was included to 'eliminate the bias of the reference genome'.

Thanks for spotting this. The similar statement was given on lines 647-649 in the manuscript, and we have added the sentence “Given the larger D-stats when an African ass is included in the donkey reference analysis is an artifact of reference genome bias, so that the horse reference genome results should be considered more reliable” in lines 649-652 according to the suggestion.

It is stated that Z-scores of {greater than or equal to}3 were considered statistically significant on L577-579. However, the Z-score data are not presented and so it is not possible to determine which of the D-statistics in Figure 3 —figure supplement 1 and 2 are significant or not (L793, L801).

We apologize for missing the legends to present the Z-score data. We have added the sentences “The nonsignificant results are shown in gray” in lines 887-888 and 895-896.

The G-PhoCS analysis suggests major introgression between the early non-caballine equid lineages. However, none of these events are recovered in the TreeMix analyses even with up to three migration edges considered. The apparently conflicting signal between these two analyses needs to be explained.

We thank the reviewer for the suggestion. Previous studies found that the TreeMix models will work best when gene flow between populations is restricted to a relatively short time period, situations of continuous migration violate this assumption and lead to unclear results (see Pickrell, Joseph K., and Pritchard, Jonathan K. (2012), https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002967). So compared with the G-PhoCS analysis, the TreeMix analyses had its limitation in an enormous range of time scales.

It will not be clear to non-specialists how the total migration rate for a single direction can be >100% (Tables S6 and S7). Please add a statement in the Methods as to why this occurs.

Thank you for this suggestion. We have now added a statement in the Methods. “The total migration rate gives an accumulated rate over a long period of time so that it can be >100%.” (lines 714-716)

More details are needed in the G-PhoCs methods to enable reproducibility. Specifically, how were non-coding RNA genes identified and removed (L599-600), and what thresholds and methods were used to detect enriched misaligned bases and high false negative rates (L602-604)? In this regard, it may be helpful if the authors made their code available for these methods.

Thanks for the suggestion. The non-coding RNA genes were identified using GFF annotation files.

We apologize for being unclear, and we have changed the sentence in lines 673-677: “Besides the various hard 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. So in this case, different individuals may be treated differently depending on the result of genotyping in 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”. Enriched misaligned bases and high false negative rates were embodied in (1) indels, (2) triallelic sites, (3) positions with depth of coverage twice the mean depth recorded for each individual, and; (4) transition sites.

There is no discussion of the discordance between the mitochondrial and nuclear genome trees, which the G-PHoCS analysis seems to shed some light upon. I invite the authors to comment on this.

This is certainly interesting suggestion. The discordance between the mitochondrial and nuclear genome trees can be caused by two reasons. First, mitochondrial DNA is maternally inherited and therefore variation in it will reflect disper-sal and history of the maternal lineage only. Second, two mitochondrial Maximum Likelihood trees based on all 6 partitions and excluding the control region were both reconstructed in Figure 2—figure supplement 5. It is not the latter but the former is discordant with nuclear genome trees, which may cause by exhibited significant increased damage in the mitochondrial control region.

I almost missed that the authors have already made their raw sequence data publicly available (included in eLife manuscript information but not in the manuscript). To ensure readers can easily find the raw data, I suggest that the authors give a link to the European Nucleotide Archive BioProject code on L456.

Thanks for the suggestion, and we have now given a link to the European Nucleotide Archive BioProject code in line 512.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed. Detailed points of required revision are noted below, but in addition I will note that in my view your manuscript revisions in response to excellent and important points made by the reviewers in their original review are too superficial in multiple instances, e.g. in cases where expanded discussion or incorporation of particular concepts into your interpretation were requested, but the points were addressed with the addition of a short sentence only rather than taking the opportunity to maximally improve the manuscript, which is what we expect.

Thus, in your response, please detail the further revisions you made to the previous set of review comments, with the above in mind, in addition to further point-by-point responses to the specific comments below. This will be the final opportunity to revise your manuscript.

1. The issue with the title is not resolved.

We have rephrased the title to indicate that it is the combination of both radiocarbon dating and phylogenomic that help reconsider the extinction/survival of Equus Sussemionus to the late Holocene. Our new title reads as follows:

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

2. Please provide expanded information on the Β Analytic 14C methods and results, if at all possible. (please contact the company for more specific information on how the samples were processed, for the sake of methodological completeness and data reproducibility).

We now provide the requested information as:

– A dedicated paragraph in the Methods section (page 29, lines 555-560: “Radiocarbon dating of the samples was performed at the Β Analytic Radiocarbon Dating Laboratory, Miami, Florida. Bone or tooth pieces about 2g 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.”);

– A sentence in the main text (lines 142-146): “Combined, these samples were radiocarbon dated to 3,456-4,460 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. 3,456-3,616 cal BP) (Supplementary File 1b).”

– A table referring to laboratory numbers, uncalibrated estimates and confidence range, and calendar years calibrated estimates (IntCal20) (see Supplementary File 1b).

3. The level of detail in the new archaeological background paragraph should be further improved. The revisions do accomplish the important goal of pointing the reader to the relevant background literature/citations, but acknowledgment and/or summary of the state of knowledge of the archaeological record of equids in the study region is incomplete. At the least this should include reference to Yuan and Flad's 2006 summary (Research on Early Horse Domestication in China. In Equids in Time and Space, ed. by Marjan Mashkour, pp. 124-131.)

We apologize for being unclear. Based on all excavated equine fossil bones found at Honghe, the Minimum Number of Individuals (NMI) was estimated to 31 individuals. And because of the preservation status, ancient DNA sequences were recovered from 20 of the 31 samples (Supplementary File 1a). This is now fully detailed at page 26 (lines 491-501): “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).”

4. The new sentence on 102-103 should be removed, and the sentence ending this paragraph on lines 104-105 needs to explain much more specifically what "no traces of domestication" means (e.g. no paleopathological problems?) and what "indicates they were hunted for food" (e.g. butchery patterns indicating meat removal? arrowheads imbedded in bone?).

We have now removed the sentence indicated, and have rephrase the following ones, appearing on lines 117-125: “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 these specimens were hunted.”

5. It is unclear where the information about minimum number of individuals is coming from. The authors state that there were at least "31 individuals in the Honghe samples" (L460), yet there are only 20 samples from this locality in Supplementary File 1a. Further, the authors need to expand on the statement "the same process is repeated for the other two sites to ensure the specimens are unique" by including comparable counts.

We apologize for being unclear. Based on all excavated equine fossil bones found at Honghe, the Minimum Number of Individuals (NMI) was estimated to 31 individuals. And because of the preservation status, ancient DNA sequences were recovered from 20 of the 31 samples (Supplementary File 1a). This is now fully detailed at page 26 (lines 491-501): “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).”

6. We thank the authors for applying some of the suggested changes to the BEAST analysis. However, the tip dates for known sample ages have still not been constrained. Contrary to the rebuttal letter, there should not be any "deviations from the known ages" as these parameters should be fixed.

We have proceeded according to the editor’s suggestion and have now used tip-calibrations (average calibrated radiocarbon dates) in our BEAST analyses (Supplementary File 1j). The full procedure is now described on lines 633-649, with the resulting tree shown on Figure 2—figure supplement 4.

7. The explanation for the discordance between the G-PhoCS and TreeMix analyses needs to be stated in the manuscript.

We have added the requested explanation in lines 787-805 (pages 39-40): “We caution that the analyses carried out using TreeMix and G-PhoCS returned partly discordant results. This may be due to TreeMix modelling 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 modelled 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.”

8. G-PhoCS burn-in: although the MCMC run settings may use a burn-in of 0, the burn-in needs to be removed during post-processing (as is stated in the G-PhoCS user manual). For example, Vershinina et al. 2021 (doi: 10.1111/mec.15977) used a burn-in of 10%. This needs to be applied.

We apologize for the unclear explanation. Although we run the MCMC with a pre-set burn-in of 0, the first 10% iterations were removed as burn-in during post-processing with Tracer v1.6 (http://tree.bio.ed.ac.uk/software/tracer/). This is shown in Author response image 1. Accordingly, we have edited the sentence in lines 742-743 (page 37): “The final results in Figure 3 are based on 500,000 MCMC iterations, considering the first 10% as burn-in.”

Author response image 1
Download asset Open asset

9. Figure 4 figure supplement 2: This is still confusing to the reader. The panel keys should be updated to reflect that all PSMC plots are based on E. a. somalicus.

We apologize for being unclear. The figure captions have now been rephrased to describe the procedure followed. It reads as follows (pages 60-61, lines 961-969): “Figure 4—figure supplement 2. Determining the uniform false-negative rate (uFNR) that was necessary for scaling PSMC results. (A) HH06D (11.30×), (B) KIA (10.68×) and (C) ONA (18.38×). The most suitable uFNR values for rescaling the PSMC profile are reported between squared brackets, to the right of the species names considered. The PSMC trajectory retrieved when considering all the sequence data available for the SOM individual is shown in blue. The green line provides the PSMC trajectory reconstructed when down-sampling these data to the average genome depth of coverage obtained for the species examined (top: Equus Sussemionus, red; centre: Equus kiang, purple, and; bottom: E. hemionus, purple).”

10. Supplementary File 1c: correct 'gender' to 'sex'

Done.

11. Supplementary File 1e: correct 'yBC'. The age given for Haringtonhippus francisci is uncalibrated.

Done.

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

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.

Senior and Reviewing Editor

  1. George H Perry, Pennsylvania State University, United States

Publication history

  1. Received: August 25, 2021
  2. Preprint posted: September 16, 2021 (view preprint)
  3. Accepted: May 11, 2022
  4. Accepted Manuscript published: May 11, 2022 (version 1)
  5. Version of Record published: May 27, 2022 (version 2)

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

  • 1,225
    Page views
  • 282
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

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

Categories and tags