A new genus of horse from Pleistocene North America.

The extinct 'New World stilt-legged', or NWSL, equids constitute a perplexing group of Pleistocene horses endemic to North America. Their slender distal limb bones resemble those of Asiatic asses, such as the Persian onager. Previous palaeogenetic studies, however, have suggested a closer relationship to caballine horses than to Asiatic asses. Here, we report complete mitochondrial and partial nuclear genomes from NWSL equids from across their geographic range. Although multiple NWSL equid species have been named, our palaeogenomic and morphometric analyses support the idea that there was only a single species of middle to late Pleistocene NWSL equid, and demonstrate that it falls outside of crown group Equus. We therefore propose a new genus, Haringtonhippus, for the sole species H. francisci. Our combined genomic and phenomic approach to resolving the systematics of extinct megafauna will allow for an improved understanding of the full extent of the terminal Pleistocene extinction event.

Introduction

The family that includes modern horses, asses, and zebras, the Equidae, is a classic model of macroevolution. The excellent fossil record of this family clearly documents its ~55 million year evolution from dog-sized hyracotheres through many intermediate forms and extinct offshoots to present-day Equus, which comprises all living equid species (MacFadden, 1992). The downside of this excellent fossil record is that many dubious fossil equid taxa have been erected, a problem especially acute within Pleistocene Equus of North America (Macdonald et al., 1992). While numerous species are described from the fossil record, molecular data suggest that most belonged to, or were closely related to, a single, highly variable stout-legged caballine species that includes the domestic horse, E. caballus (Weinstock et al., 2005). The enigmatic and extinct ‘New World stilt-legged’ (NWSL) forms, however, exhibit a perplexing mix of morphological characters, including slender, stilt-like distal limb bones with narrow hooves reminiscent of extant Eurasian hemionines, the Asiatic wild asses (E. hemionus, E. kiang) (Eisenmann, 1992; Eisenmann et al., 2008; Harington and Clulow, 1973; Lundelius and Stevens, 1970; Scott, 2004), and dentitions that have been interpreted as more consistent with either caballine horses (Lundelius and Stevens, 1970) or hemionines (MacFadden, 1992).

On the basis of their slender distal limb bones, the NWSL equids have traditionally been considered as allied to hemionines (e.g. Eisenmann et al., 2008; Guthrie, 2003; Scott, 2004; Skinner and Hibbard, 1972). Palaeogenetic analyses based on mitochondrial DNA (mtDNA) have, however, consistently placed NWSL equids closer to caballine horses (Barrón-Ortiz et al., 2017; Der Sarkissian et al., 2015; Orlando et al., 2008, 2009; Vilstrup et al., 2013; Weinstock et al., 2005). The current mtDNA-based phylogenetic model therefore suggests that the stilt-legged morphology arose independently in the New and Old Worlds (Weinstock et al., 2005) and may represent convergent adaptations to arid climates and habitats (Eisenmann, 1985). However, these models have been based on two questionable sources. The first is based on 15 short control region sequences (<1000 base pairs, bp; Barrón-Ortiz et al., 2017; Weinstock et al., 2005), a data type that can be unreliable for resolving the placement of major equid groups (Der Sarkissian et al., 2015; Orlando et al., 2009). The second consist of two mitochondrial genome sequences (Vilstrup et al., 2013) that are either incomplete or otherwise problematic (see Results). Given continuing uncertainty regarding the phylogenetic placement of NWSL equids—which impedes our understanding of Pleistocene equid evolution in general—we therefore sought to resolve their position using multiple mitochondrial and partial nuclear genomes from specimens representing as many parts of late Pleistocene North America as possible.

The earliest recognized NWSL equid fossils date to the late Pliocene/early Pleistocene (~2–3 million years ago, Ma) of New Mexico (Azzaroli and Voorhies, 1993; Eisenmann, 2003; Eisenmann et al., 2008). Middle and late Pleistocene forms tended to be smaller in stature than their early Pleistocene kin, and ranged across southern and extreme northwestern North America (i.e. eastern Beringia, which includes Alaska, USA and Yukon Territory, Canada). NWSL equids have been assigned to several named species, such as E. conversidens Owen 1869, E. tau Owen 1869, E. francisci Hay (1915), E. calobatus Troxell 1915, and E. (Asinus) cf. kiang, but there is considerable confusion and disagreement regarding their taxonomy. Consequently, some researchers have chosen to refer to them collectively as Equus (Hemionus) spp. (Guthrie, 2003; Scott, 2004), or avoid a formal taxonomic designation altogether (Der Sarkissian et al., 2015; Vilstrup et al., 2013; Weinstock et al., 2005). Using our phylogenetic framework and comparisons between specimens identified by palaeogenomics and/or morphology, we attempted to determine the taxonomy of middle-late Pleistocene NWSL equids.

Radiocarbon (14C) dates from Gypsum Cave, Nevada, confirm that NWSL equids persisted in areas south of the continental ice sheets during the last glacial maximum (LGM; ~26–19 thousand years before present (ka BP); Clark et al., 2009) until near the terminal Pleistocene, ~13 thousand radiocarbon years before present (14C ka BP) (Weinstock et al., 2005), soon after which they became extinct, along with their caballine counterparts and most other coeval species of megafauna (Koch and Barnosky, 2006). This contrasts with dates from unglaciated eastern Beringia, where NWSL equids were seemingly extirpated locally during a relatively mild interstadial interval centered on ~31 14C ka BP (Guthrie, 2003), thus prior to the LGM (Clark et al., 2009), final loss of caballine horses (Guthrie, 2003; 2006), and arrival of humans in the region (Guthrie, 2006). The apparently discrepant extirpation chronology between NWSL equids south and north of the continental ice sheets implies that their populations responded variably to demographic pressures in different parts of their range, which is consistent with results from some other megafauna (Guthrie, 2006; Zazula et al., 2014; Zazula et al., 2017). To further test this extinction chronology, we generated new radiocarbon dates from eastern Beringian NWSL equids.

We analyzed 26 full mitochondrial genomes and 17 partial nuclear genomes from late Pleistocene NWSL equids, which revealed that individuals from both eastern Beringia and southern North America form a single well-supported clade that falls outside the diversity of Equus and diverged from the lineage leading to Equus during the latest Miocene or early Pliocene. This novel and robust phylogenetic placement warrants the recognition of NWSL equids as a distinct genus, which we here name Haringtonhippus. After reviewing potential species names and conducting morphometric and anatomical comparisons, we determined that, based on the earliest-described specimen bearing diagnosable features, francisci Hay is the most well-supported species name. We therefore refer the analyzed NWSL equid specimens to H. francisci. New radiocarbon dates revealed that H. francisci was extirpated in eastern Beringia ~14 14C ka BP. In light of our analyses, we review the Plio-Pleistocene evolutionary history of equids, and the implications for the systematics of equids and other Pleistocene megafauna.

Results

Phylogeny of North American late Pleistocene and extant equids

We reconstructed whole mitochondrial genomes from 26 NWSL equids and four New World caballine Equus (two E. lambei, two E. cf. scotti). Using these and mitochondrial genomes of representatives from all extant and several late Pleistocene equids, we estimated a mitochondrial phylogeny, using a variety of outgroups (Appendix 1, Appendix 2—tables 1–2, and Supplementary file 1). The resulting phylogeny is mostly consistent with previous studies (Der Sarkissian et al., 2015; Vilstrup et al., 2013), including confirmation of NWSL equid monophyly (Weinstock et al., 2005). However, we recover a strongly supported placement of the NWSL equid clade outside of crown group diversity (Equus), but closer to Equus than to Hippidion (Figure 1, Figure 1—figure supplement 1a, Figure 1—source data 1, and Appendix 2—tables 1–2). In contrast, previous palaeogenetic studies placed the NWSL equids within crown group Equus, closer to caballine horses than to non-caballine asses and zebras (Barrón-Ortiz et al., 2017; Der Sarkissian et al., 2015; Orlando et al., 2008, 2009; Vilstrup et al., 2013; Weinstock et al., 2005). To explore possible causes for this discrepancy, we reconstructed mitochondrial genomes from previously sequenced NWSL equid specimens and used a maximum likelihood evolutionary placement algorithm (Berger et al., 2011) to place these published sequences in our phylogeny a posteriori. These analyses suggested that previous results were likely due to a combination of outgroup choice and the use of short, incomplete, or problematic mtDNA sequences (Appendix 2 and Appendix 2—table 3).

Appendix 2—table 1.

Topological shape and support values for the best supported trees.

These results are from the Bayesian and maximum likelihood (ML) analyses of mtDNA data sets 1–3, including either the all or reduced partition sets, and with Hippidion sequences either included or excluded. Topology numbers and node letters refer to those outlined in Appendix 2—figure 3. Bayesian posterior probability support of >0.99 and ML bootstrap support of >95% are in bold for nodes A and B. *support for nodes that are consistent with topology one in Appendix 2—figure 3. NCs: non-caballines.

Outgroup	Partitions	Hippidion?	Tips	Analysis method	Topology	Support
						Node A	Node B	Hippidion	NWSL	NCs	Caballines
White rhino  (Data set 1)	All	Excluded	63	Bayesian	1/2/3	0.996*	N/A	N/A	1.000	1.000	1.000
				ML	1/2/3	71*	N/A	N/A	100	99	100
		Included	69	Bayesian	2	0.751	1.000*	1.000	1.000	1.000	1.000
				ML	1	64*	96*	100	100	100	100
	Reduced	Excluded	63	Bayesian	1/2/3	1.000*	N/A	N/A	1.000	1.000	1.000
				ML	1/2/3	100*	N/A	N/A	99	100	100
		Included	69	Bayesian	2	0.948	1.000*	1.000	1.000	1.000	1.000
				ML	2	73	98*	100	99	100	100
Malayan tapir  (Data set 2)	All	Excluded	63	Bayesian	5/7	0.971	N/A	N/A	1.000	1.000	1.000
				ML	5/7	87	N/A	N/A	100	99	99
		Included	69	Bayesian	6	0.808	0.867	1.000	1.000	1.000	1.000
				ML	6	55	63	100	100	100	100
	Reduced	Excluded	63	Bayesian	1/2/3	0.675*	N/A	N/A	1.000	1.000	1.000
				ML	4/6	28	N/A	N/A	100	96	98
		Included	69	Bayesian	3	0.685	0.864*	1.000	1.000	1.000	1.000
				ML	3	70	69	100	100	100	100
Dog + ceratomorphs  (Data set 3)	All	Excluded	71	Bayesian	1/2/3	0.598*	N/A	N/A	1.000	1.000	1.000
				ML	4/6	59	N/A	N/A	100	100	100
		Included	77	Bayesian	1	1.000*	1.000*	1.000	1.000	1.000	1.000
				ML	1	94*	96*	100	100	100	100
	Reduced	Excluded	71	Bayesian	1/2/3	0.999*	N/A	N/A	1.000	1.000	1.000
				ML	1/2/3	97*	N/A	N/A	100	100	100
		Included	77	Bayesian	1	1.000*	1.000*	1.000	1.000	1.000	1.000
				ML	1	99*	100*	100	100	100	100

Appendix 2—table 2.

The a posteriori phylogenetic placement likelihood for eight ceratomorph (rhino and tapir) outgroups.

These analyses used a ML evolutionary placement algorithm, whilst varying the partition set used (all or reduced), and either including or excluding Hippidion sequences. Likelihoods >0.95 are in bold. Topology numbers refer to those outlined in Appendix 2—figure 3. Genbank accession numbers are given in parentheses after outgroup names.

Partitions	Outgroup	Hippidion?	Included				Excluded
		Topology	1	2	3	6	1/2/3	4/6	5/7
All	Tapirus terrestris (AJ428947)		0.456	0.317	0.205	0.018	0.549	0.313	0.139
	Tapirus indicus (NC023838)		0.275	0.105	0.225	0.389	0.050	0.908	0.042
	Coelodonta antiquitatis (NC012681)		0.998				0.248	0.451	0.301
	Dicerorhinus sumatrensis (NC012684)		0.981		0.009		0.155	0.553	0.292
	Rhinoceros unicornis (NC001779)		0.998				0.529	0.334	0.137
	Rhinoceros sondaicus (NC012683)		0.989	0.006			0.732	0.196	0.072
	Ceratotherium simum (NC001808)		0.448	0.499	0.053		0.949	0.018	0.033
	Diceros bicornis (NC012682)		0.917	0.065	0.018		0.851	0.073	0.076
Reduced	Tapirus terrestris (AJ428947)		0.410	0.391	0.199		0.987		0.012
	Tapirus indicus (NC023838)		0.536	0.298	0.166		0.995
	Coelodonta antiquitatis (NC012681)		0.411	0.554	0.035		1.000
	Dicerorhinus sumatrensis (NC012684)		0.983	0.015			1.000
	Rhinoceros unicornis (NC001779)		0.998				1.000
	Rhinoceros sondaicus (NC012683)		0.895	0.102			1.000
	Ceratotherium simum (NC001808)		0.296	0.704			1.000
	Diceros bicornis (NC012682)		0.996				1.000

Figure 1.

Phylogeny of extant and middle-late Pleistocene equids, as inferred from the Bayesian analysis of full mitochondrial genomes.

Purple node-bars illustrate the 95% highest posterior density of node heights and are shown for nodes with >0.99 posterior probability support. The range of divergence estimates derived from our nuclear genomic analyses is shown by the thicker, lime green node-bars ([Orlando et al., 2013]; this study). Nodes highlighted in the main text are labeled with boxed numbers. All analyses were calibrated using as prior information a caballine/non-caballine Equus divergence estimate of 4.0–4.5 Ma (Orlando et al., 2013) at node 3, and, in the mitochondrial analyses, the known ages of included ancient specimens. The thicknesses of nodes 2 and 3 represent the range between the median nuclear and mitochondrial genomic divergence estimates. Branches are coloured based on species provenance and the most parsimonious biogeographic scenario given the data, with gray indicating ambiguity. Fossil record occurrences for major represented groups (including South American Hippidion, New World stilt-legged equids, and Old World Sussemiones) are represented by the geographically coloured bars, with fade indicating uncertainty in the first appearance datum (after (Eisenmann et al., 2008; Forsten, 1992; O'Dea et al., 2016; Orlando et al., 2013) and references therein). The Asiatic ass species (E. kiang, E. hemionus) are not reciprocally monophyletic based on the analyzed mitochondrial genomes, and so the Asiatic ass clade is shown as ‘E. kiang + hemionus’. Daggers denote extinct taxa. NW: New World.

All analyses supported topology one in Appendix 2—figure 3. HPD: highest posterior density.

(A) Read mapping statistics. (B) Relative transversion frequencies for approaches 1–3. (C) Relative private transversion frequencies for approach 4. DNA extraction 1: (Rohland et al., 2010); DNA extraction 2: (Dabney et al., 2013b); library preparation 1: (Meyer and Kircher, 2010; Heintzman et al., 2015); library preparation 2: (Meyer and Kircher, 2010; Vilstrup et al., 2013). In (C), data in length bins with fewer than 200,000 called sites are italicized.

Minimum and maximum NWSL:Equus ratios between relative frequencies are in bold, and are used for the divergence estimates in Figure 1—figure supplement 3. Total and mean values are for the four longest bins only (90–99 to 120–129 bp). Mean values equally weight each length bin. bp: base pairs.

This topology resulted from the analysis of mtDNA data set 3 (see Appendix 1) with all partitions and Hippidion included, and dog and ceratomorphs as outgroup (not shown). Numbers above branches are Bayesian posterior probability support values from equivalent MrBayes and BEAST analyses, with those below indicating ML bootstrap values calculated in RAxML, and are shown for major nodes. (A) Full phylogeny of the analyzed equid sequences. (B) The Haringtonhippus (NWSL equid) clade, with tips color coded by geographic origin: east Beringia, blue; contiguous USA, red (following Figure 3). Tips in bold were included in the BEAST analysis (see also Supplementary file 1).

A comparison of relative private transversion frequencies between the nuclear genomes of a caballine Equus (horse, E. caballus; green), a non-caballine Equus (donkey, E. asinus; red), and 17 NWSL equids (=Haringtonhippus francisci; blue) at different read lengths, with reads divided into 10 base pair (bp) bins. Analyses are based on alignment to the horse (A) or donkey (B) genome coordinates. To account for bins with low data content, we only display comparisons with at least 200,000 observable sites.

Relative branch lengths are from Figure 1—source data 3. Minimum (darker blue) and maximum (lighter blue) estimates are shown for the NWSL equid branch.

Figure 1—figure supplement 1.

An example maximum likelihood (ML) phylogeny of equid mitochondrial genomes.

This topology resulted from the analysis of mtDNA data set 3 (see Appendix 1) with all partitions and Hippidion included, and dog and ceratomorphs as outgroup (not shown). Numbers above branches are Bayesian posterior probability support values from equivalent MrBayes and BEAST analyses, with those below indicating ML bootstrap values calculated in RAxML, and are shown for major nodes. (A) Full phylogeny of the analyzed equid sequences. (B) The Haringtonhippus (NWSL equid) clade, with tips color coded by geographic origin: east Beringia, blue; contiguous USA, red (following Figure 3). Tips in bold were included in the BEAST analysis (see also Supplementary file 1).

Appendix 2—table 3.

The a posteriori phylogenetic placement likelihood for 21 published equid mitochondrial sequences.

These analyses used the ML evolutionary placement algorithm, whilst varying the partition set used (all or reduced), and either including or excluding Hippidion sequences. Sample names are given in parentheses after the species or group name. Localities are given for NWSL equids only. Likelihoods >0.95 are in bold. *Equus includes only caballines and non-caballine equids (NCE). **For EQ04 from Alberta, other placement likelihood values for the Hippidion included/excluded partitions were: Within caballines: 0.003/0.002, Sister to caballines: 0.002/0.002, Within NCE: 0.246/0.245, Sister to NCE: 0.004/0.003. No placements were returned for ‘within Hippidion’. bp: base pairs.

Hippidion?	Partition	Published sample	Sequence length (bp)	Locality	Placement
					Sister to E. ovodovi	Sister to Hippidion	Within NWSL	Sister to NWSL	Sister to Equus*	Other**
Included	All	E. ovodovi (ACAD2305)	688		1.000
		E. ovodovi (ACAD2302)	688		1.000
		E. ovodovi (ACAD2303)	688		1.000
		H. devillei (ACAD3615)	476		N/A	1.000
		H. devillei (ACAD3625)	543		N/A	1.000
		H. devillei (ACAD3627)	543		N/A	1.000
		H. devillei (ACAD3628)	543		N/A	0.999
		H. devillei (ACAD3629)	476		N/A	0.999
		NWSL equid (JW125)	720	Klondike, YT	N/A		0.996
		NWSL equid (JW126)	720	Klondike, YT	N/A		0.999
Included	All	NWSL equid (EQ01)	620	Dry Cave, NM	N/A		0.735	0.256
		NWSL equid (EQ03)	117	Dry Cave, NM	N/A	0.002	0.974	0.011	0.003
		NWSL equid (EQ04)	117	Edmonton, AB	N/A	0.004	0.703	0.014	0.007	0.255
		NWSL equid (EQ09)	620	Natural Trap Cave, WY	N/A		0.981	0.014
		NWSL equid (EQ13)	620	Natural Trap Cave, WY	N/A		0.992
		NWSL equid (EQ16)	464	Dry Cave, NM	N/A		0.854	0.138
		NWSL equid (EQ22)	620	Natural Trap Cave, WY	N/A		0.999
		NWSL equid (EQ30)	393	San Josecito Cave, MX-NL	N/A		0.792	0.198
		NWSL equid (EQ41)	398	Natural Trap Cave, WY	N/A		0.997
		NWSL equid (JW328)	mitogenome	Mineral Hill Cave, NV	N/A		1.000
		NWSL equid (MS272)	mitogenome	Klondike, YT	N/A			1.000
	Reduced	NWSL equid (JW328)	mitogenome	Mineral Hill Cave, NV	N/A		0.996
		NWSL equid (MS272)	mitogenome	Klondike, YT	N/A			1.000
Excluded	All	E. ovodovi (ACAD2305)	688		1.000	N/A			N/A
		E. ovodovi (ACAD2302)	688		1.000	N/A			N/A
		E. ovodovi (ACAD2303)	688		1.000	N/A			N/A
		NWSL equid (JW125)	720	Klondike, YT	N/A	N/A	0.996		N/A
		NWSL equid (JW126)	720	Klondike, YT	N/A	N/A	0.999		N/A
		NWSL equid (EQ01)	620	Dry Cave, NM	N/A	N/A	0.731	0.259	N/A
		NWSL equid (EQ03)	117	Dry Cave, NM	N/A	N/A	0.980	0.010	N/A
		NWSL equid (EQ04)	117	Edmonton, AB	N/A	N/A	0.721	0.013	N/A	0.252
		NWSL equid (EQ09)	620	Natural Trap Cave, WY	N/A	N/A	0.987	0.008	N/A
		NWSL equid (EQ13)	620	Natural Trap Cave, WY	N/A	N/A	0.993		N/A
		NWSL equid (EQ16)	464	Dry Cave, NM	N/A	N/A	0.844	0.148	N/A
		NWSL equid (EQ22)	620	Natural Trap Cave, WY	N/A	N/A	0.999		N/A
		NWSL equid (EQ30)	393	San Josecito Cave, MX-NL	N/A	N/A	0.788	0.203	N/A
		NWSL equid (EQ41)	398	Natural Trap Cave, WY	N/A	N/A	0.995		N/A
		NWSL equid (JW328)	mitogenome	Mineral Hill Cave, NV	N/A	N/A	1.000		N/A
		NWSL equid (MS272)	mitogenome	Klondike, YT	N/A	N/A		1.000	N/A
	Reduced	NWSL equid (JW328)	mitogenome	Mineral Hill Cave, NV	N/A	N/A	0.995		N/A
		NWSL equid (MS272)	mitogenome	Klondike, YT	N/A	N/A		1.000	N/A

Bayesian time tree analysis results, with support and estimated divergence times for major nodes, and the tMRCAs for Haringtonhippus, E. asinus, and E. quagga summarized.

Appendix 2—figure 3.

Seven phylogenetic hypotheses for the four major groups of equids with sequenced mitochondrial genomes.

These major groups are Hippidion, the New World stilt-legged equids (=Haringtonhippus), non-caballine Equus (asses, zebras, and E. ovodovi) and caballine Equus (horses). (A) imbalanced and (B) balanced hypotheses. The hypotheses presented in (C) and (D) are identical to (A) and (B), except that Hippidion is excluded. Node letters are referenced in Appendix 2—tables 1–2. We only list combinations that were recovered by our palaeogenomic, or previous palaeogenetic, analyses.

Statistics from the phylogenetic inference analyses of nuclear genomes using all four approaches.

Summary of nuclear genome data from all 17 NWSL equids pooled together and analyzed using approach four.

Figure 1—figure supplement 3.

Calculation of divergence date estimates from nuclear genome data.

Relative branch lengths are from Figure 1—source data 3. Minimum (darker blue) and maximum (lighter blue) estimates are shown for the NWSL equid branch.

Figure 3.

The geographic distribution of Haringtonhippus.

Blue circles are east Beringian localities (KL: Klondike region, Yukon Territory, Canada). Red circles are contiguous USA localities (NTC: Natural Trap Cave, Wyoming, USA; GC: Gypsum Cave, Nevada, USA; MHC: Mineral Hill Cave, Nevada, USA; DC: Dry Cave, New Mexico, USA [Barrón-Ortiz et al., 2017; Weinstock et al., 2005]). Orange circles are localities with tentatively assigned Haringtonhippus specimens only (FB: Fairbanks, Alaska, USA; ED: Edmonton, Alberta, Canada, USA; SJC: San Josecito Cave, Nuevo Leon, Mexico; (Barrón-Ortiz et al., 2017; Guthrie, 2003). The green-star-labeled HT is the locality of the francisci holotype, Wharton County, Texas, USA. This figure was drawn using Simplemappr (Shorthouse, 2010).

To confirm the mtDNA result that NWSL equids fall outside of crown group equid diversity, we sequenced and compared partial nuclear genomes from 17 NWSL equids to a caballine (horse) and a non-caballine (donkey) reference genome. After controlling for reference genome and ancient DNA fragment length artifacts (Appendices 1–2), we examined differences in relative private transversion frequency between these genomes (Appendix 1—figure 1). We found that the relative private transversion frequency for NWSL equids was ~1.4–1.5 times greater than that for horse or donkey (Appendix 2, Figure 1—source data 3, Figure 1—figure supplement 2, and Figure 1—source data 2). This result supports the placement of NWSL equids as sister to the horse-donkey clade (Figure 1—figure supplement 3), the latter of which is representative of living Equus diversity (e.g. [Der Sarkissian et al., 2015; Jónsson et al., 2014]) and is therefore congruent with the mitochondrial genomic analyses.

Appendix 1—figure 1.

An overview of the nuclear genome analysis pipeline.

A first reference genome sequence (red; step 1) is divided into 150 bp pseudo-reads, tiled every 75 bp for exactly 2 × genomic coverage (step 2). These pseudo-reads are then mapped to a second reference genome (blue; step 3), and a consensus sequence of the mapped pseudo-reads is called (step 4). Regions of the second reference genome that are not covered by the pseudo-reads are masked (step 5). For each NWSL equid sample, reads (orange) are mapped independently to the first reference consensus sequence (step 6a) and masked second reference genome (step 6b). Alignments from steps 6a and 6b are then merged (step 7). For alignment coordinates that have base calls for the first reference, second reference, and NWSL equid sample genomes, the relative frequencies of private transversion substitutions (yellow stars) for each genome are calculated (step 8). The co-ordinates from the second reference genome (blue) are used for each analysis.

Figure 1—figure supplement 2.

A comparison of relative private transversion frequencies between the nuclear genomes of a horse, donkey, and 17 NWSL equids.

A comparison of relative private transversion frequencies between the nuclear genomes of a caballine Equus (horse, E. caballus; green), a non-caballine Equus (donkey, E. asinus; red), and 17 NWSL equids (=Haringtonhippus francisci; blue) at different read lengths, with reads divided into 10 base pair (bp) bins. Analyses are based on alignment to the horse (A) or donkey (B) genome coordinates. To account for bins with low data content, we only display comparisons with at least 200,000 observable sites.

Divergence times of Hippidion, NWSL equids, and Equus

We estimated the divergence times between the lineages leading to Hippidion, the NWSL equids, and Equus. We first applied a Bayesian time-tree approach to the whole mitochondrial genome data. This gave divergence estimates for the Hippidion-NWSL/Equus split (node 1) at 5.15–7.66 Ma, consistent with (Der Sarkissian et al., 2015), the NWSL-Equus split (node 2) at 4.09–5.13 Ma, and the caballine/non-caballine Equus split (node 3) at 3.77–4.40 Ma (Figure 1 and Figure 1—source data 1). These estimates suggest that the NWSL-Equus mitochondrial split occurred only ~500 thousand years (ka) prior to the caballine/non-caballine Equus split. We then estimated the NWSL-Equus divergence time using relative private transversion frequency ratios between the nuclear genomes, assuming a caballine/non-caballine Equus divergence estimate of 4–4.5 Ma (Orlando et al., 2013) and a genome-wide strict molecular clock (following [Heintzman et al., 2015]). This analysis yielded a divergence estimate of 4.87–5.69 Ma (Figure 1—figure supplement 3), which overlaps with that obtained from the relaxed clock analysis of whole mitochondrial genome data (Figure 1). These analyses suggest that the NWSL equid and Equus clades diverged during the latest Miocene or early Pliocene (4.1–5.7 Ma; late Hemphillian or earliest Blancan).

Systematic palaeontology

The genus Equus (Linnaeus, 1758) was named to include three living equid groups – horses (E. caballus), donkeys (E. asinus), and zebras (E. zebra) – whose diversity comprises all extant, or crown group, equids. Previous palaeontological and palaeogenetic studies have uniformly placed NWSL equids within the diversity of extant equids and therefore this genus (Barrón-Ortiz et al., 2017; Bennett, 1980; Der Sarkissian et al., 2015; Harington and Clulow, 1973; Orlando et al., 2008; 2009; Scott, 2004; Vilstrup et al., 2013; Weinstock et al., 2005). This, however, conflicts with the phylogenetic signal provided by palaeogenomic data, which strongly suggest that NWSL equids fall outside the confines of the equid crown group (Equus). Nor is there any morphological or genetic evidence warranting the assignment of NWSL equids to an existing extinct taxon such as Hippidion. We therefore erect a new genus for NWSL equids, Haringtonhippus, as defined and delimited below:

Order: Perissodactyla, Owen 1848

Family: Equidae, Linnaeus 1758

Subfamily: Equinae, Steinmann & Döderlein 1890

Tribe: Equini, Gray 1821

Genus: Haringtonhippus, gen. nov. urn:lsid:zoobank.org:act:35D901A7-65F8-4615-9E13-52A263412F67

Type species. Haringtonhippus francisci Hay 1915.

Etymology

The new genus is named in honor of C. Richard Harington, who first described NWSL equids from eastern Beringia (Harington and Clulow, 1973). ‘Hippus’ is from the Greek word for horse, and so Haringtonhippus is implied to mean ‘Harington’s horse’.

Holotype

A partial skeleton consisting of a complete cranium, mandible, and a stilt-legged third metatarsal (MTIII) (Figure 2a and Figure 2—figure supplement 1b), which is curated at the Texas Vertebrate Paleontology Collections at The University of Texas, Austin (TMM 34–2518). This specimen is the holotype of ‘E’. francisci, originally described by Hay (1915), and is from the middle Pleistocene Lissie Formation of Wharton County, Texas (Hay, 1915; Lundelius and Stevens, 1970).

Figure 2.

Morphological analysis of extant and middle-late Pleistocene equids.

(A) Crania of Haringtonhippus francisci, upper: LACM(CIT) 109/156450 from Nevada, lower: TMM 34–2518 from Texas. (B) From upper to lower, third metatarsals of: H. francisci (YG 401.268), E. lambei (YG 421.84), and E. cf. scotti (YG 198.1) from Yukon. Scale bar is 5 cm. (C) Principal component analysis of selected third metatarsals from extant and middle-late Pleistocene equids, showing clear clustering of stilt-legged (hemionine Equus (orange) and H. francisci (green)) from stout-legged (caballine Equus; blue) specimens (see also Figure 2—source data 1). Symbol shape denotes the specimen identification method (DNA: square, triangle: DNA/morphology, circle: morphology). The first and second principal components explain 95% of the variance.

(A) LACM(CIT) 109/156450 from Nevada, identified through mitochondrial and nuclear palaeogenomic analysis. Upper: right side (reflected for comparison), lower: left side. (B) Part of the H. francisci holotype, TMM 34–2518 from Texas.

(A) Third metatarsals from H. francisci (upper; YG 401.268), E. lambei (middle; YG 421.84), and E. cf. scotti (lower; YG 198.1). (B) Third metacarpals from H. francisci (upper; YG 404.663), E. lambei (middle; YG 109.6), and E. cf. scotti (lower; YG 378.15). (C) Proximal fragments of radii from H. francisci (left; YG 303.1085), and E. lambei (right; YG 303.325). (D) First phalanges from H. francisci (left; YG 130.3), E. lambei (middle; YG 404.22), and E. cf. scotti (right; YG 168.1).

This specimen (KU 47800; JK260) was originally referred to Equus sp., but is here identified as H. francisci on the basis of mitochondrial and nuclear genome data. We note the relative slenderness of this specimen, which is comparable to YG 404.663 (H. francisci) from Yukon in Figure 2—figure supplement 2.

This specimen (LACM(CIT)109/150708; JW277/JK166) was originally identified by Weinstock et al., 2005. List of Source data files.

Figure 2—figure supplement 1.

The two crania assigned to H. francisci.

(A) LACM(CIT) 109/156450 from Nevada, identified through mitochondrial and nuclear palaeogenomic analysis. Upper: right side (reflected for comparison), lower: left side. (B) Part of the H. francisci holotype, TMM 34–2518 from Texas.

An example equid metacarpal from Natural Trap Cave, Wyoming.

Figure 2—figure supplement 2.

Comparison between the limb bones of H. francisci, E. lambei, and E. cf. scotti from Yukon.

(A) Third metatarsals from H. francisci (upper; YG 401.268), E. lambei (middle; YG 421.84), and E. cf. scotti (lower; YG 198.1). (B) Third metacarpals from H. francisci (upper; YG 404.663), E. lambei (middle; YG 109.6), and E. cf. scotti (lower; YG 378.15). (C) Proximal fragments of radii from H. francisci (left; YG 303.1085), and E. lambei (right; YG 303.325). (D) First phalanges from H. francisci (left; YG 130.3), E. lambei (middle; YG 404.22), and E. cf. scotti (right; YG 168.1).

Referred material

On the basis of mitochondrial and nuclear genomic data, we assign the following material confidently to Haringtonhippus: a cranium, femur, and MTIII (LACM(CIT): Nevada); three MTIIIs, three third metacarpals (MCIII), three premolar teeth, and a molar tooth (KU: Wyoming); two radii, 12 MTIIIs, three MCIIIs, a metapodial, and a first phalanx (YG: Yukon Territory); and a premolar tooth (University of Texas El Paso, UTEP: New Mexico); (Figure 2—figure supplements 1–4 and Supplementary file 1; (Barrón-Ortiz et al., 2017; Weinstock et al., 2005). This material includes at least four males and at least six females (Appendix 2, Appendix 2—Table 4 and Appendix 2—Table 4—source data 1). We further assign MTIII specimens from Yukon Territory (n = 13), Wyoming (n = 57), and Nevada (n = 4) to Haringtonhippus on the basis of morphometric analysis (Figure 2c and Figure 2—source data 1). On the basis of short mitochondrial DNA sequences, we tentatively assign to Haringtonhippus a premolar tooth (LACM(CIT): Nuevo Leon); a premolar and a molar (UTEP: New Mexico); and a premolar (Royal Alberta Museum, RAM/PMA: Alberta) (Barrón-Ortiz et al., 2017). We also tentatively assign 19 NWSL equid metapodial specimens from the Fairbanks area, Alaska (Guthrie, 2003) to Haringtonhippus, but note that morphometric and/or palaeogenomic analysis would be required to confirm this designation.

Figure 2—figure supplement 4.

An example femur of H. francisci from Gypsum Cave, Nevada.

This specimen (LACM(CIT)109/150708; JW277/JK166) was originally identified by Weinstock et al., 2005. List of Source data files.

Appendix 2—table 4.

Sex determination analysis of 17 NWSL equids.

Chromosome ratio is the relative mapping frequency ratio between all autosomes and the X-chromosome. Males are inferred if the ratio is 0.45–0.55 and females if the ratio is 0.9–1.1.

Sample	Museum accession	Chromosome ratio	Inferred sex
AF037	YG 402.235	0.48	male
JK166	LACM(CIT) 109/150807	0.93	female
JK167	LACM(CIT) 109/149291	0.91	female
JK207	LACM(CIT) 109/156450	0.92	female
JK260	KU 47800	0.95	female
JK276	KU 53678	0.91	female
MS341	YG 303.1085	0.50	male
MS349	YG 130.55	0.48	male
MS439	YG 401.387	0.98	female
PH008	YG 404.205	0.90	female
PH013	YG 130.6	0.87	probable female
PH014	YG 303.371	0.46	male
PH015	YG 404.662	0.44	probable male
PH021	YG 29.169	0.83	probable female
PH023	YG 160.8	0.91	female
PH036	YG 76.2	0.81	probable female
PH047	YG 404.663	0.88	probable female

Geographic and temporal distribution

Haringtonhippus is known only from the Pleistocene of North America (Figure 3). In addition to the middle Pleistocene holotype from Texas, Haringtonhippus is confidently known from the late Pleistocene of Yukon Territory (Klondike region), Wyoming (Natural Trap Cave), Nevada (Gypsum Cave, Mineral Hill Cave), and New Mexico (Dry Cave), and is tentatively registered as present in Nuevo Leon (San Josecito Cave), Alberta (Edmonton area), and Alaska (Fairbanks area) (Appendix 2, Supplementary file 1, and Appendix 2—table 3; [Barrón-Ortiz et al., 2017; Vilstrup et al., 2013; Weinstock et al., 2005]).

To investigate the last appearance date (LAD) of Haringtonhippus in eastern Beringia, we obtained new radiocarbon dates from 17 Yukon Territory fossils (Appendix 1 and Supplementary file 1). This resulted in three statistically-indistinguishable radiocarbon dates of ~14.4 14C ka BP (derived from two independent laboratories) from a metacarpal bone (YG 401.235) of Haringtonhippus, which represents this taxon’s LAD in eastern Beringia (Supplementary file 1). The LAD for North America as a whole is based on two dates of ~13.1 14C ka BP from Gypsum Cave, Nevada (Supplementary file 1; [Weinstock et al., 2005]).

Mitogenomic diagnosis

Haringtonhippus is the sister genus to Equus (equid crown group), with Hippidion being sister to the Haringtonhippus-Equus clade (Figure 1). Haringtonhippus can be differentiated from Equus and Hippidion by 178 synapomorphic positions in the mitochondrial genome, including four insertions and 174 substitutions (Appendix 1—Table 2 and Appendix 1—table 2—source data 1). We caution that these synapomorphies are tentative and will likely be reduced in number as a greater diversity of mitochondrial genomes for extinct equids become available.

Appendix 1—table 2.

Summary of the number and type of synapomorphic bases for each of the three examined equid genera.

A full list of these substitutions, and their position relative to the E. caballus reference mitochondrial genome (NC_001640), can be found in Appendix 1—table 2-Source data 1. *total includes a further five synapomorphic sites that have unique states in each genus.

The horse reference mtDNA has Genbank accession NC_001640.1.

Substitution	Hippidion	Haringtonhippus	Equus
Transition	338	147	66
Transversion	43	22	4
Insertion	2	4	0
Deletion	3	0	0
Total*	391	178	75

Morphological comparisons of third metatarsals

We used morphometric analysis of caballine/stout-legged Equus and stilt-legged equids (hemionine/stilt-legged Equus, Haringtonhippus) MTIIIs to determine how confidently these groups can be distinguished (Figure 2c). Using logistic regression on principal components, we find a strong separation that can be correctly distinguished with 98.2% accuracy (Appendix 2; Heintzman et al., 2017). Hemionine/stilt-legged Equus MTIIIs occupy the same morphospace as H. francisci in our analysis, although given a larger sample size, it may be possible to discriminate E. hemionus from the remaining stilt-legged equids. We note that Haringtonhippus seems to exhibit a negative correlation between latitude and MTIII length, and that specimens from the same latitude occupy similar morphospace regardless of whether DNA- or morphological-based identification was used (Figure 2c and Figure 2—source data 1).

Comments

On the basis of morphology, we assign all confidently referred material of Haringtonhippus to the single species H. francisci Hay (1915) (Appendix 2). Comparison between the cranial anatomical features of LACM(CIT) 109/156450 and TMM 34–2518 reveal some minor differences, which can likely be ascribed to intraspecific variation (Figure 2a and Appendix 2 and Figure 2—figure supplement 1). Further, the MTIII of TMM 34–2518 is comparable to the MTIIIs ascribed to Haringtonhippus by palaeogenomic data, and is consistent with the observed latitudinally correlated variation in MTIII length across Haringtonhippus (Figure 2c and Appendix 2).

This action is supported indirectly by molecular evidence, namely the lack of mitochondrial phylogeographic structure and the estimated time to most recent common ancestor (tMRCA) for sampled Haringtonhippus. The mitochondrial tree topology within Haringtonhippus does not exhibit phylogeographic structure (Figure 1—figure supplement 1b), which is consistent with sampled Haringtonhippus mitochondrial genomes belonging to the same species. Using Bayesian time-tree analysis, we estimated a tMRCA for the sampled Haringtonhippus mitochondrial genomes of ~200–470 ka BP (Figure 1 and Figure 1—source data 1; Heintzman et al., 2017). The MRCA of Haringtonhippus is therefore more recent than that of other extant equid species (such as E. asinus and E. quagga, which have a combined 95% HPD range: 410–1030 ka BP; Figure 1 and Figure 1—source data 1; Heintzman et al., 2017). Although the middle Pleistocene holotype TMM 34–2518 (~125–780 ka BP) may predate our Haringtonhippus mitochondrial tMRCA, this sample has no direct date and the range of possible ages falls within the tMRCA range of other extant equid species. We therefore cannot reject the hypothesis of its conspecificity with Haringtonhippus, as defined palaeogenomically. We attempted, but were unable, to recover either collagen or genomic data from TMM 34–2518 (Appendix 2), consistent with the taphonomic, stratigraphic, and geographic context of this fossil (Hay, 1915; Lundelius and Stevens, 1970). Altogether, the molecular evidence is consistent with the assignment of H. francisci as the type and only species of Haringtonhippus.

Discussion

Reconciling the genomic and fossil records of Plio-Pleistocene equid evolution

The suggested placement of NWSL equids within a taxon (Haringtonhippus) sister to Equus is a departure from previous interpretations, which variably place the former within Equus, as sister to hemionines or caballine horses (Figure 1). According to broadly accepted palaeontological interpretations, the earliest equids exhibiting morphologies consistent with NWSL and caballine attribution appear in the fossil record only ~2–3 and ~1.9–0.7 Ma ago (Eisenmann et al., 2008; Forsten, 1992), respectively, whereas our divergence estimates suggest that these lineages to have diverged between 4.1–5.8 and 3.8–4.5 Ma, most likely in North America. Dating incongruence might be attributed to an incomplete fossil record, but this seems unlikely given the density of the record for late Neogene and Pleistocene horses. Conversely, incongruence might be attributed to problems with estimating divergence using genomic evidence. However, we emphasize that the NWSL-Equus split is robustly calibrated to the caballine/non-caballine Equus divergence at 4.0–4.5 Ma, which is in turn derived from a direct molecular clock calibration using a middle Pleistocene horse genome (Orlando et al., 2013).

Other possibilities to explain the incongruence include discordance between the timing of species divergence and the evolution of diagnostic anatomical characteristics, or failure to detect or account for homoplasy (Forsten, 1992). For example, Pliocene Equus generally exhibits a primitive (‘plesippine’ in North America, ‘stenonid’ in the Old World) morphology that presages living zebras and asses (Forsten, 1988, 1992), with more derived caballine (stout-legged) and hemionine (stilt-legged) forms evolving in the early Pleistocene. The stilt-legged morphology appears to have evolved independently at least once in each of the Old and New Worlds, yielding the Asiatic wild asses and Haringtonhippus, respectively. We include the middle-late Pleistocene Eurasian E. hydruntinus within the Asiatic wild asses (following [Bennett et al., 2017; Burke et al., 2003; Orlando et al., 2006]), and note that the Old World sussemione E. ovodovi may represent another instance of independent stilt-legged origin, but its relation to Asiatic wild asses and other non-caballine Equus is currently unresolved (as depicted in Der Sarkissian et al., 2015; Orlando et al., 2009; Vilstrup et al., 2013; and Figure 1). It is plausible that features at the plesiomorphous end of the spectrum, such as those associated with Hippidion, survived after the early to middle Pleistocene at lower latitudes (South America, Africa; Figure 1). By contrast, the more derived hemionine and caballine morphologies evolved from, and replaced, their antecedents in higher latitude North America and Eurasia, perhaps as adaptations to the extreme ecological pressures perpetuated by the advance and retreat of continental ice sheets and correlated climate oscillations during the Pleistocene (Forsten, 1992, Forsten, 1996Forsten, 1996). We note that this high-latitude replacement model is consistent with the turnover observed in regional fossil records for Pleistocene equids in North America (Azzaroli, 1992; Azzaroli and Voorhies, 1993) and Eurasia (Forsten, 1988, 1992, Forsten, 1996). By contrast, in South America Hippidion co-existed with caballine horses until they both succumbed to extinction, together with much of the New World megafauna near the end of the Pleistocene (Forsten, 1996; Koch and Barnosky, 2006; O'Dea et al., 2016). This model helps to explain the discordance between the timings of the appearance of the caballine and hemionine morphologies in the fossil record and the divergence of lineages leading to these forms as estimated from palaeogenomic data.

Although we can offer no solution to the general problem of mismatches between molecular and morphological divergence estimators–an issue scarcely unique to equid systematics–this model predicts that some previously described North American Pliocene and early Pleistocene Equus species (e.g. E. simplicidens, E. idahoensis; [Azzaroli and Voorhies, 1993]), or specimens thereof, may be ancestral to extant Equus and/or late Pleistocene Haringtonhippus.

Temporal and geographic range overlap of Pleistocene equids in North America

Three new radiocarbon dates of ~14.4 14C ka BP from a Yukon Haringtonhippus fossil greatly extends the known temporal range of this genus in eastern Beringia. This result demonstrates, contrary to its previous LAD of 31,400 ± 1200 14C years ago (AA 26780; [Guthrie, 2003]), that Haringtonhippus survived throughout the last glacial maximum in eastern Beringia (Clark et al., 2009) and may have come into contact with humans near the end of the Pleistocene (Goebel et al., 2008; Guthrie, 2006). These data suggest that populations of stilt-legged Haringtonhippus and stout-legged caballine Equus were sympatric, both north and south of the continental ice sheets, through the late Pleistocene and became extinct at roughly the same time. The near synchronous extinction of both horse groups across their entire range in North America suggests that similar causal mechanisms may have led each to their demise.

The sympatric nature of these equids raises questions of whether they managed to live within the same community without hybridizing or competing for resources. Extant members of the genus Equus vary considerably in the sequence of Prdm9, a gene involved in the speciation process, and chromosome number (karyotype) (Ryder et al., 1978; Steiner and Ryder, 2013), and extant caballine and non-caballine Equus rarely produce fertile offspring (Allen and Short, 1997; Steiner and Ryder, 2013). It is unlikely, therefore, that the more deeply diverged Haringtonhippus and caballine Equus would have been able to hybridize. Future analysis of high coverage nuclear genomes, ideally including an outgroup such as Hippidion, will make it possible to test for admixture that may have occurred soon after the lineages leading to Haringtonhippus and Equus diverged, as occurred between the early caballine and non-caballine Equus lineages (Jónsson et al., 2014). It may also be possible to use isotopic and/or tooth mesowear analyses to assess the potential of resource partitioning between Haringtonhippus and caballine Equus in the New World.

Fossil systematics in the palaeogenomics and proteomics era: concluding remarks

Fossils of NWSL equids have been known for more than a century, but until the present study their systematic position within Plio-Pleistocene Equidae was poorly characterized. This was not because of a lack of interest on the part of earlier workers, whose detailed anatomical studies strongly indicated that what we now call Haringtonhippus was related to Asiatic wild asses, such as Tibetan khulan and Persian onagers, rather than to caballine horses (Eisenmann et al., 2008; Guthrie, 2003; Scott, 2004; Skinner and Hibbard, 1972). That the cues of morphology have turned out to be misleading in this case underlines a recurrent problem in systematic biology, which is how best to discriminate authentic relationships within groups, such as Neogene equids, that were prone to rampant convergence. The solution we adopted here was to utilize both palaeogenomic and morphometric information in reframing the position of Haringtonhippus, which now clearly emerges as the closest known outgroup to all living Equus.

Our success in this regard demonstrates that an approach which incorporates phenomics with molecular methods (palaeogenomic as well as palaeoproteomic, e.g. [Welker et al., 2015]) is likely to offer a means for securely detecting relationships within speciose groups that are highly diverse ecomorphologically. All methods have their limits, with taphonomic degradation being the critical one for molecular approaches. However, proteins may persist significantly longer than ancient DNA (e.g. [Rybczynski et al., 2013]), and collagen proteomics may come to play a key role in characterizing affinities, as well as the reality, of several proposed Neogene equine taxa (e.g. Dinohippus, Pliohippus, Protohippus, Calippus, and Astrohippus; [MacFadden, 1998]) whose distinctiveness and relationships are far from settled (Azzaroli and Voorhies, 1993; Forsten, 1992). A reciprocally informative approach like the one taken here holds much promise for lessening the amount of systematic noise, due to oversplitting, that hampers our understanding of the evolutionary biology of other major late Pleistocene megafaunal groups such as bison and mammoths (Enk et al., 2016; Froese et al., 2017). This approach is clearly capable of providing new insights into just how extensive megafaunal losses were at the end of the Pleistocene, in what might be justifiably called the opening act of the Sixth Mass Extinction in North America.

Materials and methods

We provide an overview of methods here; full details can be found in Appendix 1.

Sample collection and radiocarbon dating

We recovered Yukon fossil material (17 Haringtonhippus francisci, two Equus cf. scotti, and two E. lambei; Supplementary file 1) from active placer mines in the Klondike goldfields near Dawson City. We further sampled seven H. francisci fossils from the contiguous USA that are housed in collections at the University of Kansas Biodiversity Institute (KU; n = 4), Los Angeles County Museum of Natural History (LACM(CIT); n = 2), and the Texas Vertebrate Paleontology Collections at The University of Texas (TMM; n = 1). We radiocarbon dated the Klondike fossils and the H. francisci cranium from the LACM(CIT) (Supplementary file 1).

Morphometric analysis of third metatarsals

For morphometric analysis, we took measurements of third metatarsals (MTIII) and other elements. We used a reduced data set of four MTIII variables for principal components analysis and performed logistic regression on the first three principal components, computed in R (R Development Core Team, 2008) (Source code 1).

DNA extraction, library preparation, target enrichment, and sequencing

We conducted all molecular biology methods prior to indexing PCR in the dedicated palaeogenomics laboratory facilities at either the UC Santa Cruz or Pennsylvania State University. We extracted DNA from between 100 and 250 mg of bone powder following either Rohland et al. (2010) or Dabney et al. (2013a). We then converted DNA extracts to libraries following the Meyer and Kircher protocol (Meyer and Kircher, 2010), as modified by (Heintzman et al., 2015) or the PSU method of (Vilstrup et al., 2013). We enriched libraries for equid mitochondrial DNA. We then sequenced all enriched libraries and unenriched libraries from 17 samples using Illumina platforms.

Mitochondrial genome reconstruction and analysis

We prepared raw sequence data for alignment and mapped the filtered reads to the horse reference mitochondrial genome (Genbank: NC_001640.1) and a H. francisci reference mtDNA genome (Genbank: KT168321), resulting in mitogenomic coverage ranging from 5.8× to 110.7× (Supplementary file 1). We were unable to recover equid mtDNA from TMM 34–2518 (the francisci holotype) using this approach (Appendix 2). We supplemented our mtDNA genome sequences with 38 previously published complete equid mtDNA genomes. We constructed six alignment data sets and selected models of molecular evolution for the analyses described below (Appendix 1—table 1, and Supplementary file 1; Heintzman et al., 2017).

Appendix 1—table 1.

Selected models of molecular evolution for partitions of the first five mtDNA genome alignment data sets.

All lengths are in base pairs. Reduced length excludes the Coding3 and CR partitions. For all RAxML analyses the GTR model was implemented. *The TrN model was selected, but this cannot be implemented in MrBayes and so the HKY model was used. EPA: evolutionary placement algorithm; CR: control region.

Data set		Partition						Total length
		Coding1	Coding2	Coding3	rRNAs	tRNAs	CR	All	Reduced
1. White rhino outgroup	Length	3803	3803	3803	2579	1529	1066	16583	11714
	Model	GTR + I + G	HKY + I + G	GTR + I + G	GTR + I + G	HKY + I + G	HKY*+I + G
2. Malayan tapir outgroup	Length	3803	3803	3803	2585	1530	1065	16589	11721
	Model	GTR + I + G	HKY + I + G	GTR + I + G	GTR + I + G	HKY + I + G	HKY*+G
3. Dog + ceratomorphs outgroups	Length	3803	3803	3803	2615	1540	N/A	15564	11761
	Model	GTR + I + G	HKY + I + G	GTR + I + G	GTR + I + G	HKY + I + G	N/A
4. EPA	Length	3803	3803	3803	2601	1534	1118	16662	11741
	Model	GTR + I + G	TrN + I + G	GTR + I + G	GTR + I + G	HKY + I + G	HKY + I + G
5. Equids only	Length	3802	3802	3802	2571	1528	971	16476	11703
	Model	TrN + I + G	TrN + I + G	GTR + G	TrN + I + G	HKY + I	HKY + G

We tested the phylogenetic position of the NWSL equids (=H. francisci) using mtDNA data sets 1–3 and applying Bayesian (Ronquist et al., 2012) and maximum likelihood (ML; [Stamatakis, 2014]) analyses. We varied the outgroup, the inclusion or exclusion of the fast-evolving partitions, and the inclusion or exclusion of Hippidion sequences. Due to the lack of a globally supported topology across the Bayesian and ML phylogenetic analyses, we used an Evolutionary Placement Algorithm (EPA; [Berger et al., 2011]) to determine the a posteriori likelihood of phylogenetic placements for candidate equid outgroups using mtDNA data set four. We also used the same approach to assess the placement of previously published equid sequences (Appendix 2). To infer divergence times between the four major equid groups (Hippidion, NWSL equids, caballine Equus, and non-caballine Equus), we ran Bayesian timetree analyses (Drummond et al., 2012) using mtDNA data set five. We varied these analyses by including or excluding fast-evolving partitions, constrained the root height or not, and including or excluding the E. ovodovi sequence.

To facilitate future identification of equid mtDNA sequences, we constructed, using data set six, a list of putative synapomorphic base states, including indels and substitutions, that define the genera Hippidion, Haringtonhippus, and Equus at sites across the mtDNA genome.

Phylogenetic inference, divergence date estimation, and sex determination from nuclear genomes

To test whether our mtDNA genome-based phylogenetic hypothesis truly reflects the species tree, we compared the nuclear genomes of a horse (EquCab2), donkey (Orlando et al., 2013), and the shotgun sequence data from 17 of our NWSL equid samples (Figure 1—source data 2, Appendix 1, Appendix 1—figure 1, and Supplementary file 1). We applied four successive approaches, which controlled for reference genome and DNA fragment length biases (Appendix 1).

We estimated the divergence between the NWSL equids and Equus (horse and donkey) by fitting the branch length, or relative private transversion frequency, ratio between horse/donkey and NWSL equids into a simple phylogenetic scenario (Figure 1—figure supplement 3). We then multiplied the NWSL equid branch length by a previous horse-donkey divergence estimate (4.0–4.5 Ma; [Orlando et al., 2013]) to give the estimated NWSL equid-Equus divergence date, following (Heintzman et al., 2015) and assuming a strict genome-wide molecular clock (Figure 1—figure supplement 3).

We determined the sex of the 17 NWSL equid samples by comparing the relative mapping frequency of the autosomes to the X chromosome.

DNA damage analysis

We assessed the prevalence of mitochondrial and nuclear DNA damage in a subset of the equid samples using mapDamage (Jónsson et al., 2013).

Data availability

Repository details and associated metadata for curated samples can be found in Supplementary file 1. MTIII and other element measurement data are in Figure 2—source data 1, and the Rscript used for morphometric analysis is in the DRYAD database (Heintzman et al., 2017). MtDNA genome sequences have been deposited in Genbank under accessions KT168317-KT168336, MF134655-MF134663, and an updated version of JX312727. All mtDNA genome alignments (in NEXUS format) and associated XML and TREE files are in the DRYAD database (Heintzman et al., 2017). Raw shotgun sequence data used for the nuclear genomic analyses and raw shotgun and target enrichment sequence data for TMM 34–2518 (francisci holotype) have been deposited in the Short Read Archive (BioProject: PRJNA384940).

Nomenclatural act

The electronic edition of this article conforms to the requirements of the amended International Code of Zoological Nomenclature, and hence the new name contained herein is available under that Code from the electronic edition of this article. This published work and the nomenclatural act it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix ‘http://zoobank.org/'. The LSID for this publication is: urn:lsid:zoobank.org:pub:8D270E0A-9148-4089-920C-724F07D8DC0B. The electronic edition of this work was published in a journal with an ISSN, and has been archived and is available from the following digital repositories: PubMed Central and LOCKSS.

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "A new genus of horse from Pleistocene North America" for consideration by eLife. Your article has been reviewed by four peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Diethard Tautz as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Jose Luis Prado (Reviewer #1); Anna Linderholm (Reviewer #2); Maria Teresa Alberdi (Reviewer #4).

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

Summary:

This study uses paleogenomic data (from both the mitochondrial and nuclear genome) to convincingly revise the taxonomic diversity and position of the extinct New World stilt-legged (NWSL) equids. Based on their results, the authors propose a new genus, Haringtonhippus, for what they demonstrate is a single NWSL species, H. francisci.

Reviewers praised the thoroughness of the ancient DNA and phylogenomic analysis; the taxonomic revision itself is convincing. There are two essential revisions that the authors must address before the paper can be considered further for publication in eLife.

Essential revisions:

1) There were questions about whether the authors are justified to name a new genus based on the results as presented; demonstration that the NWSL are a sister taxon to non-NWSL Equus spp. is insufficient on its own. Of course there is not a straightforward cutoff available, but can the authors provide a comparative sense of scale based on the range of phylogenetic divergence (e.g. nucleotide divergence) for other genus-genus pairs in the same order? Showing comparability to conventional genus pairs in the order would help to justify the decision by the authors to name a new genus; otherwise the conservative choice may be to not do so.

2) Somewhat related to the above point, the estimated divergence time for the NWSL is much older than the actual fossil record seems to suggest. This could easily be a result of poor fossil records or it could be that the genetic analysis is not foolproof (several other studies has shown that using molecular clocks as part of divergence estimates can be problematic). The authors have some treatment of this discordance in the Discussion, but they should consider more fully the potential sources of inaccuracy in the genomic-based divergence estimates.

Essential revisions:

1) There were questions about whether the authors are justified to name a new genus based on the results as presented; demonstration that the NWSL are a sister taxon to non-NWSL Equus spp. is insufficient on its own. Of course there is not a straightforward cutoff available, but can the authors provide a comparative sense of scale based on the range of phylogenetic divergence (e.g. nucleotide divergence) for other genus-genus pairs in the same order? Showing comparability to conventional genus pairs in the order would help to justify the decision by the authors to name a new genus; otherwise the conservative choice may be to not do so.

We appreciate that the rules governing taxonomy are somewhat fluid, and also that our decision to create a new genus may be seen as not particularly conservative. We argue, however, that our decision is robust given recognized taxonomy. Contrary to all previous phylogenetic work on the NWSL equids, which suggests that they fall within genus Equus, our work establishes that it does not. To be so included, the last common ancestor of all recognized species of Equus would have to also be the ancestor of NWSL equids. Our palaeontological and phylogenetic analyses strongly indicate that any such ancestor would likely also have been ancestral to many other taxa not normally included in Equus, resulting in undesirable paraphyly. According to nomenclatural rules, because the species-level taxon comprising NWSL equids falls outside Equus as narrowly defined, it must be included in another, cognate clade at the genus level. In the absence of a pre-existing, valid name for this newly-recognized clade, we have fashioned our own, Haringtonhippus. To clarify this in the manuscript, we have refined the first paragraph of the systematic palaeontology section.

In addition, and to address this suggestion directly, we investigated during the course of this work the use of phylogenetic/nucleotide divergence to set a genus-genus cutoff. Unfortunately, the distributions of within-genus (species-species) and within-family (genus-genus) nucleotide divergences of mammals are strongly overlapping, including the within-Equus nucleotide divergences (Johns and Avise 1998). Within the order Perissodactyla and based on our full/reduced mitochondrial genome alignments (dataset four), species in the tapir genus (Tapirus) have a greater mean nucleotide divergence (10.6%/4.9%) than between-genus rhinoceros pairs (e.g. Ceratotherium-Diceros: 7.1%/3.3%; and Dicerorhinus-Coelodonta: 8.9%/3.5%). This suggests that mean nucleotide divergence is not consistent with taxonomic rank in our study group, and so would not be a conservative approach to discriminate genera here.

Johns, G. C. & Avise, J. C. (1998) A comparative summary of genetic distances in the vertebrates from the mitochondrial cytochrome b gene. Molecular Biology and Evolution, 15(11), 1481961490.

2) Somewhat related to the above point, the estimated divergence time for the NWSL is much older than the actual fossil record seems to suggest. This could easily be a result of poor fossil records or it could be that the genetic analysis is not foolproof (several other studies has shown that using molecular clocks as part of divergence estimates can be problematic). The authors have some treatment of this discordance in the Discussion, but they should consider more fully the potential sources of inaccuracy in the genomic-based divergence estimates.

This observation is correct - the fossil record suggests that NWSL equids appear ~2-3 Ma and caballine horses ~1.9-0.7 Ma, and our molecular clock-based estimates suggest that these lineages are at least twice as old as this. Differences between fossil and genetic divergence estimates are not uncommon, and this can be due to a variety of issues including an incomplete fossil record or discordance between the appearance of morphologically distinctive traits and genetically distinct lineages, poor or insufficient calibration of the molecular clock or variations in the rate of molecular evolution over time, errors in phylogenetic estimates, and a combination of these. In this specific case, it seems unlikely that the incompleteness of the fossil record could explain this discrepancy, as the equid fossil record is well described in North America. However, our molecular clock-based estimates are based on a previously estimated divergence between caballine horses and donkeys of 4.0-4.5 Ma, which is a range that was estimated by Orlando et al., 2013 in their analyses of living and a very old (560-780,000-year-old) horse genome. In fact, the observation that these divergence estimates were approximately twice as old as generally accepted was first explored in this 2013 manuscript, where the authors also noted that their result was in agreement both with previous molecular data (Vilstrup et al., 2013) and with the age of a fossil horse from Mexico, Dinohippus mexicanus, which is considered a direct cladogenetic ancestor of early Equus (MacFadden and Carranza-Castaneda 2002).

To explore this incongruence more thoroughly in our manuscript, we have revised the text in the Discussion section 'Reconciling the genomic and fossil records of Plio-Pleistocene equid evolution'. We note that this type of problem is not unique to equids, and certainly one into which genomic data from increasingly diverse lineages may help to improve, in particular as methods to recover genomic data from increasingly old and poorly preserved remains improve.

Orlando, L. et al. Recalibrating Equus evolution using the genome sequence of an early Middle Pleistocene horse. Nature, 499(7456), 749678 (2013).

Vilstrup, J. T. et al. Mitochondrial phylogenomics of modern and ancient equids. PLoS ONE 8, e55950 (2013).

MacFadden, B. J. & Carranza-Castaneda, O. Cranium of Dinohippus mexicanus (Mammalia Equidae) from the early Pliocene (latest Hemphillian) of central Mexico and the origin of Equus. Bulletin of the Florida Museum Natural History 43, 16396185 (2002).