# Detecting taxonomic and phylogenetic signals in equid cheek teeth: towards new palaeontological and archaeological proxies

T. Cucchi, A. Mohaseb, S. Peigné, K. Debue, L. Orlando, M. Mashkour

► To cite this version:

T. Cucchi, A. Mohaseb, S. Peigné, K. Debue, L. Orlando, et al.. Detecting taxonomic and phylogenetic signals in equid cheek teeth: towards new palaeontological and archaeological proxies. Royal Society Open Science, 2017, 4 (4), pp.160997. 10.1098/rsos.160997. hal-01526821

HAL Id: hal-01526821

https://hal.sorbonne-universite.fr/hal-01526821

Submitted on 23 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ROYAL SOCIETY
OPEN SCIENCE

rsos.royalsocietypublishing.org

Research

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Check for
updates

**Cite this article:** Cucchi T, Mohaseb A, Peigné S, Debue K, Orlando L, Mashkour M. 2017

Detecting taxonomic and phylogenetic signals
in equid cheek teeth: towards new
palaeontological and archaeological proxies.

*R. Soc. open sci.* **4**: 160997.

http://dx.doi.org/10.1098/rsos.160997

Received: 5 December 2016

Accepted: 7 March 2017

**Subject Category:**

Biology (whole organism)

**Subject Areas:**

palaeontology/biometrics/evolution

**Keywords:**

*Equus*, phylogenetic signal, shape, tooth
geometric morphometrics, equid evolutionary
history, fossil record

**Author for correspondence:**

T. Cucchi

e-mail: cucchi@mnhn.fr

†These authors contributed equally to this study.

Electronic supplementary material is available
online at https://dx.doi.org/10.6084/m9.
figshare.c.3722728.

Detecting taxonomic and
phylogenetic signals in
equid cheek teeth: towards
new palaeontological and
archaeological proxies

T. Cucchi<sup>1,2,†</sup>, A. Mohaseb<sup>1,†</sup>, S. Peigne<sup>3</sup>, K. Debue<sup>1</sup>,
L. Orlando<sup>4,5</sup> and M. Mashkour<sup>1</sup>

<figure><img src="image_4.png" /></figure>

TC, 0000-0001-6021-5001

The Plio-Pleistocene evolution of *Equus* and the subsequent
domestication of horses and donkeys remains poorly
understood, due to the lack of phenotypic markers
capable of tracing this evolutionary process in the
palaeontological/archaeological record. Using images from
345 specimens, encompassing 15 extant taxa of equids,
we quantified the occlusal enamel folding pattern in four
mandibular cheek teeth with a single geometric morphometric
protocol. We initially investigated the protocol accuracy by
assigning each tooth to its correct anatomical position and
taxonomic group. We then contrasted the phylogenetic signal
present in each tooth shape with an exome-wide phylogeny
from 10 extant equine species. We estimated the strength of
the phylogenetic signal using a Brownian motion model of
evolution with multivariate *K* statistic, and mapped the dental
shape along the molecular phylogeny using an approach based
on squared-change parsimony. We found clear evidence for
the relevance of dental phenotypes to accurately discriminate
all modern members of the genus *Equus* and capture their
phylogenetic relationships. These results are valuable for both
palaeontologists and zooarchaeologists exploring the spatial

<footnote n="1"><sup>1</sup>CNRS, Muséum national d'Histoire naturelle, Sorbonne Universités, UMR 7209, Archéozoologie, Archéobotanique: Sociétés, Pratiques et Environnements, 75005 Paris, France</footnote>

<footnote n="2"><sup>2</sup>Department of Archaeology, University of Aberdeen, St Mary's, Aberdeen, UK</footnote>

<footnote n="3"><sup>3</sup>UMR 7207 Centre de recherche sur la paléobiodiversité et les paléoenvironnements (CR2P), MNHN/CNRS/Univ. Paris 06, CP/38, 8 rue Buffon, 75005 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Centre for GeoGenetics, Natural History Museum of Denmark, University of Copenhagen, Øster Voldgade 5-7, 1350 K Copenhagen, Denmark</footnote>

<footnote n="5"><sup>5</sup>Laboratoire d'Anthropobiologie Moléculaire et d'Imagerie de Synthèse, Université de Toulouse, University Paul Sabatier, CNRS UMR 5288, 31000 Toulouse, France</footnote>
and temporal dynamics of the evolutionary history of the horse family, up to the latest domestication trajectories of horses and donkeys.

# 1. Introduction

Extant living species of equids are assigned to the genus *Equus* and are spread across large and diverse environments in Africa and Eurasia. The *Equus* genus is made up of two lineages: the caballine and non-caballine horses, and three phylogenetic clades. The first clade encompasses the caballine lineage that includes domesticated horses (*Equus ferus caballus*) and the endangered Przewalski's horse (*Equus ferus przewalskii*), which is the only remaining population of truly wild horses on the planet [1]. Zebras and wild asses form the two other remaining clades and belong to the non-caballine lineage. Zebras include three species plus a number of subspecies/morphotypes. Wild asses include the African wild ass (*Equus africanus*) and its domestic form the donkey (*Equus africanus asinus*), the Asian wild asses, *Equus kiang* and *E. hemionus*, the latter comprising several recognized subspecies [2–6]. This extant diversity is also made up of hybrids, as equine species can be crossbred with each other but remain sterile (although with some rare exceptions [7]). The most common hybrids are mules, a cross between a male donkey and a female horse, and to a lesser extent hinnies, the reciprocal mating.

The relatively limited biodiversity present in contemporary equine species appears in striking contrast with the past diversity of the Equidae family, in which several dozen genera of species have been identified from a fossil record spanning the last 55 Myr [8]. Most of this diversity emerged during an adaptive diversification stage which took place in the second half of the Miocene, following the increased range of new biomes and the transition from browsing to grazing [9]. The *Equus* genus, which represents the only genus of living equids, emerged approximately 4.0–4.5 Myr ago in Northern America, and then spread into the Old World 2.1–3.4 Myr ago [5,10] where the zebra and ass clades emerged. Caballine horses crossed the Bering Strait much later before expanding into the Old World; they became extinct in Northern America some 10 000 years ago [11].

Members of the *Equus* genus were abundant and diverse during the European Pleistocene [12], but almost became extinct at the beginning of the Holocene due to a combination of several possible factors, including climate change and overhunting. However, the *Equus* species have persisted in small ranges until now; with two of these being domesticated around 6000 years ago—in the west side of the Eurasian Steppe for the horse [13] and in Northern Africa for the donkey [14]—although the domestication processes of both [14–19] remain poorly understood.

Although the evolutionary history of the equid family provides one classical textbook example in palaeontology, our knowledge of a number of key evolutionary stages, such as the Plio-Pleistocene evolution of *Equus*, especially in America [20,21], is incomplete. Understanding how equine species spread through space and time requires the development of methods enabling robust identification from material preserved in the fossil record [22,23]. However, such identifications remain challenging for palaeontologists [8,24] and archaeologists [25], mostly due to the fragmentary nature of the material, convergent evolution and environmentally driven morphological changes, all of which limit the recoverable phylogenetic signal.

So far, ancient DNA analyses represent the ultimate method to secure an accurate taxonomic identification of ancient equids [1,5,6,20]. However, this approach is restricted to areas and/or periods conducive to the preservation of ancient DNA [26]; this would usually—but not always—exclude warm environments such as the deserts of the Middle East, Western Asia and Africa [27,28], though there are some exceptions [29,30]. Clearly, an integrative approach combining molecular and morphological proxies [31,32] would largely improve our understanding of the macro and micro evolutionary dynamics underlying equine speciation and domestication.

Teeth provide standard material to palaeontologists [33] and zooarchaeologists [34] for the reconstruction of past evolutionary changes. Teeth are well preserved in the fossil record and present relatively limited ecophenotypic plasticity, due to their genetically controlled development and absence of tissue remodelling once mineralized [35,36]. Crown morphology and the enamel pattern of mandibular cheek tooth occlusal surfaces are among the most popular morphological markers allowing palaeontologists and archaeologists to identify the taxonomy of Plio-Pleistocene equines [2,10,15,37,38], and reconstruct phylogenies and population histories [9,39]. The characters based on enamel patterns have long since been considered inappropriate to distinguish Pleistocene and Holocene horses [40], given that both age and wearing represent limiting factors for evolutionary [17] or taxonomic [25]
reconstructions. However, with recent advances in geometric morphometrics (GMM) in the field of biometry [41–44], it is now possible to capture and graphically visualize the complexity of tooth morphology. When applied to the cheek teeth of modern horses, GMM has demonstrated statistically that age and wearing have no significant effects on the overall occlusal enamel shape pattern when excluding juvenile and senile specimens [45]. In addition, the basis of the shape of this occlusal enamel pattern was found to support a reliable assignment of the various source populations [45], which pleads for further investigations aimed at reconstructing the history of horse breeding.

Here, we assess the relevance of occlusal enamel patterns in the lower cheek teeth of extant Old World equids, collected to provide reliable interspecific taxonomic and phylogenetic signals. A phylogenetic signal can be observed when related species show close phenotypic relationships [46–48]. However, to test this signal we have chosen to map the underlying morphospace onto the whole-exome based phylogeny presented by Jónsson and co-workers [5], following previous approaches for the projection of a phylogenetic tree into geometric morphometric space [49,50], and estimate the strength of the phylogenetic signal in the morphometric dataset using multivariate K statistic [51]. Our first objective was to propose a single GMM protocol to capture the form of the occlusal enamel patterns of the four mandibular cheek teeth (P/3, P/4, M/1, M/2) used in species identification of equid remains [15,37]. Key to our GMM protocol was to perform a reliable assessment of the anatomical tooth position in order to enable future applications to the vast collections of isolated teeth available from the fossil record and museum collections. Our second objective was to find which mandibular cheek tooth captured the strongest interspecific differences and phylogenetic signal. Using dental phenotypes as reliable systematic and phylogenetic markers in equids would open up an incredibly large range of future investigations. Firstly, of the diversity of Plio–Pleistocene horses and, secondly, of the Holocene diversity of southwest Asian equids; where all three extant clades plus their hybrids and extinct forms could have occurred [29,52,53], and from which the wild ancestors of domesticated horses and donkeys emerged.

## 2. Material and methods

### 2.1. Material

We studied a total of 345 specimens from 15 equine taxa covering the current distribution of all living members of the *Equus* genus (figure 1), including two taxa of horses, two taxa of African asses, six taxa of Asian asses, four taxa of zebras and one group of ass/horse hybrids (mules) (table 1; electronic supplementary material, S1). For each individual specimen, the occlusal view of the entire mandibular tooth row was photographed (electronic supplementary material, S1, details can be found also on www.vera-eisenmann.com). The occlusal surface of the mandibular cheek teeth were photographed using several film cameras (i.e. Hasselblad, Canon and Nikon), ensuring that the occlusal surface was parallel to the camera focal place. All the argentic photographs were scanned with a 300 dpi resolution, using a Toshiba e-studio 2330c scanner. Some of the specimens photographed within this project, were from the comparative anatomy collection of the National Museum of Natural History of Paris, using a Canon EOS 70D with Canon macro lens EF-S 60 mm.

### 2.2. Geometric morphometrics

From the picture of the mandibular tooth row from each specimen, the occlusal enamel patterns of P/3, P/4, M/1 and M/2 were recorded using the same novel two-dimensional protocol. This approach represents a significant increment compared with our previous study of the enamel pattern in equine mandibular teeth, which was based on a basic approach for capturing the complexity of the enamel folds with few landmarks on the maximum curvature of the enamel loops [45]. This study aimed to capture the entire complexity of the enamel folding by extracting as much biological information as possible (figure 2). More specifically, we used eight landmarks and eight curves, between each landmark, with a total of 178 equidistant semi-landmarks. Each curve was drawn in the centre of the enamel thickness (see details of the number of points per curves in figure 2). Only the left mandibles and teeth were used; all pictures of right mandibles were mirrored left.

Landmarks and semi-landmarks were digitized from images with a TPSdig2 v. 2.17 [54].

To be statistically comparable, enamel outlines must be standardized by a generalized Procrustes superimposition (GPS), which aligns their coordinates with a least squared approach in order to remove all information of position, scale and orientation from the initial points' configurations. In contrast with landmarks, semi-landmarks do not have an exact correspondence on the enamel curve. During
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Distribution of extant clades of equids in Eurasia and Africa.</figcaption></figure>

<table><caption>**Table 1.** Equid taxa studied with their abbreviation code and their sample size per each mandibular cheek tooth.</caption><thead><tr><th>taxa</th><th>vernacular names</th><th>code</th><th>P/3</th><th>P/4</th><th>M/1</th><th>M/2</th></tr></thead><tbody><tr><td><i>E. ferus caballus</i></td><td>domestic horse</td><td>CBL</td><td>15</td><td>15</td><td>15</td><td>15</td></tr><tr><td><i>E. f. przewalskii</i></td><td>Przewalski's horse</td><td>PRZ</td><td>9</td><td>8</td><td>8</td><td>8</td></tr><tr><td><i>E. africanus somaliensis</i></td><td>Somali wild ass</td><td>AFR</td><td>7</td><td>7</td><td>7</td><td>7</td></tr><tr><td><i>E. a. asinus</i></td><td>domestic donkey</td><td>ASN</td><td>11</td><td>11</td><td>10</td><td>10</td></tr><tr><td><i>E. kiang</i></td><td>Tibetan kiang</td><td>KNG</td><td>12</td><td>13</td><td>13</td><td>13</td></tr><tr><td><i>E. hemionus hemionus</i></td><td>Mongolian wild ass</td><td>MNG</td><td>13</td><td>10</td><td>13</td><td>13</td></tr><tr><td><i>E. h. onager</i></td><td>onager</td><td>ONG</td><td>7</td><td>8</td><td>7</td><td>7</td></tr><tr><td><i>E. h. khur</i></td><td>Indian wild ass</td><td>KHR</td><td>7</td><td>8</td><td>8</td><td>6</td></tr><tr><td><i>E. h. kulan</i></td><td>kulan</td><td>KLN</td><td>9</td><td>12</td><td>10</td><td>11</td></tr><tr><td><i>E. h. hemippus</i></td><td>Syrian wild ass</td><td>HMP</td><td>2</td><td>3</td><td>2</td><td>3</td></tr><tr><td><i>E. grevyi</i></td><td>Grévy's zebra</td><td>GRV</td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><td><i>E. zebra hartmannae</i></td><td>Hartmanns' mountain zebra</td><td>ZBR</td><td>4</td><td>4</td><td>4</td><td>5</td></tr><tr><td><i>E. quagga quagga</i></td><td>Quagga plain zebra</td><td>QGA</td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><td><i>E. q. burchellii</i></td><td>Burchell's zebra</td><td>BRC</td><td>8</td><td>8</td><td>8</td><td>8</td></tr><tr><td>hybrid</td><td>donkey × horse</td><td>HBR</td><td>4</td><td>4</td><td>4</td><td>4</td></tr></tbody></table>

the superimposition process, they are slid along the tangent of the curves using the bending energy algorithm [55] to minimize the difference between the shape configurations [56]. The by-product of this superimposition is a new set of Procrustes coordinates (or shape coordinates) used for statistical analyses. The GPS was performed with a TPSrelw v. 1.54 [57]. We provide the full morphometric dataset as the electronic supplementary material, including specimen ID, grouping factors and Procrustes coordinates (electronic supplementary material, S2).
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (a) Position of the four cheek teeth in the *Equus* mandible with (b) their anatomical terminology in occlusal view from [15] and (c) their common GMM protocol for the quantification of the enamel folding with the location of the 8 landmarks (dots with grey filling) and the number of semilandmarks (table) sampling the curves of the enamel folds between each landmark.</figcaption></figure>

## 2.3. Statistical analyses

### 2.3.1. Anatomical and taxonomic comparisons

To assess whether the shape difference among the four mandibular teeth and all equine taxa were significant, we initially used a factorial MANOVA to account for possible interactions between both anatomical and taxonomic factors. Detecting significant interaction, between the anatomical position of the molar in the mandibular sequence and the taxonomic identification of the specimens, would mean that the taxonomic differences varied according to the anatomical position of the tooth; therefore, discriminant analysis within taxonomic groups had to be performed separately for each tooth.

To assess the assignment accuracy among the four mandibular teeth and the 15 taxa, considered using the shape variation of the enamel outline, we performed a canonical variate analysis (CVA) with a twofold cross-validated classification. This approach calculates the most discriminant axes of variation among the taxonomic grouping factors, following normalization of the intra-group variance by creating a discriminant shape space. Note that the latter is different from the original shape space obtained from a principal component analysis (PCA) which calculates the main axis of variation in the shape coordinates without normalization.

To estimate the performance of the anatomical and taxonomic discriminant models in practice, and their ability to identify the right anatomical position or the right taxonomic group of an isolated tooth, we calculated the percentage of correct classification of each specimen into its known group, be it anatomical or taxonomic, following a twofold cross-validation approach. This consisted of partitioning the dataset into two subsets: the training set on which the analysis was performed, and the testing set that was used to validate the analysis. This procedure is done over multiple rounds (10 000) of cross-validation, changing the dataset partitioning at every round. This provided an average of correct classification for each specimen and a percentage of correct classification, or correct cross-validation (CCV).

To reduce the dimensionality of our shape dataset, we performed a PCA on the shape coordinates before each CVA analysis, and retained the PC scores that explained 99% of the total shape variance.

To visualize shape changes along the main axis of variation (PCA) or along the discriminant axes (CVA), we used a multivariate regression approach [58].

All these analyses were performed with R [54], using the libraries ade4 [59] and geomorph [60]. Graphical displays were performed with MorphoJ [61].
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Lineage divergence in equids including the numbering of the internal nodes based on the exome-wide phylogeny [5].</figcaption></figure>

### 2.3.2. Phylogenetic signal in enamel outline shape

We used the published molecular phylogeny of equids (figure 3; electronic supplementary material, S3) built for 10 taxa: two horses (E. ferus przewalskii, E. f. caballus), four zebras (E. quagga quagga, E. q. burchellii, E. grevyi, E. zebra hartmannae), two asses (E. africanus somaliensis, E. a. asinus) and two hemiones (E. h. onager, E. kiang). This phylogeny was obtained applying a maximum-likelihood method on a super-matrix of 20374 protein-coding genes, each partitioned in 1st + 2nd and 3rd codon position (representing a total of 31219584 sites). Node robustness was assessed using 100 bootstrap pseudo-replicates. The analysis was performed using the phylogenetic module of the automated PALEOMIX pipeline [62].

To assess the association between the molecular phylogeny and the dental shape space, we created a phylomorphospace which fitted the shape changes of each tooth from the 10 species onto the tree topology of the molecular phylogeny (including branch length) [49,50]. A dental shape space was constructed for each mandibular tooth of the 10 species using a GPS, followed by the computation of the species mean shape. A PCA was subsequently performed on the matrix of variance covariance of the averaged Procrustes coordinates. To fit the dental shape space onto the phylogeny, we used the squared-change parsimony criterion [63], which finds the PC scores of the shape of the ancestral (root) and internal nodes, for which the changes along the branches of the entire phylogeny (in Procrustes distances units) are minimal [49]. After which, the branches of phylogeny—connecting the ancestral, internal and terminal nodes—can be drawn to estimate the path of the evolutionary lineages in the shape space. The squared-change parsimony approach, which produces the same ancestral shape estimates as maximum-likelihood phylogenetic reconstruction and phylogenetic generalized linear models (PGLM) [64], allows the reconstruction of the ancestral shape at the root of the phylogeny from the shape average of the terminal nodes. It is then possible to graphically visualize the evolutionary changes of the teeth shape from the common ancestor of the *Equus* species, and each internal node and terminal tip of the phylogeny.

To assess the phylogenetic signal in the shape space of the four teeth, we used the *Kmult* method [51]; which is more suited to high-dimensional multivariate data and less affected by the difference of traits variation among species [65]. This method is a multivariate generalization of the *K* statistic which estimates the strength of the phylogenetic signal in a univariate trait relative to a Brownian motion model of evolution [65]. Under a Brownian motion model of evolution the expected *Kmult* value is 1.0. Observed *Kmult* values are evaluated by comparison with the permutated *K* value from randomized shape data relative to the tree topology of the molecular phylogeny.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Canonical variate analyses (CVA) comparing shapes P/3, P/4, M/1 and M/2 in equids with visualization of shape differences along the canonical axes and shape reconstruction at the extreme values of each axis.</figcaption></figure>

<table><caption>**Table 2.** Results of the factorial MANOVA which tested the differences between the four mandibular cheek teeth (anatomical position) and the 15 taxa of our dataset (table 1), as well as the interaction between the two factors (anatomical position : taxa).</caption><thead><tr><th></th><th>d.f.</th><th>SS</th><th>MS</th><th>Rsq</th><th>F</th><th>Z</th><th>Pr (&gt; F)</th></tr></thead><tbody><tr><td>anatomical position</td><td>3</td><td>1.3784</td><td>0.45948</td><td>0.13564</td><td>41.2</td><td>18.452</td><td>0.001**</td></tr><tr><td>taxa</td><td>14</td><td>2.2202</td><td>0.15859</td><td>0.21848</td><td>14.2198</td><td>7.7823</td><td>0.001**</td></tr><tr><td>anatomical position : taxa</td><td>42</td><td>1.4222</td><td>0.03386</td><td>0.13995</td><td>3.0362</td><td>1.7319</td><td>0.001**</td></tr></tbody></table>

## 3. Results

The factorial MANOVA found significant differences in shape among the four cheek teeth and the 15 equine taxa investigated (table 2). Although the taxonomic signal was stronger than the anatomical signal, with 21.8% against 13.5% of shape variation explained, respectively, significant interaction between both factors was found to explain 14% of the overall variation. This indicates that the taxonomic differences vary across the four cheek teeth and therefore a separate investigation of the taxonomic resolution for each tooth is required.

### 3.1. Anatomical distinction of equine mandibular cheek teeth excluding P/2 and M/3

The discriminant model (figure 4) found correct anatomical identification for each mandibular tooth across the 15 taxa dataset for 93% of the classifications after cross-validation. This shows that this model can accurately identify the anatomical position of any isolated mandibular tooth of the genus *Equus* from the archaeological record with a false-positive probability of only 7%. Misclassifications were systematically between M/1 and M/2 (37%) and between P/3 and P/4 (49%). The remaining 14% were random.

The visualization of the shape changes along canonical axes 1 and 2 (figure 4) showed that the main differences between premolars and molars were a more penetrating ectoflexid, and consequently, a less wide postflexid in the latter. In addition, the premolars displayed flatter proto- and hypoconids as well as a rounder hypoconulid than the molars.

### 3.2. Taxonomic distinction among equids

The discriminant model of the CVA shows that all four mandibular teeth could distinguish the 15 equine taxa of our dataset (table 1), with more than 95% of correct classification with cross-validation (table 3). Premolars performed better than molars, with P/3 being the most efficient with a 100% cross-validated classification (CCV).

The pooled within-group shape space of the 15 taxa studied (figure 5) shows a clear taxonomic signal separating the three clades of zebras, asses and horses. However, the hemippe, also known as
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Morphospace for the four mandibular equine cheek teeth. The shape changes of the enamel folding along the major axes of variations (PC1 and PC2) are displayed for the extreme values of each axis. The colour of the dots corresponds to the molecular clades (blue: caballine horses, red: asses, orange: hemiones, green: zebras) and hybrids (grey filling). The three letter abbreviations correspond to the taxonomic groups described in table 1.</figcaption></figure>

<table><caption>**Table 3.** Results using MANOVA of the differences between the 15 taxa for each mandibular cheek tooth and the percentage of correct cross-validated classification (CCV) for each mandibular tooth.</caption><thead><tr><th></th><th>Pillai</th><th>Approx F</th><th>num d.f.</th><th>den d.f.</th><th>Pr</th><th>CCV</th></tr></thead><tbody><tr><td>P/3</td><td>12.439</td><td>8.09</td><td>896</td><td>910</td><td>&lt;0.0001</td><td>100</td></tr><tr><td>P/4</td><td>12.581</td><td>9.14</td><td>896</td><td>924</td><td>&lt;0.0001</td><td>98.47</td></tr><tr><td>M/1</td><td>12.186</td><td>8.08</td><td>826</td><td>994</td><td>&lt;0.0001</td><td>95.41</td></tr><tr><td>M/2</td><td>12.188</td><td>7.38</td><td>854</td><td>938</td><td>&lt;0.0001</td><td>96.12</td></tr></tbody></table>

the Syrian wild ass (*E. hemionus hemippus*), is divergent from other hemiones in all teeth but the M/2. Mules (donkey × horse hybrids) have a divergent phenotype from horses and asses, except for the P/4 where mules are within the asinine dental variation. The divergence of mules is almost intermediate
<table><caption>**Table 4.** Phylogenetic signal from the four mandibular cheek teeth of equids using the *Kmult* method [55].</caption><thead><tr><th></th><th><em>Kmult</em> obs</th><th>Permut <em>p</em>-value</th></tr></thead><tbody><tr><td>P/3</td><td>0.432</td><td>0.002</td></tr><tr><td>P/4</td><td>0.4253</td><td>0.001</td></tr><tr><td>M/1</td><td>0.4905</td><td>0.001</td></tr><tr><td>M/2</td><td>0.5711</td><td>0.001</td></tr></tbody></table>

between horses and wild asses in molars, which strongly suggests that mules could be discriminated in the archaeological record with a statistical predictive approach, should greater sample sizes of mules be available. Such methodological approach could greatly contribute to document the role played by equine hybrids in human history, especially during the Roman and medieval periods [66].

The shape deformation along the main axis of shape variation illustrates how different the morphological criteria are along the location sequence in the tooth row, as previously evidenced by the factorial MANOVA.

(1) For the P/3, the width and shape of the postflexid separates zebras and horses from hemiones and the African wild ass and donkey; while the smoother V-shape of the linguaflexid (double knot) separates horses from zebras, except for *E. z. hartmannae*, which shows a more asinine shape.

(2) For the P/4, the penetration of the ectoflexid is deeper in zebras than in wild and domestic asses and, to a lesser extent, in horses; except for the *E. z. hartmannae*, which again shows a more asinine shape.

(3) For the M/1, the penetration of the ectoflexid up to the linguaflexid fold and flatter protoconid and hypoconid separates zebras and horses from asses; except for the *E. z. hartmannae*, which again shows a more asinine shape. Along the second PC, the kiang (*E. kiang*), the khur (*E. hemionus khur*) and the Mongolian hemione (*E. h. hemionus*) differ from asinines by slender protoconids and hypoconids.

(4) For the M/2, clearer criteria emerge due to a greater structuration of the shape space among the taxa, due to the lack of divergence of the hemippe. Similar to the M/1, the combination of the penetration of the ectoflexid up to the linguaflexid fold and flatter protoconid and hypoconid separates zebras from asses. The latter also shows a wider postflexid up to the preflexid, clearly separating it from zebras and horses. Confirming previous observations from Véra Eisenmann [15], the shape of the linguaflexid can provide criteria to distinguish zebras, horses, asses and hemiones. Asses and hemiones display a V-shaped linguaflexid compared to the more U-shape of horse. Zebras are more intermediate with a smoother V-shape. Horses can be distinguished from zebras due to their less penetrating ectoflexid and their wider postflexid.

### 3.3. Phylogenetic signal in equine mandibular cheek teeth (P/3 to M/2)

The *K* statistic found significant phylogenetic signals in all four mandibular teeth (table 4 and figure 6a). The M/2 has the highest *Kmult* close to 0.6, which suggests that the M/2 shape variation provides the greatest degree of phylogenetic structure among the 10 *Equus* taxa studied. However, the level of phylogenetic signal is still below the *K* value of 1, which suggests that the M/2 shape evolution does not fit an unconstrained Brownian model of evolutionary process.

Nevertheless, the trend of phylogenetic signals from *Equus* mandibular teeth is supported by the phylomorphospace (figure 6b); with related species showing phenotypic similarity to their dental mean shapes, while more distant species are more morphologically divergent. In addition, the dental shape variation seems to diversify from the hypothesized ancestral shape at the root; with extant species occupying a distant region of the phylomorphospace compared with the central ancestral shape, especially for the M/2. Among the four cheek teeth, the M/2 clearly provides the tightest association with the main phylogenetic clades closely related in the shape space with the least discrepancy from the phylogenetic tree (figure 6b). Finally, the lack of overlapping branches in the phylomorphospace suggests that convergent evolution did not occur during the evolution of the mandibular teeth. However, aside from the M/2, all the cheek teeth of the *E. z. hartmannae* cluster towards the asinine group; though some signs of convergence among asses and hemiones are observable in the P/3 and M/1.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

**Figure 6.** (a) Histograms of the *Kmult* values obtained from 1000 permutations of cheek tooth shapes from the tips of the molecular phylogeny with the position of the observed value identified by a red arrow. (b) Plot of the phylomorphospace viewed from the first two PCs for the mean shapes of 10 equine taxa onto the tree topology of the molecular phylogeny with a squared-change parsimony approach.

The first two PCs on the M/2 shape account for almost 75% of the total variation for the 10 species. The reconstruction of the evolutionary changes of the M/2 shape on the exome phylogeny (figure 7) shows a clear divergence between the three clades. The main divergence is in the direction of the PC1 that separates zebras and asinines through the penetration of the ectoflexid of zebras, and the V-shaped linguaflexid and the shoe-shaped entoflexid in asinines. The second divergence observed along the PC2 is the separation of the horse and the diversification among zebras and asses; with the divergence of *E. z. hartmannae* from its sister species and the separation of the African and Asian asinine clades, respectively. The diversification within asinines is marked by a more pronounced V-shaped linguaflexid
<figure><img src="image_18.png" /><figcaption>**Figure 7.** Reconstruction of the evolutionary changes in the shape of M/2 equine enamel folding patterns. The phylogenetic tree from the exome-wide phylogeny has been projected onto the shape space made up of the first two principal components computed from the matrix of variance covariance among the 10 species means. The tips of the branches are the position of the species means in the shape space. The diagrams at the tip of the branches display the changes in the shape of a species (black outline) respective to the reconstructed common ancestor at the root (grey outline), with a scale factor of 1.0.</figcaption></figure>

in asses. Horses are close to the reconstructed common ancestral shape at the root of the tree; their divergence is marked mainly by the greater U-shape of the linguaflexid.

## 4. Discussion

Traditional morphological analyses have so far provided invaluable information to palaeontologists and zooarchaeologists for documenting the evolutionary history of equids, and remain key for future investigation. However, these approaches face two main limitations in their ability to capture the entire geometry of equine teeth and to disentangle the size from the shape component of their form variation. As the environmental impacts on the size variation of the skeleton in mammals [67], it can also potentially confound the phylogenetic and population signals retrieved from morphological data [48,68]. GMM approaches have contributed to resolve these long-lasting caveats by enabling the separate extraction of the size and shape components of biological forms such as teeth, redeeming the use of morphological variation in synergy with palaeogenetic analyses to investigate macro- and microevolution processes in the palaeontological and archaeological records [32,34,46].

Our previous GMM approach of horse upper cheek tooth morphology had already shown that the shape of the occlusal enamel pattern was a reliable character to discriminate populations [45]. It also clearly evidenced the absence of an age-related effect in the shape variation of the occlusal enamel pattern, except for senile and juvenile specimens [45], contradicting previous morphometric studies which considered teeth to be too much affected by wear patterns to provide a reliable discriminant marker between species [40,69]. Here, we improved the quantification of the complex enamel folding patterns in equine mandibular cheek teeth in order to fully explore their potential as taxonomic and phylogenetic markers for future palaeontological and archaeological investigations.

Our GMM approach takes into account the heterogeneous nature of the fossil record, which often consists of isolated teeth. Our discriminant model computed over 15 equine taxa and successfully separated the dental sequence of the four mandibular teeth, with a probability above 90%. This means that the anatomical position in the dental sequence of any archaeological isolated mandibular tooth, except the P/2 and M/3 can now be assessed accurately. In practice, this can be done through a CVA with a predictive approach [70], which aims to assign a set of undetermined objects to a series of four
predefined groups (P/3, P/4, M/1, M/2). This approach will compute, initially, the multiple discriminant or canonical functions from the samples defined by the four anatomical groups and will, secondly, apply these canonical functions to the undetermined isolated mandibular teeth from the fossil record. The assignment to their most likely anatomical groups can then be assessed based on their Mahalanobis distances and the *ad hoc* probabilities associated with the four anatomical groups' centroid.

All four mandibular cheek teeth investigated in our study have shown great taxonomic accuracy at the species level, especially the P/3. This suggests that there is a strong taxonomic pattern in the shape of the enamel folding of equids. However, we found that anatomical and taxonomic factors interact, which clearly stresses the fact that each location in the dental sequence will give different phenotypic criteria to identify the species (figure 4). For example, the criterion of the penetration of the ectoflexid up to the linguaflexid folds to separate zebras from asses and horses is only applicable to the molars and not the premolars, supporting previous morphoscopic observations [37]. For the P/3 and P/4, the main criteria separating zebras from asses and horses is the shape of the entoflexid, which appears longer and thinner in zebras and horses than in asses. Again, the U-shape of the linguaflexid to separate the *E. ferus caballus* from zebras, asses and hemiones is mainly visible on the M/2; and appears closer to a smooth V-shape than a true U-shape, which is not obvious for the other mandibular teeth. This dental-specific taxonomic signature probably contributes to the difficulty to perform accurate species assignment when dealing with isolated mandibular teeth [25]. This significant result is specifically useful for archaeozoologists, as it can be concluded that the criteria used previously with the shape of the double knot is valid on the M/2 but not for the premolars.

Among the four mandibular cheek teeth, the P/3 provides complete accuracy in the identification of the 15 equine taxa. However, the M/2 morphospace provides the clearest taxonomic clustering of extant equids, probably because most of the shape space of the other teeth is driven by the divergence of the Syrian wild ass (*E. hemionus hemippus*). The divergent position of hybrid (mules) dental phenotype, between horse/zebra and the asinines phenotypes, has been observed in molars but not premolars. This suggests that phenotypic divergence in equid molars is more congruent with molecular divergence than premolars. The M/2 also provides the strongest phylogenetic signal and the best inference of the equine phylogeny, based on exome sequence data [5]. This demonstrates for the first time that the shape of enamel folds in mandibular teeth encapsulates the phylogenetic relationships within the genus *Equus*, and that related species of *Equus* share similar dental phenotype. However, the evolutionary processes underlying this phylogenetic signal are potentially manifold [71]. The diversification of the related *Equus* taxa from a central common ancestor in a different corner of the phylomorphospace of the M/2 (figure 6) suggests that the evolutionary changes of the enamel fold shape from the 4.5 Myr- old most recent common ancestor of modern *Equus* until the present day, can chiefly be explained by a vicariant speciation during the radiation from the Americas into Eurasia and Africa. It also suggests that most of the morpho-functional changes in the tooth morphology may have happened earlier in equine evolution, probably during the Early–Middle Miocene (20–15 Ma) when equids adapted their tooth morphology to grazing in extensive grasslands of open biomes [9]. The power of GMM for taxonomic discrimination of early equids based on cheek teeth remains, however, to be demonstrated. GMM has proved to be suitable to address non-dental evolution of early taxa of the family [72], but when applied to cheek teeth of crown equids such an approach may be less discriminant due to their less divergent morphology.

However, even the highest *K* value found in an M/2 is below the threshold of a Brownian motion model of evolution; which suggests that genetic drift and randomly fluctuating natural selection over time does not fully account for the dental shape diversification in the *Equus* and that evolutionary conservatism [73,74] in dental variation should be rejected. This absence of complete hierarchical phylogenetic relationships in the dental shape variation suggests that other factors, such as morphological convergence and interspecific hybridization, could have occurred [65], as also suggested from patterns of genomic variation [9]. The morphological proximities observed among distant species such as asses and hemiones or between the mountain zebra (*E. z. hartmannea*) and asinines on all but the M/2 cheek teeth could support this possibility, without disentangling convergence from interspecific hybridization as a causal factor. But the latest study on the speciation process in *Equus* is based on an extensive genome-wide dataset and has provided strong evidence for a speciation gene flow in equids [10]. These authors found four, yet undated, events of gene flow among equine lineages following their divergence: between horses and zebra/ass clades, between kiang and donkey, between the African wild ass and Grévy's zebra, and between asses and mountain zebra. This clearly supports gene flow as one of the factors explaining the incomplete hierarchical phylogenetic relationships in *Equus* dental shape variation.
## 5. Conclusion

This comparative study provides clear evidence for the relevance of dental phenotypes to accurately discriminate between all modern members of the genus *Equus* and capture their phylogenetic relationships. Our work redeems the use of dental morphology to explore, in the fossil record, the spatial and temporal dynamics of the equine evolutionary history during the Plio-Pleistocene; and opens up total evidence approaches integrating fossils with molecular phylogeny to improve our understanding of their phenotypic evolution. In the Late Pleistocene and Holocene fossil record of the Old World, the implementation of such a comparative approach will provide an unprecedented opportunity to explore the diversity of the wild ancestors of the domestic horse and donkey, and how these relate to the earliest archaeological evidence of domestic forms.

Data accessibility. The data supporting this article are available in the electronic supplementary material.

Authors' contributions. T.C. and M.M. conceived the study together. A.M. and M.M. performed the data collection in collaboration with K.D. T.C. performed the data analysis and wrote the manuscript with input from A.M., L.O., S.P. and M.M. All authors gave final approval for publication.

Competing interests. The authors have no competing interests.

Funding. The Morphometric Geometric work was supported by the National Museum of Natural History 2014 ATM grant 'Formes possibles, Formes réalisées' led by M.M. and T.C. with contributions from A.M. The genetic work was supported by the Danish National Research Foundation (grant no. DNRF94); the ERC-Consolidator Grant PEGASUS (Grant 681605) and the 'Initiative d'Excellence Chaires d'attractivité', Université de Toulouse (OURASI).

Acknowledgements. This work would not have been possible without the generous help of Dr Véra Eisenmann, who authorized the use of her photographic documentation and offered insightful comments on the initial steps of this study. We thank Shiva Sheikhi (UMR 7209) for the scanning of this documentation. We also thank J. Lesur, curator of the Mammalian osteological collections at the National Museum of Natural History of Paris (MNHN) for access to the equid collections. We are very grateful to David Polly and three anonymous reviewers for their constructive critiques that have greatly improved the quality of the manuscript.

## References

1. Der Sarkissian C et al. 2015 Mitochondrial genomes reveal the extinct *Hippidion* as an outgroup to all living equids. *Biol. Lett.* **11**, 20141058. (doi:10.1098/rsl.2014.1058)

2. Eisenmann V. 1981 Etude des dents jugales inférieures des *Equus* (Mammalia, Perissodactyla) actuels et fossiles. *Palaeovertebrata* **10**, 127–226.

3. Eisenmann V, Turlot J-C. 1978 Sur la taxinomie du genre *Equus* (équidés). Description et discrimination des espèces d'après les données crâniométriques. *Les Cahiers de l'Analyse des Données*, **3**, 179–201.

4. Eisenmann V, Mashkour M. 1999 The small equids of Binagady (Azerbaidjan) and Qazvin (Iran): *E. hemionus binagadensis* nov. subsp. and *E. hydruinus. Geobios* **32**, 105–122. (doi:10.1016/S0016-6995(99)80085-0)

5. Jónsson H et al. 2014 Speciation with gene flow in equids despite extensive chromosomal plasticity. *Proc. Natl Acad. Sci. USA.* **111**, 18 655–18 660. (doi:10.1073/pnas.1412627111)

6. Vilstrup JT et al. 2013 Mitochondrial phylogenomics of modern and ancient equids. *PLoS ONE* **8**, e55950. (doi:10.1371/journal.pone.0055950)

7. Steiner CC, Ryder OA. 2013 Characterization of *Prdm9* in equids and sterility in mules. *PLoS ONE* **8**, e61746. (doi:10.1371/journal.pone.0061746)

8. MacFadden BJ. 1995 Fossil horses: systematics, paleobiology, and evolution of the family Equidae. *Terra Nova* **7**, 102–103. (doi:10.1111/j.1365-3121.1995.tb00674.x)

9. Macfadden BJ. 2005 Evolution. Fossil horses—evidence for evolution. *Science* **307**, 1728–1730. (doi:10.1126/science.1105458)

10. Eisenmann V. 1979 Caractères évolutifs et phylogénie du genre *Equus* (Mammalia, Perissodactyla). *C. R. Acad. Sci. Paris* **288**, 497–500.

11. Haile J et al. 2009 Ancient DNA reveals late survival of mammoth and horse in interior Alaska. *Proc. Natl Acad. Sci. USA* **106**, 22 352–22 357. (doi:10.1073/pnas.0912510106)

12. Forsten A. 1989 Horse diversity through the ages. *Biol. Rev.* **64**, 279–304. (doi:10.1111/j.1469-185X.1989.tb00677.x)

13. Outram AK, Stear NA, Bendrey R, Olsen S, Kasparov A, Zaibert V, Thorpe N, Evershed RP. 2009 The earliest horse harnessing and milking. *Science* **323**, 1332–1335. (doi:10.1126/science.1168594)

14. Kimura B, Marshall F, Beja-Pereira A, Mulligan C. 2013 Donkey domestication. *Afr. Archaeol. Rev.* **30**, 83–95. (doi:10.1007/s10437-012-9126-8)

15. Eisenmann V. 1986 Comparative osteology of modern and fossil horses, half-asses and asses. In *Equids in the ancient world* (ed. RHM Uerpmann), pp. 67–106. Wiesbaden, Germany: Dr. Ludwig Reichert Verlag.

16. Mashkour M. (ed.). 2006 Equids in time and space. In *Papers in honour of Vera Eisenmann (Proceedings of the 9th ICAZ Conference, Durham, 2002)*, 240. Oxford, UK: Oxbow Books.

17. Olsen SL (ed.). 1996 *Horses through time*. Boulder, CO: Roberts Rinehart Publishers for Carnegie Museum of Natural History.

18. Beja-Pereira A et al. 2004 African origins of the domestic donkey. *Science* **304**, 1781. (doi:10.1126/science.1096008)

19. Bendrey R. 2012 From wild horses to domestic horses: a European perspective. *World Archaeol.*

44, 135–157. (doi:10.1080/00438243.2012.647571)

20. Weinstock J et al. 2005 Evolution, systematics, and phylogeography of pleistocene horses in the new world: a molecular perspective. *PLoS Biol.* **3**, e241. (doi:10.1371/journal.pbio.0030241)

21. Eisenmann V, Howe J, Pichardo M. 2008 Old world hemiones and new world slender species (Mammalia, Equidae). *Palaeovertebrata* **36**, 159–233. (doi:10.18563/pv.36.1-4.159-233)

22. Eisenmann V. 1996 Quaternary horses: possible candidates to domestication. In *Proc. XII Cong. Forlì-Italia 1996*, 8–14 September (ed. Collectif), pp. 27–36. Forlì, Italy: A.B.A.C.O. Edizioni.

23. Eisenmann V, Alberdi MT, Giuli (de) C, Staesche U. 1988 Volume I: Methodology. In *Studying Fossil Horses: New York Int. Hipparion Conf.*, 1981 (eds MO Woodburne, P Sondaar), pp. 1–71. Leiden, The Netherlands: E.J. Brill.

24. Eisenmann V, Baylac M. 2000 Extant and fossil *Equus* (Mammalia, Perissodactyla) skulls: a morphometric definition of the subgenus *Equus*. *Zool. Scripta* **29**, 89–100. (doi:10.1046/j.1463-6409.2000.00034.x)

25. Twiss KC, Wolfhagen J, Madgwick R, Foster H, Demirergi GA, Russell N, Everhart JL, Pearson J, Mulville J. 2016 Horses, hemiones, hydruntines? Assessing the reliability of dental criteria for assigning species to southwest Asian equid remains. *Int J Osteoarchaeol*. (Early view) (doi:10.1002/oa.2524)

26. Orlando L et al. 2013 Recalibrating *Equus* evolution using the genome sequence of an early Middle
Pleistocene horse. *Nature* **499**, 74–78. (doi:10.1038/
nature12323)

27. Smith CI, Chamberlain AT, Riley MS, Stringer C, Collins MJ. 2003 The thermal history of human fossils and the likelihood of successful DNA amplification. *J. Hum. Evol.* **45**, 203–217. (doi:10.1016/S0047-2484(03)00106-4)

28. Bollongino R, Vigne J-D. 2008 Temperature monitoring in archaeological animal bone samples in the Near East arid area, before, during and after excavation. *J. Archaeol. Sci.* **35**, 873–881. (doi:10.1016/j.jas.2007.06.023)

29. Orlando L, Mashkour M, Burke A, Douady C, Eisenmann V, Hänni C. 2006 Geographic distribution of an extinct equid (*Equus hydruntinus*: Mammalia, Equidae) revealed by morphological and genetical analyses of fossils. *Mol. Ecol.* **15**, 2083–2093. (doi:10.1111/j.1365-294X.2006.02922.x)

30. Broushaki F *et al.* 2016 Early Neolithic genomes from the eastern Fertile Crescent. *Science* **353**, 499–503. (doi:10.1126/science.aaf7943)

31. Slater GJ, Harmon LJ, Alfaro ME. 2012 Integrating fossils with molecular phylogenies improves inference of trait evolution. *Evolution* **66**, 3931–3944. (doi:10.1111/j.1558-5646.2012.01723.x)

32. Jablonski D, Shubin NH. 2015 The future of the fossil record: paleontology in the 21st century. *Proc. Natl Acad. Sci. USA* **112**, 4852–4858. (doi:10.1073/pnas. 1505146112)

33. Polly PD. 2003 Paleophylogeography: the tempo of geographic differentiation in marmots (*Marmota*). *J. Mammal.* **84**, 369–384. (doi:10.1644/1545-1542 (2003)084–0369:PTT0GD>–2.0.CO;2)

34. Cucchi T *et al.* 2014 The changing pace of insular life: 5000 years of microevolution in the Orkney vole (*Microtus arvalis* orca denis). *Evolution* **68**, 2804–2820. (doi:10.1111/evo.12476)

35. Jernvall J, Kangas A, Kavanagh K, Salazar-Ciudad I. 2004 Mammalian dental diversity: so many shapes, so few genes. *Integr. Comp. Biol.* **44**, 577.

36. Jernvall J, Keränen SVE, Thesleff I. 2000 Evolutionary modification of development in mammalian teeth: quantifying gene expression patterns and topography. *Proc. Natl Acad. Sci. USA* **97**, 1444–1448. (doi:10.1073/pnas.97.26. 14444)

37. Davis S. 1980 Late Pleistocene and Holocene equid remains from Israel. *Zool. J. Linn. Soc.* **70**, 289–312. (doi:10.1111/j.1096-3642.1980.tb00854.x)

38. Hopwood AT. 1936 The former distribution of caballine and zebrine horses in Europe and Asia. *Proc. Zool. Soc. Lond.* **106**, 897–912. (doi:10.1111/ j.1469-7998.1936.tb06291.x)

39. McFadden BJ. 1994 *Fossil horses: systematics, paleobiology, and evolution of the family Equidae*. Cambridge, UK: Cambridge University Press.

40. Forsten A. 1988 The small caballoid horse of the upper Pleistocene and Holocene. *J. Anim. Breed. Genet.* **105**, 161–176. (doi:10.1111/j.1439-0388.1988. tb00288.x)

41. Adams DC, Rohlf FJ, Slice DE. 2003 Geometric morphometrics: ten years of progress following the 'revolution'. *Ital. J. Zool.* **71**, 5–16. (doi:10.1080/ 11250000409356545)

42. Adams DC, Rohlf FJ, Slice DE. 2013 A field comes of age: geometric morphometrics in the 21st century. *Hystrix Ital. J. Mammal* **24**, 7–14.

43. Zelditch ML, Swiderski DL, Sheets AD, Fink WL. 2004 *Geometric morphometrics for biologists. A primer*. Berlin, Germany: Elsevier.

44. Mitteroecker P, Gunz P. 2009 Advances in geometric morphometrics. *Evol. Biol.* **36**, 235–247. (doi:10.1007/s11692-009-9055-x)

45. Seetah K, Cucchi T, Dobney K, Barker G. 2014 A geometric morphometric re-evaluation of the use of dental form to explore differences in horse (*Equus caballus*) populations and its potential zooarchaeological application. *J. Archaeol. Sci.* **41**, 904–910. (doi:10.1016/j.jas.2013.10.022)

46. Polly PD. 2001 On morphological clocks and paleophylogeography: towards a timescale for *Sorex* hybrid zones. In *Microevolution rate, pattern, process* (eds AP Hendry, MT Kinnison), pp. 339–357. Dordrecht, The Netherlands: Springer.

47. Cardini A, Thorington RW, Polly PD. 2007 Evolutionary acceleration in the most endangered mammal of Canada: speciation and divergence in the Vancouver Island marmot (Rodentia, Sciuridae). *J. Evol. Biol.* **20**, 1833–1846. (doi:10.1111/j.1420-9101. 2007.01398.x)

48. Caumul R, Polly PD. 2005 Phylogenetic and environmental components of morphological variation: skull, mandible and molar shape in marmots (*Marmota*, Rodentia). *Evolution* **59**, 2460–2472. (doi:10.1111/j.0014-3820.2005. tb00955.x)

49. Klingenberg CP, Gidaszewski NA. 2010 Testing and quantifying phylogenetic signals and homoplasy in morphometric data. *Syst. Biol.* **59**, 245–261. (doi:10.1093/sysbio/syp106)

50. Rohlf FJ. 2002 Geometric morphometrics and phylogeny. In *Morphology, shape, and phylogeny* (eds N MacLeod, PL Forey), pp. 175–193. London, UK: Taylor & Francis.

51. Adams DC. 2014 A generalized K statistic for estimating phylogenetic signal from shape and other high-dimensional multivariate data. *Syst. Biol.* **63**, 685–697. (doi:10.1093/sysbio/syu030)

52. Mashkour M. 2014 Equids in the northern part of the Iranian central plateau from the Neolithic to Iron Age: new zoogeographic evidence. In *Prehistoric steppe adaptation and the horse* (eds M Levine, C Renfrew, K Boyle), pp. 129–138. Cambridge, UK: University of Cambridge. (McDonald Institute Monographs).

53. Vila E. 2006 Data on equids from late fourth and third millennium sites in Northern Syria. In *Equids in time and space: papers in honor of Vera Eisenman* (ed. M Mashkour), pp. 101–123. Oxford, UK: Oxbow.

54. R Core Team. 2016 *R: a language and environment for statistical computing*. Vienna, Austria: R Foundation for Statistical Computing.

55. Bookstein FL. 1997 Landmark methods for forms without landmarks: morphometrics of group differences in outline shape. *Med. Image Anal.* **1**, 225–243. (doi:10.1016/S1361-8415(97)85012-8)

56. Gunz P, Mitteroecker P. 2013 Semilandmarks: a method for quantifying curves and surfaces. *Hystrix Ital. J. Mammal* **24**, 103–109.

57. Rohlf FJ. 2005 *TpsRelw* 1.4-Thin *Plate Spline Relative Warp*. New York, NY: Ecology & Evolution, State University at Stony Brook.

58. Monteiro LR. 1999 Multivariate regression models and geometric morphometrics: the search for

causal factors in the analysis of shape. *Syst. Biol.* **48**, 192–199. (doi:10.1080/106351599260526)

59. Dray S, Dufour AB. 2007 The ade4 package: implementing the duality diagram for ecologists. *J. Stat. Softw.* **22**, 1–20. (doi:10.18637/jss.v022.i04)

60. Adams DC, Otarola-Castillo E. 2013 geomorph: an R package for the collection and analysis of geometric morphometric shape data. *Methods Ecol. Evol.* **4**, 393–399. (doi:10.1111/2041-210X.12035)

61. Klingenberg CP. 2010 Morphol: an integrated software package for geometric morphometrics. *Mol. Ecol. Resour.* **11**, 353–357. (doi:10.1111/j.1755- 0998.2010.02924.x)

62. Schubert M *et al.* 2014 Characterization of ancient and modern genomes by SNP detection and phylogenomic and metagenomic analysis using PALEOMIX. *Nat. Protocols* **9**, 1056–1082. (doi:10.1038/nprot.2014.063)

63. Maddison WP. 1991 Squared-change parsimony reconstructions of ancestral states for continuous-valued characters on a phylogenetic tree. *Syst. Zool.* **40**, 304–314. (doi:10.2307/ 2992324)

64. Martins EP, Hansen TF. 1997 Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data. *Am. Nat.* **149**, 646–667. (doi:10.1086/286013)

65. Blomberg SP, Garland T, Ives AR. 2003 Testing for phylogenetic signal in comparative data: behavioral traits are more labile. *Evolution* **57**, 717–745. (doi:10.1111/j.0014-3820.2003.tb00285.x)

66. Armitage PL, Chapman H. 1979 Roman mules. *Lond. Archaeol.* **3**, 339–346.

67. Freckleton RP, Harvey PH, Pagel M, Losos AEJB. 2003 Bergmann's rule and body size in mammals. *Am. Nat.* **161**, 821–825. (doi:10.1086/374346)

68. Nova Delgado M, Galbany J, Pérez-Pérez A. 2016 Molar shape variability in platyrhine primates. *J. Hum. Evol.* **99**, 79–92. (doi:10.1016/j.jhevol.2016. 07.006)

69. Gidley JW. 1901 Tooth characters and revision of the North American species of the genus *Equus*. *Bull. AMNH* **14**, 91–142.

70. Dryden IL, Mardia KV. 2016 *Statistical shape analysis: with applications in R*, 2nd edn, 496. Oxford, UK: Wiley.

71. Revell LJ, Harmon LJ, Collar DC. 2008 Phylogenetic signal, evolutionary process, and rate. *Syst. Biol.* **57**, 591–601. (doi:10.1080/10635150802 302427)

72. Jones KE. 2016 New insights on equid locomotor evolution from the lumbar region of fossil horses. *Proc. R. Soc. B* **283**, 20152947. (doi:10.1098/rspb. 2015.2947)

73. Swenson NG, Enquist BJ, Thompson J, Zimmerman JK. 2007 The influence of spatial and size scale on phylogenetic relatedness in tropical forest communities. *Ecology.* **88**, 1770–1780. (doi:10.1890/ 06-1499.1)

74. Losos JB. 2008 Phylogenetic niche conservatism, phylogenetic signal and the relationship between phylogenetic relatedness and ecological similarity among species. *Ecol. Lett.* **11**, 995–1003. (doi:10.1111/j.1461-0248.2008.01229.x)