As a second part of this study, we examine modularity of the equid skull. The concept of modularity [17, 18] has attracted much attention in recent years (e.g. [19–26]), having emerged as a quantitative framework for exploring questions relating to facilitation and constraint in morphological evolution, with the goal of understanding how (and by how much) the direction in which variation is generated is biased [27–29]. Many studies have quantified patterns of modularity in the cranium using inter-trait correlations extracted from geometric morphometric data (see [30] and references therein) and, taken together, their results have supported a common pattern in therian mammals, with some variability in the magnitude of integration among species (e.g. [20, 31]). In contrast, comparatively little is known about the lability of modular patterning and integration magnitudes on relatively short time scales and under selective breeding regimes, although changes in magnitude, rather than patterning, have been implicated as the target for selection [32]. Providing examples of selective breeding for features to suit human activities, the study of domestication events offers an opportunity to empirically examine the role of modularity and integration in the generation of cranial disparity over short evolutionary timescales. A modular structure of the skull is expected to be uncovered for horses, as has been found across a wide spectrum of mammals (e.g. [31]), and we assess the fit of our shape data to functional and developmental hypotheses for modular patterning [33] that have been previously tested in the mammalian cranium. According to Darwin’s hypothesis, domesticated horses should show more variation in shape than the wild equid species. If this hypothesis is supported, then we should find differences in integration and disparity measures for cranial modules between wild equids and domesticated horses. To do so, we assess whether cranial modules display a) higher or lower magnitudes of integration and b) high or low disparity for domesticated horses and wild equids, and c) we investigate whether there is a relationship between module integration/disparity and regions of the cranium showing most variability in shape among domesticated horses.

Our aim is to characterize and quantify the patterning and magnitude of shape variation in the skulls and teeth of domesticated horses compared to the wild species of Equus. We use geometric morphometric methods to: a) test Darwin’s hypothesis that the magnitude of intraspecific disparity in horses is larger than the interspecific disparity in equids, b) examine the extent to which domestication influenced tooth shape in equids, c) investigate whether the patterning of shape variation in horse skulls reflects a modular structure, specifically identifying the model best supported for the landmark data by evaluating four modular hypotheses that reflect developmental and functional trait interactions in the cranium, and d) quantify differences in the magnitude of modularity and integration between domesticated horses and wild equids.

Given the well-documented palaeontological record of horses [34], these animals offer the possibility to compare diversification in macroevolutionary and microevolutionary scale like few others. In fact, previous classic studies by Radinsky [35, 36] investigated some of the cranial transformations with morphometric approaches typical of that time. Our study expands the studies of the extant species using newer morphometric approaches and provides the basis for future works comparing also the fossil record of the clade.

The results of our study on skull shape variation support Darwin’s hypothesis that the intraspecific disparity in horses is larger than the interspecific disparity in equids. The shape variation of domesticated horses is not only larger than that of the closest relative (Przewalski’s horses) but similar to the shape variation of all the wild equid species existing today. Horses do not only dominate the Procrustes shape space when comparing crania, but also comparing mandibles or teeth - showing higher shape variation in all tested elements.

The overall classification among domesticated horses (H), Przewalski’s horses (P), donkeys (D), and zebras (Z) had an average accuracy of 44.5% (range: 34.5% - 78.4%), which lies significantly above the random chance accuracy of 25%. When we excluded the Przewalski’s specimens due to their small sample size the average accuracy increased to 89.2% (range: 63.9% - 98.2%). The separation of caballine from non-caballine taxa and the clustering within these sister clades of H & P, and D & Z by the Mahalanobis distances, is in accordance with the phylogenetic relationship of equids. All four groups descend from a common ancestor around 4–4.5 myr BP, with zebras and donkeys splitting around 2.8 myr BP and Przewalski’s horses and the wild ancestor of today’s domesticated horses splitting 38–72 kyr BP [2, 3, 37]. The first PC of the cranial shape data tended to separate caballine from non-caballine taxa. The accompanying shape differences are dominated by an elongation of the occipital part of the cranium in zebras and donkeys, which previously have been shown to be a distinct character to separate these two taxa from the domesticated horses (for a more detailed morphological description of the skulls of each group see Table 4, Fig. 3, [38, 39]).

Type	Domesticated equids				Wild equids
	Draft horses (n = 37)	Ponies including Falabella (n = 40)	Medium horses (n = 38)	Light horses (n = 14)	Przewalski’s horses (n = 5)	Donkeys (n = 33)	Zebras (n = 47)
Trait
Average skull length based on studied sample (±SD)	56 (±2.5) cm	43 (±3.6) cm	52 (±4.5) cm	51 (±2.3) cm	46 (±2.1) cm	41 (±4.6) cm	48 (±4.1) cm
General skull shape	Skull breadth is half of skull length	Skull breadth is half of skull length	Elongated	Elongated	Skull breadth is half of skull length	Skull breadth is almost skull length	Very elongated in the Grevy Zebra
Frontals	Broad	Broad in some breeds	Broad in some breeds	Narrow	Narrow	Narrow	Narrow
Degree of cranial flexion (after Bennett, 1980)	High	High in smaller breeds	Medium	High	Medium	High	Low
Nasals	Convex; End at ventral end of premaxillary-maxillary suture	Concave in smallest breed, straight in larger breeds; End before ventral end of premaxillary-maxillary suture	Straight with exceptions in some breeds; End at ventral end of premaxillary-maxillary suture	Concave; End at ventral end of premaxillary-maxillary suture	Straight or slightly concave; End at ventral end of premaxillary-maxillary suture	Straight or slightly concave; End at ventral end of premaxillary-maxillary suture	Straight or slightly concave; End at ventral end of premaxillary-maxillary suture
Angle between premaxilla and maxilla	Obtuse	Obtuse, almost straight	Straight	Straight	Right to obtuse	Obtuse to straight	Obtuse to straight
Facial crest	Prominent	Prominent	Prominent	Less prominent	Prominent	Prominent	Prominent
Zygomatic, temporal and postorbital bars	Very broad, Horizontal to skull midline	Broad to normal; Horizontal to skull midline	Thin, with the temporal being broader; Horizontal to skull midline	Very thin; Horizontal to skull midline	Normal; Angled to skull midline	Normal; Angled to skull midline	Very broad; Angled to skull midline
Orbit shape	Roundish with a more squared part at the dorso-posterior edge	Almost round	Egg-shaped	Oval	Oval	Oval	Round to egg-shaped
Zygomatic process of the temporal bone	Flat	Flat	Curved	Curved	Flat	Flat	Flat
Length of the lambdoidal crest	Does not overlap the condyles in lateral view	Does not overlap the condyles in lateral view	Does not overlap the condyles in lateral view	Does not overlap the condyles in lateral view	Overlaps the condyles in lateral view	Overlaps the condyles in lateral view	Overlaps the condyles in lateral view
Cranial vault	Flat	Rounded in smaller breeds	Flat	Flat	Flat	Rounded	Flat

In contrast to the results of the cranial analysis, we did not find clear group-distinguishing shape differences in the mandible or in any of the investigated teeth. We found the highest Procrustes variance in the mandible data set, very closely followed by the cranium. All four teeth showed a lower Procrustes variance pointing towards less variability, a result also reflected in the low magnitude of disparity for landmarks belonging to the molar (MR) module, probably related to dietary constraints. Our findings on tooth shape differences are congruent with the “long-fuse” hypothesis on teeth by Seetah et al. [14], stating that “shape changes in equids have been modest […] until the development of modern breeds in recent centuries”. Cucchi et al. [16], however, found a strong taxonomic pattern in the shape of the enamel folding, allowing for a more distinct taxonomic separation at the species level. These resulting differences are most likely due to the different choice of teeth ([14]: UP4; [16]: LP3 – LM2), as suggested by one of the articles [16], different wear stages, and/or the different methods used (landmarks vs. semilandmarks).

Modular patterning of the cranium is well-supported by empirical evidence, representing shared development and functional associations between parts of the cranium, resulting in their parcellation into semi-autonomous units. Consistent with patterns recovered for other placental mammals [20], our results indicate that trait variation in the equid cranium is best supported by a six-module hypothesis. These six modules reflect functional groups: the anterior oral-nasal (AON) and molar (MR) modules represents the primary masticatory apparatus, the zygomatic-pterygoid (ZP) module includes the region of attachment for masticatory muscles, the orbit (ORB) module contains the visual sensory organs, the cranial vault (CV) supports and protects the brain, and the cranial base (CB) supports the braincase and is the point of attachment between the skull and axial skeleton. Previous analyses of module disparity and integration for these six modules in a sample of carnivorans provided some support for highly integrated modules showing low disparity, particularly the basicranium (CB), and weakly integrated modules showing high disparity (ORB and ZP) [20, 28]. Our results are broadly consistent with this trend (Fig. 2), which is suggestive of strong integration acting to limit trait variation or the direction of response to selection. The ZP module in horses is also found to be weakly integrated and showing high disparity. Among the carnivoran sample, the molar region showed an unusual pattern of high disparity and high levels of integration [28], in our sample the molar module is also recovered as the most highly integrated module but displays the lowest levels of disparity for both wild and domesticated forms. The discrepancy between these results is likely explained by the diversity of dietary habits represented by the carnivoran sample (e.g. hypercarnivores, insectivores, frugivores) in that earlier study, which resulted in a high disparity among the landmarks captured in the molar module.

Selection acting on shared developmental and functional processes can result in an uncoupling of trait associations at different levels [40–43], providing evidence for complex interactions between modularity and selection. Following, it might therefore be expected that domestication events, as examples of selective breeding regimes, could alter patterns of modularity and integration, and these alterations may differ among breeds, acknowledging that the features targeted for selection (e.g. gait, conformation) are likely to differ for some breeds. There exist few empirical tests of this hypothesis and results on dogs are inconclusive, with reports that patterns of integration have remained stable despite the morphological diversification associated with domestication [9, 44], but also that high module disparity is associated with greater cranial shape variation in dogs compared to wolves [45]. In contrast to Parr et al. [45] our results show highly similar magnitudes of module disparity among wild and domesticated forms, and instead lower magnitudes of module integration are recovered in domesticated horses compared to their wild relatives. Variability in integration magnitude, as recovered here, rather than patterning has been proposed to underlie cranial diversity in mammals [31, 32], such that general conservatism in patterning across mammals may be explained as a product of stabilizing selection on basic developmental processes whereas directional selection could act by altering magnitudes of integration. A recent study conducted simulations to test the role of integration in generating morphological disparity and noted that integration may not affect disparity in morphospace in the way that it is usually measured (as a volume of occupied morphospace or as a measure of dissimilarity), making the relationship between morphospace occupation and modularity results potentially difficult to interpret [46]. That study did not compare shape variation and its partitioning into modules, however our PCA plots indicate that the main axes of shape variation in the equid sample are spread across landmarks located in at least three modules (CV, AON and ZP). Similarly, Parr et al. [45] found shape variance in wild and domesticated dogs to be spread across modules with different magnitudes of integration. It has been suggested that modularity may repartition variance along new directions in morphospace, thereby exploring a greater volume, however the so-far limited empirical evidence appears to raise the question of the extent to which those new directions may be aligned with the axes recovered by eigen-decomposition of shape variables into mathematically orthogonal axes, as happens in ordination techniques such as PCA.

A total of 216 crania, 224 mandibles, 198 upper and 141 lower second premolars (U2P and L2P respectively), and 192 upper and 144 lower third molars (U3M and L3M respectively) were analysed (Table 2).

We examined specimens from the following collections: Museum für Naturkunde Berlin (MfN Berlin, Germany), Institut für Haustierkunde (Christian-Albrechts-Universität of Kiel, Germany), Museum für Haustierkunde "Julius Kühn" (University of Halle, Germany), Naturhistorisches Museum Wien (NHW Vienna, Austria), and Museo de la Plata (MLP La Plata, Argentina). The dataset includes all recent species of the genus Equus [37]. Due to inconsistent species assignment within zebras and donkeys across museums, we analysed all zebra (cranium n = 47; mandible n = 48) and donkey (cranium n = 31; mandible n = 33) species as one group, respectively. We further included five crania and two mandibles of Przewalski’s horses. The largest number of specimens in our data set belongs to the domesticated horses (cranium n = 133; mandible n = 141) including the following breeds: Ancient Breed (Roman period), Anglo-Norman, Arab, Birkenfelder, Belgian Draft, Bosnian Pony, Clydesdale, Exmoor Pony, Falabella, Galician Farm Horse, Grisons (Graubündner), German Riding Pony, Hannoverian, Hackney, Holstein, Hungarian, Huzule, Icelandic Horse, Indian Pony, Kladrubian, Konik, Kosarian, Lipizzan, Mongolian, Norik, Oldenburgian, Pinzgau, Polish Farm Horse, Scottish Pony, Seneca Sarajevo, Shetland Pony, Shire, Styrian, Suffolk, English Thoroughbred, Togo Pony, Trakehner, and Welsh (Table 2).

Analyses of cranial, mandibular and teeth size and shape were performed using landmark-based geometric morphometric (GMM) approaches. The crania and mandibles were measured in three-dimension (3D) using a MicroScribe ® MLX6 (Revware, Inc., Raleigh, North Carolina, USA), while the teeth were measured in two-dimension (2D) using high resolution photographs. A total of 62 type I and type II landmarks [47] were collected on the cranium (Table 5, Additional file 3). The dorsal and ventral sides of the crania (Fig. 4) were measured separately and were subsequently combined using three reference landmarks (numbered 1, 2, and 33, Table 5). For the mandible 24 type II landmark coordinates were measured (Table 5, Fig. 4, Additional file 3).

	Position	Type
Cranium
1–2	Posterior tip of the upper third incisor	II
3–4	Posterior most point of the nasal-premaxilla suture	I
5–6	Premaxillary-maxillary-nasal suture	I
7–8	Dorsoposterior tip of the infraorbital foramen	II
9–10	Anterior tip of the facial crest	II
11	Nasion, nasal-frontal suture, midline	I
12–13	Junction of the lacrimal, maxilla, and nasal sutures	I
14–15	Zygo-lacrimal suture on the orbital margin	I
16–17	Lacrimal-frontal suture on the orbital margin	I
18–19	Supraorbital foramen	II
20–21	Anterior tip of the zygo-temporal suture	I
22–23	Posterior tip of the zygo-temporal suture	I
24–25	Dorsal tip of the frontal-temporal suture	I
26–27	Ventroposterior tip of the zygomatic process	II
28–29	Dorsalmost point of the vertically orientated posterior margin of the zygomatic process	II
30–31	Ventrolateralmost point of squamous part of temporal bone	II
32	Anterior tip of the occipital triangle	I
33	Posterior tip of the nuchal crest	II
34–35	Dorsolateral tip of the nuchal crest	II
36	Dorsalmost point on the margin of the foramen magnum	II
37	Point between first incisors from ventral side	II
38–39	Posteriormost tip of the premaxillary-maxillary suture, ventral	I
40–42	Anterior tip of the second premolar	II
41–43	Posterior tip of the third molar	II
44	Posteriormost point of the incisive canal	II
45	Posterior tip of the palatine process of the incisive bone	I
46	Posterior tip of the palatine-palatine suture	I
47–48	Distal tip of the pterygoid hamulus	II
49–51	Anterior tip of the caudal alar foramen	II
50	Posterior tip of the vomer on the midline	II
52–53	Medial tip of the mandibular fossa	II
54–55	Canal for hypoglossal nerve	II
56–57	Fossa medial of the paracondylar process	II
58–59	Distal tip of the paracondylar process	II
60	Ventral tip of the foramen magnum	II
61–62	Posteriormost tip of the occipital condyle	II
Mandible
U1	Posterior point between first incisors	II
U2–U3	Posterior tip of the third lower incisor	II
U4–U5	Posterior tip of the canine	II
U6	Posterior tip of the mandible on the midline	II
U7–U8	Anterior tip of the second premolar	II
U9–U10	Posterior tip of the third molar	II
U11–U12	Junction of the bases of the coronoid and condylar processes	II
U13–U14	Lateral tip of the condylar process	II
U15–U16	Medial tip of the condylar process	II
U17–U18	Posterior tip of the mandibular mental foramen	II
U19–U20	Vascular notch of the mandible	II
U21–U22	Maximum curvature of the angle of the mandible right behind the vascular notch	II
U23–U24	Maximum curvature of the angle of the mandible	II
Cranium 2P
1	Maximum curvature of the metastyle	II
2	Maximum curvature of the mesostyle, distal side	II
3	Maximum curvature of the mesostyle, medial side	II
4	Maximum curvature of the anterior accessory rib	II
5	Maximum curvature of the parastyle	II
6	Maximum curvature of the protocone, mesial/labial side	II
7	Maximum curvature of the protocone, buccal side	II
8	Maximum curvature of the protocone, distal/labial side	II
9	Maximum curvature of the hypocone	II
Cranium 3 M
1	Maximum curvature of the metastyle	II
2	Maximum curvature of the mesostyle	II
3	Maximum curvature of the parastyle	II
4	Maximum curvature of the preprotoconal groove	II
5	Maximum curvature of the protocone, mesial side	II
6	Maximum curvature of the protocone, distal side	II
7	Maximum curvature of the post protoconal valley, buccal side	II
8	Maximum curvature of the post protoconal valley, labial side	II
9	Maximum curvature of the hypocone, buccal side	II
10	Maximum curvature of the hypocone, labial side	II
Mandible 2P
1	Maximum curvature of the protoconid, anterior/ lingual side	II
2	Maximum curvature of the preflexid, anterior side	II
3	Maximum curvature of the preflexid, posterior side	II
4	Maximum curvature of the metaconid, anterior side	II
5	Maximum curvature of the metastylid, posterior side	II
6	Maximum curvature of the postflexid, anterior side	II
7	Maximum curvature of the postflexid, posterior side	II
8	Maximum curvature of the hypoconulid, posterior side	II
9	Maximum curvature of the hypoconid, anterior/ buccal side	II
10	Maximum curvature of the ectoflexid, lingual side	II
11	Maximum curvature of the protoconid, posterior/ buccal side	II
12	Maximum curvature of the protoconid, anterior/ buccal side	II
Mandible 3 M
1	Maximum curvature of the hypoconulid, posterior side	II
2	Maximum curvature of the entoconid, anterior side	II
3	Maximum curvature of the entoflexid, posterior side	II
4	Maximum curvature of the entoflexid, anterior side	II
5	Maximum curvature of the metastylid, posterior side	II
6	Maximum curvature of the metaconid, anterior side	II
7	Maximum curvature of the metaflexid, posterior side	II
8	Maximum curvature of the protoconid, anterior side	II
9	Maximum curvature of the protoconid, posterior side	II
10	Maximum curvature of the ectoflexid, lingual side	II
11	Maximum curvature of the hypoconid, anterior/ buccal side	II
12	Maximum curvature between the hypoconid and hypoconulid, lingual side	II

Phenotypic variation of the four studied teeth was assessed using 9 to 12 two-dimensional landmarks (Type II) (Table 5, Fig. 4, Additional file 4, 5, 6 and 7) following Seetah et al. [48] for the upper teeth that was adapted for the lower teeth. The landmark coordinates were collected on high resolution photographs using TPSDig2 [49]. The photographs were all taken in a standardized manner using a Canon Eos 600d with a Canon EF 24-105 mm f/4 L S USM lens from lateral and dorsal view with a scale bar for size reference.

We described and compared shape changes in various skeletal features among extant equid species using geometric morphometrics. Our results support Darwin`s hypothesis that shape variation in the skull of domesticated horses is similar to the shape variation of all wild extant equid species. Our study further shows that lower magnitudes of integration among six cranial modules are found in domesticated horses compared to their wild relatives. Future research could address the relation between integration and disparity, investigating the relation between the two during the domestication process of diverse species.