Results and Discussion

Authenticity of aDNA Sequences.

We successfully recovered aDNA sequences from 22 of the 35 specimens (Fig. S1 and Table S1A): one horse (E. caballus from the Scladina cave, Belgium); 10 Hippidion specimens (five each from Patagonia and Peru); six E. (cf.) hydruntinus (two from Europe, three from Russia, and one from the Middle East); one E. hemionus (from the Middle East); and four E. capensis (from South Africa). Thirteen fossils did not yield any aDNA information despite extensive extraction and PCR attempts (Table S1B), and by default, also acted as cross-contamination controls in different extraction experiments. The poor DNA yields from these fossils are in agreement with their preservation state [high thermal age (26); see SI Text]. Thirteen of the samples were successfully amplified and sequenced in two independent laboratories [Lyon and Australian Centre for Ancient DNA (ACAD)], using different extraction methods, PCR conditions, and sequencing strategies (cloning/sequencing of different amplicons in Lyon versus direct sequencing at ACAD). It is notable that of the 6,569 bp duplicated at ACAD (from 64 amplicons), comparison with the final consensus sequences revealed only minor discrepancies (67 differences; average number of differences per base per fragment = 0.884%). As expected, a large majority of the variants (74.6%) consisted of GC→AT substitutions and are likely to derive from postmortem cytosine deamination characteristic of ancient templates (27, 28). The base recovered from the majority of independent amplicons was accepted as authentic, in agreement with standard recommendations (27, 29), and importantly, the extent of damage-induced errors was well below the intraspecific levels of genetic diversity among equids (Table S3). It is noteworthy that no modern equid DNA was stored in the aDNA laboratories, and though modern Grevy's zebra hairs were analyzed in the Lyon post-PCR laboratory, this haplotype was not found in any ancient sample. Of the new sequences, the mtDNA HVR-I sequences of the Patagonian Hippidion saldiasi specimens were highly similar to previously reported Hippidion sequences (15–17) (Fig. 1). The mtDNA HVR-I and/or cyt b sequences of the E. (cf.) hydruntinus from Scladina (Belgium) and Zagheh (Iran) showed high similarity with the two E. hydruntinus haplotypes reported in ref. 17 (Fig. 1 and Fig. S2 B and C), as did the new Iranian E. hemionus sample analyzed here (CH28). These are consistent with the ancient sequences being authentic.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Phylogenetic relationships among equids. Unrooted phylogenetic tree based on maximum likelihood analysis of the HVR-1 dataset. Sequences generated in this study are highlighted in gray. Numbers above/below branches correspond to node bootstrap values, approximate likelihood ratio test (aLRT SH-like), and Bayesian posterior probabilities. Δ28: 28-bp deletion in the mitochondrial HVR-1 (positions 15533–15560 according to the complete horse mitochondrial genome; GenBank accession no. X79547) (Fig. 2A). The potential locations of the root are reported with gray (alternative)/black (most likely) triangles according to likelihood-based topological tests (Fig. S4 and Table S5) and previous reports (15, 16). C1 and C2 nodes have been used as calibration points for molecular dating. Two radiocarbon dates were retrieved at the University of Colorado Laboratory for AMS Radiocarbon Preparation and Research (CURL) and the University of California-Irvine Accelerator Mass Spectrometry (UCIAMS) Facility for sussemione samples ACAD2303 and ACAD2305 (see SI Text).

Phylogenetic Analyses.

South America.

Both the maximum likelihood and Bayesian phylogenetic analyses (Fig. 1) of the mtDNA HVR sequences identified three monophyletic New World groups: caballine horses (node G); NWSL horses (node F); and hippidiforms (node C2), in agreement with earlier studies (15, 16, 18). The exact branching order of the three groups (relative to the Old World noncaballine horses) could not be resolved with either the HVR-1 dataset alone (node B shows weak supports in Fig. 1; Table S5, topologies L–Q) or various combinations of HVR-1 plus cyt b/HVR-2 data (Table S5, topologies R–V). Similarly, although the strongest support was obtained for the placement of hippidiforms within a larger paraphyletic group of Equus (with merged datasets, Table S2, and with marginal support from SH tests: P values = 0.065–0.083; Table S5, topologies R–V), likelihood-based topological tests cannot significantly reject alternative topologies and the exact position of the root, previously reported between nodes A and H (Fig. 1) (15, 16) should be considered preliminary. The lack of resolution is complicated by the short divergence time among caballines and New World horses (circa 0.5 MY; nodes A and B/B1/B2; Table S4) and the lack of a close outgroup, as has been noted with mammoths (30, 31). When the rhino was used as an outgroup, the data were RY coded to reduce possible mutation saturation artifacts resulting from this the deep divergence (55 MYA), but this removed support for most nodes (Table S2).

According to our molecular dating estimates, the different equid lineages (hippidiforms, NWSL, caballines, and noncaballines) originated 3.7–4.3 MYA (95% confidence range: 2.8–6.2 MYA; Table S4). This directly contrasts with classical palaeontological models of hippidiform origins as descendants of the Pliohippines (divergence time with the Equus lineage >10 MYA) (19) or as a lineage diverging from a (Dinohippus, Astrohippus, and Equus) clade ≈7–8 MYA (32), and considerably reduces the time gap between the supposed divergence of the hippidiform lineages and their first appearance in the fossil record 2.5 MYA (20). A deep split within hippidiforms separates one cluster (Fig. 1, node D) from Patagonia and Argentina (Hippidion saldiasi and Hippidion principale, respectively) from a second (Fig. 1, node E) consisting of specimens from a high-altitude cave in Peru. The genetic differentiation between these clusters is greater than intraspecific levels of genetic diversity in any current equid species (Fig. 2 and Table S3) but comparable to that between indisputable species from the same genus (Equus), such as hemiones and donkeys (i.e., two noncaballine Old World equids), and less than that between hippidions and different Equus species (such as caballines horses or donkeys). We therefore support the Peruvian specimens as a distinct species within the same genus as hippidiforms from Patagonian/Chile. The extinct Peruvian taxa is referable to Onohippidium devillei according to the mean length and width of metatarsals, which are similar to the equids from the type locality of Tarija, Bolivia (33), and should be renamed Hippidion devillei. Our analysis therefore shows that Onohippidium is a junior synonym genus of Hippidion in agreement with Alberdi and Prado (34, 35) and in contrast to a reanalysis of facial characters, dental patterns, and metatarsal proportions (19). The molecular dating estimates (Table S4) suggest that the hippidiforms (from South America) separated from the North American taxa (caballines and NWSL) after 3.2–3.6 MYA, consistent with the final formation of the Panamanian isthmus ≈3.0–3.7 MYA (36). Though we have only limited sample coverage, especially in Peru, we speculate that after this first vicariant event, the Peruvian Hippidion species may have specialized in higher altitude Andean habitats (almost all specimens from H. devillei come from high altitude environments) (37), and was genetically isolated from a more generalist species (H. saldiasi/principale) that extended from Brazil to Patagonia and showed large size plasticity (the adult weight of H. principale was twice that of H. saldiasi) (see refs. 35–37).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Intra- and interspecific corrected (GTR+Γ+I) pairwise distances among equids. The dataset consisted of 1,544 sequences encompassing positions 15518 and 15818 from the complete horse mitochondrial genome (GenBank accession no. X79547). Numeric values are reported on Table S3A. A complete list of accession numbers used for defining each taxonomic group is provided as SI Text. Parameters and distances were estimated using PhyML 3.0 and PAUP* 4.0, respectively. Bar, median; box, 50% quantiles; bars, 75% quantiles; CAB, E. caballus; ASI, E. asinus; QUA, E. quagga (all Plain zebras, including E. quagga quagga); QUA_J2, node J2 is defined on Fig. 1 and includes samples ACAD226 and ACAD230; QUA_ALL-J2, QUA but excluding samples for node J2; GRE, E. grevyi; ZEB, E. hartmannae; HEM, E. hemionus; KIA, E. kiang; HIP, Hippidion saldiasi/principale; ONO, Hippidion devillei (Peruvian hippidions; ACAD3615, ACAD3625, ACAD3627, ACAD3628, and ACAD3629); STL, New World Stilt-Legged horses; SUS, Sussemiones (Khakassia, SW Siberia; ACAD2302, ACAD2303, and ACAD2305). Only two distributions (GRE/ZEB and HEM/KIA) are shown for the intergroup genetic distances, given that all others are included within these two.

Africa.

Four specimens identified as giant Cape zebra (E. capensis) were found to have good-quality DNA preservation. Three specimens were from Wonderwerk Cave, south of Kuruman in central southern Africa, and one from Glen Craig Shelter near Port Elizabeth on the southeast coast. Fragments up to 477-bp long could be amplified, and the six target HVR-1 and cyt b gene fragments were characterized (687 bp in total) apart from one Wonderwerk specimen (ACAD227), which yielded only the barcode fragment (see below and Table S1A). Maximum likelihood and Bayesian analyses placed the sequences amongst the plains zebra group (Fig. 1, nodes C1 and J; Fig. S2 B and C), which consists of six subspecies that lack population genetic structure (9) and includes the extinct quagga (8) (E. q. quagga). The Cape zebra sequences clustered amongst the two southern subspecies (E. q. quagga, E. q. burchellii), with the Wonderwerk specimens ACAD226 and ACAD230 exhibiting haplotypes outside the previously known variation (node J2), but inside that of the northern subspecies (E. q. boehmi, node K, with the exception of specimen AF220918). Although only limited sequence information was retrieved (Table S1A), sample ACAD227 showed minimum genetic distance with node K E. q. boehmi specimens (i.e., over 87 nucleotides, one substitution is observed versus four or more for J1/J2 haplotypes). In addition, the Glen Craig specimen (ACAD236) grouped within the southern clade, suggesting limited reproductive isolation (if any) amongst Cape quaggas and plains zebras, or potential issues with the morphological identification of the fossil material. In support of the latter possibility, radiocarbon dating of one Wonderwerk specimen (ACAD230) yielded a recent date (cal. 148 ± 120 BP) in contrast to an expected Pleistocene age. Though further specimens and dates are obviously required, all specimens morphologically identified as Cape zebras yielded sequences that grouped phylogenetically within plains zebras. The genetic distance between the giant Cape zebra and extant plains zebras is within the range of the intraspecific diversity found within other zebra species (11) (Fig. 2 and Table S3), and so we do not consider separate species (or subspecies) status is warranted (until independent information from nuclear loci is available). In contrast, Churcher (24) used the large size of giant Cape zebra fossils (estimated 400 kg and 150 cm height at the withers) to support conspecific status with Grevy's zebra, the largest extant wild equid. However, skulls are considered to be the best taxonomic indicators for equids at the morphological level (38), and the skull proportions of a giant Cape zebra from Elandsfontein (South Africa) showed a close affinity with true Cape quaggas (23), in accordance with our genetic findings. Together, these findings suggest that the giant Cape zebra formed part of the same diverse taxonomic group as plains zebras and quaggas (E. q. quagga), and that this group possessed a marked plasticity in both morphological size (reminiscent of Hippidion) and coat color (8).

Eurasia.

We analyzed two Middle Eastern hemiones from Iran (CH28, CH30) and six Equus (cf.) hydruntinus samples (one Middle Eastern, TZ9; two European, CH561 and CH1069; and three from Khakassia, Russia, southwest Siberia) to clarify the taxonomic relationships and geographic range of these groups. Surprisingly, the Siberian specimens did not cluster with E. hydruntinus nor hemiones but formed a new monophyletic group with strong support (Fig. 1, node I) and no specific affinity with other noncaballine horses (Table S5, topologies J–K4). Cyt b sequences further support this arrangement (Fig. S2 B and C). Interestingly, this new clade does not possess a 28-bp deletion in the first HVR-1 PCR fragment that unite E. hydruntinus specimens from Belgium (CH561) and Crimea (17), and the other Asian and Middle Eastern hemiones [Figs. 1 (CH28) and 2A] relative to other equids. The genetic distances separating the three Russian specimens from hemiones are of a similar magnitude to other species (e.g., those separating donkey and hemione; Fig. 2 and Table S3), and together this evidence suggests that clade I should be considered a distinct taxonomic group, representing a new species, distinct from E. hydruntinus, with no extant relative. Isolated bones that are morphologically similar to the three Khakassian specimens have been recovered from cave sites in the Altai Mountains (Denisova, Kaminnaya, Logovo Gieny, Okladnikov, and Strashnaya) as well as from alluvial sites in the Altai piedmont plain (39). In contrast with our molecular findings, previous morphological analyses suggested similar metapodial proportions to classical E. hydruntinus from Western Europe, the Crimea, the Middle East, and Caucasus, although the former was larger and more robust. Our analysis of tooth morphological data (Fig. S4) shows more consistency with the molecular phylogeny in that there are poor affinities with hemiones and E. hydruntinus, and instead there are strong similarities with the archaic Middle Pleistocene Sussemiones from Süssenborn, Germany (40). The age to the most recent common ancestor of the three samples analyzed here provided a conservative estimate of 1.1–1.3 MYA (95% confidence range: 0.4–2.3 MYA; Table S4) for the emergence of this new species. Sussemiones possibly differentiated much earlier (circa 2 MYA) in Beringia, as suggested by the lower cheek series found at Lost Chicken, Alaska, and were extremely successful and dispersed across all Eurasia and Africa (41). Interestingly, Sussemiones were supposedly extinct well before the Late Pleistocene (40), but radiocarbon dating of two Khakassian specimens revealed that the Siberian species survived at least until circa 45–50 KYBP (SI Text).

To further examine the extinction date of E. hydruntinus, we examined a specimen from 17th-century Portugal (CH1069), purported to have been the last representative of this species (42). However, HVR-1 (median spanning network; Fig. S2 B and C), cyt b (Fig. S2 D and E), and cranial morphological data (Fig. S4G) indicate that this specimen showed a strong relationship to donkeys. This implies that E. hydruntinus, which first appeared in the fossil record 350 KYA in Lunel-Viel (43) (southern France), did not survive until the Middle Ages but most probably became extinct in the Iron Age (44). Interestingly, the genetic distance between E. hydruntinus and E. hemionus is within the diversity observed within hemiones, in agreement with previous studies of short (144–288 bp) HVR-1 sequences (17). Although the combination of the relative average dimensions of the metapodials and upper-cheek teeth parameters of E. hydruntinus differed from any other hemione (45), the genetic data suggests that E. hydruntinus was a subspecies of E. hemionus. However, it is worth noting that the kiang and onager/kulan are currently classified as separate species based on differences in coat color, morphology, geographic distribution, and the number of chromosomes, and yet show poor mitochondrial differentiation (7, 8, 46). Therefore, the data are still compatible with E. hydruntinus being a separate species, which perhaps experienced multiple introgression events from different stocks (from onagers in Zagheh, Iran, sample TZ9; and from kulan in Sagzabad, Iran, sequence DQ464012 reported in ref. 17). Consequently, we suggest that the taxonomic status of Equus hydruntinus as a separate species or a subspecies of E. hemionus should be clarified with genetic information from nuclear loci.

Genetic Taxonomy of Late Pleistocene and Holocene Equids.

The short (<89 bp) first fragment of HVR-1 appears to be a potentially useful barcode for the higher-level taxonomic identification of fossil equids (Fig. S2A). This fragment can be amplified in caballine horses, zebras, asses, and hemiones using our PCR conditions, and, importantly, contains the informative 28-bp deletion characteristic of the (E. hydruntinus/E. hemionus/E. kiang) lineage. It shows higher genetic distances between, rather than within, extant taxonomic groups (Table S3C) and could potentially be used to assess the taxonomic status of the many paleontological equid species. For instance, it could be used to test the suggestion that the many species of Late Pleistocene New World hemione-like equids (e.g., E. francisci, E. tau, E. quinni, E. cf. hemionus, E. cf. kiang) are conspecific within a single species of NWSL horse (16). This rapid assessment method could be a useful tool to disentangle the temporal and regional body-size variation and marked anatomical convergence observed within Late Pleistocene and Holocene equids. Of course, an mtDNA HVR genetic survey would be only a first step in identifying major taxonomic groups, which would need to be examined with morphological and other molecular studies. Indeed, the relatively short sequences we have generated cannot clearly resolve the exact branching order of major taxonomic groups and the bootstrap values, posterior probabilities, and approximate likelihood ratio tests (aLRT, SH-like) for these nodes were generally not strong. Recently, the use of whole mitochondrial genome sequences from Pleistocene fossils has helped resolve the rapid radiation events within elephantids (31) and ursids (47). Such a mitogenomic approach is likely to uncover the phylogenetic relationships within caballines and New World horses (Fig. 1, nodes A–G), and potentially within extant noncaballine horses as well, which are still poorly understood at the molecular level (Fig. 1, node H). For example, in Fig. 1, Grevyi's zebras nest within paraphyletic hemionids, and mountain zebras do not cluster with plains zebras but with donkeys, in striking contrast with traditionally accepted taxonomy based on coat patterns, behavior, morphology, geographic separation (6, 7, 10, 48), and the 28-bp deletion (Fig. 1) (17). However, current datasets provide only moderate (if any) bootstrap and aLRT supports and posterior probabilities (Table S2) for these branches, and they should be regarded with caution. Besides, likelihood-based topogical tests show significant support for the standard topology (in accordance with the current taxonomy) (KH and SH tests P values ≤0.003 for model E7, and KH test P value = 0.008 for model E7 vs. E6; Table S5). It seems likely that a (mito)genomic approach will be needed to solve major phylogenetic splits within equids and rapid radiation events (Table S4).

In this study, two species of extinct equids have been identified by mtDNA. The first is related to hippidiforms and corresponds to the paleontological genus Onohippidium, which should be reclassified within the Hippidion genus as Hippidion devillei. The second is a unique basal lineage of Old World equid, which appears morphologically related to the Middle Pleistocene Sussemiones, yet survived in southwestern Siberia at least until circa 45–50 KYBP. In addition, we propose to synonymize a variety of other species [Cape zebras, quaggas, and plains zebras (11, 13) and potentially E. hydruntinus and E. hemionus]. This reinforces similar suggestions for hippidiforms, NWSL, and caballids (15, 16, 18) and E. hydruntinus (17). This pattern of taxonomic oversplitting does not appear to be restricted to equids but is widespread amongst other Quaternary megafauna [e.g., Late Pleistocene bison (49); Holarctic cave lions (50); New World brown bears (51), and ratite moas (52, 53)]. Together, these findings suggest that the morphological plasticity of large terrestrial vertebrates across space and time has generally been underestimated, opening the way to detailed studies of the environmental, ecological, and epigenetic factors involved. Interestingly, in this regard the human lineage shows a rich fossil record over the last 6 MY, spreading over seven possible genera and 22 species (54). The exact number of taxonomic groups that should be recognized is still debated, even within our own genus (55), and in this context it is pertinent to consider the degree of taxonomic oversplitting, from species to generic levels, that aDNA has revealed amongst Late Pleistocene equids and other megafauna. A further important implication of this finding is that the number of megafaunal extinctions and loss of taxonomic diversity from the Pleistocene to modern day may not have been nearly as large as previously thought, at least at the species or subspecies level. Conversely, at the molecular level, aDNA studies on a wide range of large mammal taxa (49, 50, 56, 57) have revealed that the loss of genetic diversity over this time period has been much larger than previously recognized with major implications for the conservation biology of surviving populations (58).