Skip to main content

Cranial shape diversification in horses: variation and covariation patterns under the impact of artificial selection


The potential of artificial selection to dramatically impact phenotypic diversity is well known. Large-scale morphological changes in domestic species, emerging over short timescales, offer an accelerated perspective on evolutionary processes. The domestic horse (Equus caballus) provides a striking example of rapid evolution, with major changes in morphology and size likely stemming from artificial selection. However, the microevolutionary mechanisms allowing to generate this variation in a short time interval remain little known. Here, we use 3D geometric morphometrics to quantify skull morphological diversity in the horse, and investigate modularity and integration patterns to understand how morphological associations contribute to cranial evolvability in this taxon. We find that changes in the magnitude of cranial integration contribute to the diversification of the skull morphology in horse breeds. Our results demonstrate that a conserved pattern of modularity does not constrain large-scale morphological variations in horses and that artificial selection has impacted mechanisms underlying phenotypic diversity to facilitate rapid shape changes. More broadly, this study demonstrates that studying microevolutionary processes in domestic species produces important insights into extant phenotypic diversity.

Peer Review reports


The phenotypic diversification of domestic species provides a unique and accelerated perspective on evolutionary processes. Artificial selection has proven able to strongly impact the phenotype of domestic taxa over short time frames, producing great amount of morphological disparity often exceeding that of wild counterparts [1,2,3,4,5,6]. Indeed, sustained selection by breeders (e.g. for specific morphological, functional or behavioral features) can generate novel shape variation and contribute to large-scale phenotypic diversification in a few generations [3]. Among domestic taxa, the morphological diversification in domestic horses (Equus caballus) appears as particularly suitable for investigating rapid evolutionary processes having produced substantial shape variation in a few short centuries [7]. Indeed, in terms of both breeding practices and genomic makeup, domestic livestock such as extant horse breeds largely has its origins in the eighteenth century [8,9,10].

This ability of artificial selection to strongly impact the morphological features of domestic animals raises the issue of the existence of microevolutionary mechanisms facilitating rapid shape changes [11]. Phenotypic diversification is underpinned by several mechanisms which determine the variation that is available for selection to act upon. Notably, the developmental and functional relationships between the different component parts of organisms are known to influence patterns of morphological variation [12,13,14]. This tendency of morphological traits to covary, or “morphological integration”, is thus a key factor influencing morphological diversification under selection [15,16,17,18,19]. A set of highly correlated morphological traits, acting in a semi-autonomous way, is referred as a module [20]. Morphological modularity and integration are tightly related to evolution as they are thought to influence the “evolvability” (i.e. capacity to evolve) [19] of morphological traits, by constraining the variation of individual traits or facilitating evolution through coordinated shape changes [21,22,23]. Selective processes may cause changes in modularity patterns or in magnitude of integration, and therefore these can be examined as a way to understand how interactions among traits drive or limit the generation of variation in evolution [15, 19].

Cranial structures are commonly used as a model for studying morphological modularity as they are functionally and developmentally well known, providing thus hypotheses of modular patterning [15, 24, 25]. Previous studies demonstrated the conservation of cranial modularity patterns (i.e. relationships between traits) across placental mammals [24, 26, 27], including equids [4]. Conversely, the magnitude of morphological integration (i.e. intensity of the association between traits) has been shown to vary considerably across taxa, which could have consequences on evolvability [26, 27] and thus facilitate cranial diversification in domestic taxa [1, 3, 6, 28]. Whether modularity constrains or facilitates evolution is still a subject of debate, and no clear relationship between degree of modularity and shape disparity has yet been demonstrated [1, 22, 28]. In the genus Equus, a previous study demonstrated a lower magnitude of morphological integration in domestic than in wild taxa. This suggests that artificial selection would be associated with reduced inter-trait relationships, thus potentially contributing to increase flexibility and enhance shape diversification [4]. This could also suggest that variable intensities of alteration in the magnitude of integration could potentially be observed across horses, according to the degree to which they have been submitted to artificial selection.

The wild ancestor of domestic horses no longer exists, and the last surviving population of wild horses is the Przewalski’s horses (Equus przewalskii) [29]. It constitutes a distinct species, only other representative of the caballine lineage. Although Przewalski’s horse went “extinct in the wild” in the 1960s [30], they have survived in captivity [31, 32] and, since the late twentieth century, have been progressively reintroduced into the wild [33]. As a different species whose morphology has also likely been impacted by modern captivity, inbreeding depression [29, 34] and potentially human management in the past [35], Przewalski’s horses cannot provide a direct analogue for pre-domesticated horses in studying domestication processes. However, they do nonetheless provide a population of closely-related horses that is not subjected to artificial selection.

The influence of selective breeding may also be assessed within populations of domestic horses, as different evolutionary pathways may explain the formation of different horse breeds. Some current standardized breeds have been subject to high levels of artificial selection and have been forged by reproductive isolation imposed by breeders [36]. This is true of most racehorse breeds, whose breeding is aimed at particular morphological features or athletic performance, and of draft horses, on which considerable selective pressures were exerted, mainly on overall size or body mass [37]. In contrast to breeds formed through deliberate human choice for specific features (e.g. conformation, performance), other breeds might be better characterized as “landraces” (deriving their shared genetic and morphological traits from natural conditions due to long isolation within a specific environment and having been mostly shaped without deliberated breeding decisions) [36]. These horses are generally free-ranging and are breed under conditions of minimal human intervention [38, 39]. Finally, feralization, which is the process by which domestic animals return to the wild, constitutes a third kind of evolutionary pathway. In this case, phenotypic traits of the rewilded animals may have been impacted by natural selection despite their ancestral state of domestication [40,41,42,43].

In the present study, we contribute to assess the impact of artificial selection on the cranial morphology of domestic horses using 3D geometric morphometrics. We explore variation in extant groups to determine whether morphological differences between horse populations reflect divergent evolutionary mechanisms implicated in their formation. We also investigate the impact of artificial selection on modularity and integration patterns, to gain insight into underpinning evolutionary mechanisms having allowed the dramatic shape diversification in domestic horses. Hypothesizing that the varying degree of artificial selection to which they have been subject could have differently impacted their morphological traits, we compare the shape variation and covariation patterns among different groups or breeds, known to having been submitted to varying degrees of artificial selection: highly standardized breeds (i.e. race- and draft horses), landraces (i.e. Mongolian, Icelandic, Shetland, Pottok), domestic breeds returned to the wild since several generations (i.e. American feral horses) and last surviving species of wild horses, supposed to have not been subject to artificial selection (i.e. Przewalski’s horses). We explore the shape variation to better understand how artificial selection would have impacted morphological diversification in horses. We then investigate potential differences among these groups in shape covariation. The aim here is to detect potential changes in the structure of modularity, in magnitude or in patterns of integration, and to try to relate them to the varying intensities of artificial selection to which the groups have been submitted.

Material and methods


Our analyzed dataset includes a total of 91 skulls from both domestic (Equus caballus, n = 74) and Przewalski’s (Equus przewalskii, n = 17) horses housed in the collections of several institutions (see Additional file 1). The domestic horses include 21 breeds or landraces, selected to be representative of a large range of diversity in morphology and size: draft (n = 20) and racehorses (n = 21) of various breeds, Mongolian (n = 15), Icelandic (n = 3), Shetland (n = 4), Pottok (n = 3) and American feral horses (n = 8). Due to the small sample size linked to individual breeds, the racing and draft breeds, respectively, were grouped together in most analyses, according the classification of the International Federation for Equestrian Sports, on the basis of functional and genetic criteria. The total sample consists of both males and females. Only adult specimens with all permanent teeth were used (older than 4 years; see Additional file 1).

Acquisition of data

Skulls were digitized in three dimensions using several devices (an Artec Space Spider for n = 21, Artec Eva for n = 22, NextEngine laser scanner for n = 22 and photogrammetry for n = 26; see Additional file 1). Bone shape was quantified using a set of anatomical landmarks and sliding semilandmarks on curves and surfaces. We defined a total of 1482 landmarks, including: 69 anatomical landmarks (from Hanot et al. [44]), 162 sliding semilandmarks placed on 20 curves and constrained by anatomical landmarks, and 1250 surface sliding semilandmarks (Fig. 1). Anatomical landmarks and curves were placed on the three-dimensional bone models using the software IDAV Landmark v. 3.0 [45]. We manually digitized surface sliding semilandmarks on a template and then semi-automatically projected these onto each mesh via Thin-Plate Spline (TPS) deformation using the “placePatch” function from the R package Morpho [46]. Semilandmarks on curves and surfaces were slid along their tangent vectors/planes to minimize bending energy using the “slider3d” function from Morpho package [46]. Symmetrization of the landmark coordinates along the median plan was performed using the “symmetrize” function from Morpho package [46].

Fig. 1.
figure 1

3D view of horse skull showing the location of the 69 anatomical landmarks (in red), 162 sliding semilandmarks placed on curves (in blue) and 1250 surface sliding semilandmarks (in green). See Additional file 2 for landmark definition

Shape analyses

A generalized Procrustes Analysis (GPA) was implemented on the landmark data to remove the effects of location, scale, and orientation of the configurations [47]. We then performed a Principal Component Analysis (PCA) on the Procrustes residuals to reduce the dimensionality of the multivariate datasets [48,49,50] producing new independent variables (Principal Components, PC) maximizing the variance within the data. The distribution of the data in shape space was displayed by plotting the two first PCs. Visualizations of the shapes associated with extreme parts of the PCs were produced using a TPS deformation of the consensus surface.

We tested differences in shape and size between groups using respectively a Multivariate Analysis of Variance (MANOVA), on PC accounting for more than 95% of the shape variability, and an Analysis of Variance (ANOVA), with Benjamini–Hochberg correction for multiple comparisons [51]. The effect of allometry was assessed by regressing shape against the log10-transformed centroid size (CS). To determine if the different groups have a common allometry, we performed a Procrustes ANOVA to test the homogeneity of allometric slopes. Allometry-free shapes were extracted from the residuals of the multivariate regression models [52]. The analyses below were then performed on both normal and allometry-free shape variables.

We assessed the impact of sexual dimorphism on shape by testing the difference between males and females using a MANOVA on shape data accounting for more than 95% of the shape variability. A two-way MANOVA was also used to assess the interaction between sexual dimorphism and difference between groups (of breeds). Due to the small number of male specimens with known castration status in our study sample, we did not test the potential impact of gelding on shape variation.

We performed Canonical Variate Analyses (CVA) on the first PCs (accounting for more than 95% of the shape variability) to describe the differences among groups. Visualizations of the shapes associated with extreme parts of the CVs were produced using a TPS deformation of the consensus surface from a projection of the CVs into the original space.

We computed the Procrustes variance to assess morphological disparity [53] within each of the main groups (racehorses, Mongolian, draft and Przewalski’s horses). To make disparity values comparable across modules, values were scaled by being divided by the number of landmarks and semilandmarks included in each module [54]. Pairwise comparisons between groups were also performed with Benjamini–Hochberg correction for multiple comparisons. We carried out these analyses using the “morphol.disparity” function from the Geomorph package [55].

Integration and modularity

For a better grasp of the processes underlying the generation of variation, we explored shape covariation patterns. The structure of modularity within the skull and its conservation across horses were examined, as well as that of patterns of inter-module morphological integration. Magnitude of between and within-module integration was also assessed, and then, related to module variance to evaluate how integration influences evolvability. Comparisons between the different groups were finally conducted to question the impact of artificial on modularity and integration.

Modular patterning

We examined the general structure of modularity in the horse skull using hypotheses of functional and developmental influences. Five alternative partitions of landmarks into modules were defined (Fig. 2) according to models previously proposed for skull modularity: (1) absence of distinct modules; (2) ossification model (dermal/endochondral); (3) tissue origin model (neural crest/paraxial mesoderm); (4) mammalian model (anterior oral-nasal/orbital/molar/basicranium/zygomatic-pterygoid/cranial vault; Goswami [24]); (5) functional model (oral/nasal/orbital/masticatory/basicranium/vault; Cheverud [15]). Hypotheses for modularity were investigated with two different approaches: using the EMMLi (“Evaluating Modularity with Maximum Likelihood”) analysis and the Covariance Ratio (CR) [56]. EMMLi is a maximum likelihood approach, implemented in the EMMLi package [57], which allows to compare different models of modularity. Because it is not exhaustive in its comparison of models and has been demonstrated as favoring the most-parameterized models [54, 58], we used it coupled with the CR to see if both methods support similar models of modularity. The CR uses the pairwise covariances between variables to quantify modular structure, with modularity in the structure considered as significant when the CR is small relative to the distribution of values obtained under the null hypothesis of random associations of variables. To compute it, we used the “modularity.test” function from Geomorph library [55]. This procedure was carried out on each group separately to confirm the hypothesized stability in modular patterning across therian mammals [24]. Because the use of surface semilandmarks may impact modularity patterns (by exaggerating within-module correlations and thus increasing global modularity) [54, 59,60,61], the procedure was also computed on landmarks and curves only, for matters of comparison. The supported model was then used in further analyses of integration and modularity.

Fig. 2
figure 2

Alternative partitions of the horse skull with model 1: no module; model 2: ossification (DER dermal, END endochondral); model 3: tissue origin (NC neural crest, PM paraxial mesoderm); model 4: mammalian (AON anterior oral-nasal, ORB orbital, MOL molar, CB basicranium, ZP zygomatic-pterygoid, CV cranial vault); model 5: functional (ORA oral, NAS nasal, OB orbital, MAS masticatory, BAS basicranium, VAU vault)

For each module, we performed separate Procrustes superimpositions to examine their shape variation irrespective of their position in the skull [47, 62]. A CVA was conducted on shape data and visualizations of the shape changes were produced using a TPS deformation of the consensus surface. We assessed morphological disparity computing the Procrustes variance [53] with pairwise comparisons between groups.

Covariation patterns

To investigate covariation between the different modules, we performed Partial Least Squares (PLS) analyses on the adjacent modules. The two block-PLS (2B-PLS) extracts the principal axes of covariation from a covariance matrix of the two shape datasets [63, 64]. The first PLS axes were plotted and associated shape deformations visualized using a TPS deformation of the consensus surface. We performed 2B-PLS using the “two.b.pls” function from the Geomorph package [55]. To compare covariation patterns among groups, we computed the angular difference of the first PLS axis for each group and tested the null hypothesis that the axes are no more similar than vectors having random directions using MorphoJ [65].

Magnitude of morphological integration

The degree of cranial modularity was assessed using the CR [56]. The CR value and associated effect size (Z-score which provides a standardized measure) [58] describe the degree of modularity within the structure, with low values corresponding to a high degree of modularity. Being insensitive to sample size or number of variables, this measure is well adapted to small sample groups. To compare the effect sizes among groups, we performed two-sample z-tests (with Benjamini–Hochberg correction for multiple comparisons) using the “compare.CR” function from Geomorph library [55].

The eigenvalue dispersion of the covariance matrices, computed from standard deviation, was used as a measure of integration within each module [66, 67]. To allow comparison between different matrices, their respective size was taken into account by dividing the observed standard deviation of eigenvalues by the theoretical maximum, producing a relative standard deviation [66, 68]. To obtain a result invariant to sample size, the expected value for the particular sample size (with no integration) was computed to calculate the deviation from the expected value [67]. The range of the relative eigenvalue variance is from zero to one, a value of zero corresponding to an absence of integration (i.e. all eigenvalues being equal). These integration indexes were computed using the “CalcEigenVar” function from the evolqg library [69].

Finally, the magnitude of morphological integration between adjacent modules was assessed by 2B-PLS using the “integration.test” function from the Geomorph library [55].

For statistical matters, the magnitude of integration was assessed only in the main groups (racehorses, Mongolian, draft and Przewalski’s horses). Moreover, the procedure was computed on landmarks and curves only, for matters of comparison. For all the analyses previously described, we considered test results as significant when p-values (p) were below 0.05. All the plots were performed using the ggplot2 library [70].


Size and shape variation

MANOVAs on shape data reveal significant pairwise differences among all the main groups (i.e. racehorses, Mongolian horses, draft horses and Przewalski’s horses; p < 0.05) but no significant difference between males and females (p > 0.05). The ANOVA on CS indicates significant pairwise differences between several breed groups (Table 1) with draft horses displaying higher bone CS than all the other groups, racehorses displaying higher bone CS than most of the groups and Shetland horses displaying lower bone CS than all the other groups (Fig. 3).

Table 1 Pairwise comparisons of centroid size between groups (significant differences with *; p < 0.05)
Fig. 3
figure 3

Boxplots of the variation in log-transformed centroid size of the skull for the different groups

The multivariate regression of the shape against the log10-transformed CS indicates significant and strong impact of size on shape (R2 = 0.93). This strong influence of allometry on shape variation is reflected by the distribution of the specimens according to size along PC1 (which accounts for 25.7% of the total shape variation; Fig. 4). Anatomically, smaller skulls are broader and characterized by a rounder braincase, a concave nasal bone and larger orbits comparatively to total skull size.

Fig. 4
figure 4

Scatter plot of the two first PCs of the PCA performed on the cranial shape data (crosses represent the centroid value for each group) with visualization of the shape changes along the axes (a). Anatomical location and intensity of the shape deformation associated with PC1 (b) and PC2 (c) using distances from the shapes at the negative to positive extreme of the axis (from blue to red)

The result of the Procrustes ANOVA shows that the different breeds share a common allometry (p > 0.05). A multivariate regression on the whole sample was thus performed to obtain allometry-free shapes.

Allometry-free shape variation

The result of the MANOVA on allometry-free shapes indicates pairwise differences among all the main groups (i.e. racehorses, Mongolian horses, draft horses and Przewalski’s horses; p < 0.05). On allometry-free shapes, differences between males and females are significant (p < 0.05; see Additional file 3). However, the two-way MANOVA demonstrates the absence of interaction between breed groups and sexual differences, thus allowing to consider that sexual dimorphism does not bias our between-group results.

The distribution of the specimens along the first axes of the PCA on allometry-free shapes does not reveal clear differentiation between individual race and draft horse breed groups, supporting their bundling into larger groups (Fig. 5). Similarly, we observe an important overlap between the main groups. Globally, the two first axes of the PCA express slight differentiation between racehorses and the three other main breed groups (i.e. Przewalski’s, Mongolian and draft horses). Other landraces (i.e. Icelandic, Pottok and Shetland horses) and to a lesser extent, American feral horses, tend to occupy an intermediate position.

Fig. 5
figure 5

Scatter plot of the two first PCs of the PCA performed on the allometry-free cranial shape data (crosses represent the centroid value for each group) with visualization of the shape changes along the axes (a). Anatomical location and intensity of the shape deformation associated with PC1 (b) and PC2 (c) using distances from the shapes at the negative to positive extreme of the axis (from blue to red)

To simplify the description of differences among breed groups, we performed a CVA on allometry-free shape data. The distribution of specimens along the first CV (55.6%) shows that shape differentiation is mainly driven by the species difference, as the axis clearly distinguishes Przewalski’s horses from domestic ones (Fig. 6). Anatomical changes along CV1 mainly concern the occipital region of the skull: condyles are more developed posteriorly in domestic horses, and exhibit a less extended nuchal crest. This analysis also reveals a difference in the general width of the skull, with Przewalski’s horses displaying a broader skull from incisive to orbital areas. Finally, Przewalski’s horses exhibit a slightly straighter nasal bone and a more robust incisive area.

Fig. 6
figure 6

Scatter plot of the two first CVs of the CVA performed on the allometry-free cranial shape data (42 PCs) with visualization of the shape changes along the axes. a Anatomical location and intensity of the shape deformation associated with CV1 (b) and CV2 (c) using distances from the shapes at the negative to positive extreme of the axis (from blue to red)

Domestic horses are differentiated along the second CV (30.6%), with racehorses occupying the positive side of the axis. Shetland ponies, Pottok and Icelandic horses are in an intermediate position, followed by draft and American feral horses. Finally, Mongolian horses occupying the negative extreme. Anatomically, CV2 mainly expresses the variation in the width and height of the skull (from incisive to orbital areas), with negative part of the axis corresponding to broader and higher skulls which present a rounder braincase and less developed occipital condyles. Variation in the shape of the incisive bone can be noticed as well, with its anterio-ventral extension at the positive extreme of the axis.

Comparisons of the Procrustes variance among the main groups show several significant differences. Przewalski’s and racehorses display the highest variances, whereas draft horses appear as the least morphologically variable (Table 2).

Table 2 Procrustes variance and p-values obtained from the pairwise comparisons between the main groups (significant differences with *; p < 0.05)

Allometry-free shape variation of the modules

Results obtained from the EMMLi approach indicate that the best-supported model of cranial modularity in horses is the Model 4 (mammalian model) [24], with both within-module and between-modules distinct correlation values (see Additional file 4: Table S1). This same model is also supported in each of the main groups independently (see Additional file 4: Tables S2 to S5). In accordance with EMMLi results, the lowest CR value is obtained for the Model 4 (CR = 0.61/p < 0.05; see Additional file 4: Table S7) with an effect size (Z-score = − 22.2) significantly lower than that of other models (see Additional file 4: Table S7). Similar results were obtained from EMMLi analyses computed on curves and landmarks only Additional file 4: Table S6). The lowest CR value for analysis performed on curves and landmarks only is obtained for the Model 4, with a lowest Z-score obtained for the Model 5 but not significantly different from that of the Model 4 (see Additional file 4: Table S7). The Model 4 was thus retained for further analyses.

To describe the shape diversity within modules, CVAs were performed separately on each of the six modules (anterior oral-nasal/AON, orbital/ORB, molar/MOL, basicranium/CB, zygomatic-pterygoid/ZP, cranial vault/CV). For all modules, the first axis distinguishes domestic from Przewalski’s horses (as already observed on the complete skull; Fig. 7). Concerning the second axis, the patterning we observed in the whole-skull dataset is again replicated within the modules AON and ZP (with racehorses occupying one morphological extreme, Mongolian horses occupying the other, and draft horses pooling with feral horses and landraces in an intermediate position). A quite similar pattern can be observed for ORB and MOL, although the groups are less clearly distinguishable. For these four modules, it should however be noted that the two Icelandic horses from our sample pool with racehorses, which differs from the results obtained on the complete skull. The only exceptions to this pattern are the modules CV and CB, for which Shetland ponies occupy the negative extreme of the axis, with draft horses the only other distinguishable group in CB.

Fig. 7
figure 7

Scatter plots of the two first CVs of the CVA performed for each module (a AON; b ORB; c MOL; d ZP; e CV; f CB) on 95% of the allometry-free shape data (corresponding respectively to 29, 23, 27, 33, 25 and 33 PCs) and visualization of the shape changes along CV1 (g) and CV2 (h) (+: extreme positive; −: extreme negative). Crosses represent the centroid value for each group. AON anterior oral-nasal, ORB orbital, MOL molar, CB basicranium, ZP zygomatic-pterygoid, CV cranial vault

Comparisons among modules of the Procrustes variance indicate that CB displays the highest adjusted variance value (3.4e−05) followed by CV (1.4e−05) and ZP (1.1e−05), whereas the lowest variances are obtained for AON (5.6e−06), ORB (6.5e−06) and MOL (7.5e−06).

Comparisons among groups show no significant difference in Procrustes variance for ORB, ZP, CB and CV (Table 3). Concerning other modules, we observe the same tendency found on the complete skull, with the highest variance values obtained in Przewalski’s horses, and to a lesser extent in racehorses, and the lowest variance occurring in draft horses. A unique pattern can be seen in the MOL module, wherein Mongolian horses are the only group to display a variance value as high as that observed in Przewalski’s horses.

Table 3 Procrustes variance within modules and p-values obtained from the pairwise comparisons between the main groups (significant differences with *; p < 0.05)

Modularity and integration

The highest CR effect size (corresponding to lower modular signal) computed on the complete skull was obtained on draft horses (CR = 0.71, effect size = − 9.833), followed by Przewalski’s horses (CR = 0.67, effect size = − 9.835), Mongolian horses (CR = 0.73, effect size = − 9.837) and racehorses (CR = 0.67, effect size = − 9.839). Pairwise comparisons reveal significant differences in the degree of modularity among all the groups. Comparable results are observed on the dataset including landmarks and curves only, with the lowest Z-score obtained for racehorses and the highest one for Przewalski’s horses (see Additional file 5). Concerning the within-module magnitude of integration, the eigenvalue dispersion of the covariance matrix shows the highest degree of morphological integration in CV (cranial vault) and the lowest in CB (basicranium; Table 4). Similar results were also obtained looking at each group separately, except in two cases: in Przewalski’s horses, for which CB is the most highly integrated module, followed by CV; and in Mongolian horses, for which ORB (orbital) and MOL (molar) display a stronger degree of integration than CV. Similar results are obtained when considering landmarks and curves only (see Additional file 5).

Table 4 Eigenvalue dispersion of the covariance matrice indicating the degree of morphological integration within each module

We quantified the magnitude of covariation between adjacent modules using 2B-PLS. Some significant pairwise differences in PLS effect size can be observed, with the strongest degree of integration obtained for the AON/MOL pair, and, to a lesser extent, ZP/CV (Table 5). Similar results are obtained looking at each group separately (see Additional file 6) as well as considering landmarks and curves only (see Additional file 5).

Table 5 Pairwise comparisons of the effect sizes of PLS analyses indicating the degree of morphological integration between the adjacent modules (significant differences with *; p < 0.05)

To visualize integration patterns, we produced scatter plots of the first PLS axes describing covariation between the adjacent modules (Fig. 8). For better visibility, 90% data ellipses were drawn [71]. A common trend in the distribution of the specimens stands out in pairs of modules from the anterior part of the skull (AON/MOL, AON/ORB, MOL/ZP and, to a lesser extent, MOL/ORB) with Mongolian and racehorses occupying distinct positions along the PLS axis. Other groups exhibit more intermediate positions but with draft horses tending to pool with Mongolian horses, and feral horses tending to pool with racehorses. Przewalski’s horses generally occupy almost all of the axis. A different tendency can be observed on the first PLS axis between the ZP and CB modules, with racehorses differing from Przewalski’s horses, whereas the other groups occupy intermediate positions. Finally, the distribution of the specimens along the axes of covariation between posterior modules (ORB/CV, ZP/CV, CB/CV) does not show clearly structured patterns.

Fig. 8
figure 8

Scatter plot of the first PLS axes describing covariation between adjacent cranial modules. rPLS PLS correlation coefficient, TC total covariation. Crosses represent the centroid value for each group. AON anterior oral-nasal, ORB orbital, MOL molar, CB basicranium, ZP zygomatic-pterygoid, CV cranial vault

Angular comparisons between PLS axes reveal no significant differences between groups in most cases, indicating that all the groups share similar covariation patterns in most cases. The only exception concern the Przewalski’s horses, for which the main PLS axis between ORB and CV differs in direction from that of Mongolian horses (88.9°), and for which the main PLS axis between ZP and CB differ from that Mongolian (87.0°) and draft horses (85.8°). The covariation between CV and CB should also be mentioned as the plot reveals the parallel distribution of Przewalski’s and racehorses along the first PLS axis, suggesting different covariation patterns without angular difference.


Cranial shape variation in the horse: the role of allometry and artificial selection in the morphological diversification

Results from this study first demonstrate that allometry strongly contributes to shape diversification in domestic horses, among which considerable size disparity exists. Allometry accounts for a large part of the shape variation in the horse’s skull, with 93% of the shape changes explained by size. Our approach using surface semilandmarks allows for the precise depiction of these changes, insofar as they are mainly located on regions deprived from anatomical landmarks (e.g. cranial vault, nasal bone). Shetland ponies can be differentiated from the other breeds based on their broad skulls displaying a round braincase, concave face and large orbits. These characteristics, which correspond to features generally described as juvenile [72,73,74], have been imputed to differences in the ontogenetic trajectory lengths and slopes among horse breeds, producing comparable skull shapes in adult ponies than in taller horses of younger age [73].

When considering the CVA based on the allometry-free shape, species variation produces a clear differentiation between domestic and Przewalski’s horses which appears largely related to the skull width, and to the shape of the occipital bone. Differences in posture and motion of the neck related to activity could contribute to explain why shape changes appear mainly concentrated in the occipital region within our sample of wild and domestic horses. Indeed, traction or riding have for instance been suggested to cause nuchal enthesopathy [75, 76], as well as differences in the captivity conditions which could have notably involved differences in the feed intake posture [76]. These variations being related to individuals’ life history, further studies comparing free-ranging and captive horses are needed to explore this question more deeply, in order to disentangle the impact of artificial selection and the potential effects of captivity on the shape of the occipital bone.

Among domestic horses, the second axis separates racehorses, which experience the strongest degree of artificial selection in our sample, from Mongolian horses, which live free-ranging lives with a lower degree of human intervention in reproduction and management in comparison to other domestic horses. Mongolian horses display broad skulls with relatively weakly developed occipital condyles, as observed in Przewalski’s horses. Importantly, their skull exhibits a round braincase, suggesting a potential impact of the brain size on the cranial vault shape. Domestication is generally considered to have led to brain size reduction in various domestic animals, notably horses [77]. If Mongolian horses have indeed been subjected to a lower degree of artificial selection, this consideration may explain the relatively larger braincases and shape differences among this group. Draft horses occupy an intermediate position on this axis, which could be related to the lower degree of artificial selection they have undergone on specific morphological features than on body size [78]. Interestingly, American feral horses also occupy an intermediate position. Considering that most contemporary American feral horses stem from light racehorse-types (imported by European settlers from the sixteenth century) [79], this result suggests that the return to wild conditions and the relaxation of artificial selection are accompanied by morphological changes. This would help explain historical observations of various feral horse populations worldwide, from Australian brumbies to West Africa and the Americas, especially concerning skull breadth [80,81,82]. The ability of feral organisms to revert to wild-type features has been observed in several taxa [6, 40,41,42, 83]. Further research including a larger number of feral populations is now needed to assess the response of domestic captive-bred horses to natural selection and free-roaming lifestyle. The possibility to consider feral populations as proxies of wild-type ancestral populations (no longer existing for domestic horses) could improve our understanding of domestication mechanisms by allowing to observe a potential reversal of the effects of domestication [84,85,86]. In the case of dogs, for which the wild ancestor is currently also extinct [87, 88], the interest given to dingoes (as a unique model of feralization because having lived isolated from domestic dogs during around 8000 years [89]) exemplifies the growing attention to the topic of feralization in research on domestication.

Although relatively small sample sizes among our individual draft and racehorse breeds compelled us to group individual draft and racehorse breeds together, the fact that we observed no apparent evidence for structuring of breeds within each of these groups suggests that breed definitions are not the most relevant level of morphological division among domestic horses. Having been subject to crossbreeding and refinement crossing to varying extents, horse breeds exhibit a wide range of within-breed homogeneity [90] and at least some of these breeds are arbitrary from a genetic point of view [36].

The shape variation within racehorses, a domestic group known to experience strong degrees of artificial selection, is as large as the intra-species variation observed within Przewalski’s horses. This demonstrates the ability of artificial selection to produce massive shape diversification over relatively short time frames of a few hundred years [3, 6]. We should caution that, although the sample of Przewalski’s horses included in this study was selected to include a diverse range of life history conditions, with specimens from both captivity and free-roaming conditions, the modern-day representatives of this species originate from a single, small group of founders [29, 31, 32] that have been bred in captivity since their extinction in their original range in the 1960’s [30]. Undoubtedly, this recent history has impacted this group’s morphology. Nonetheless, the fact that we observe a similar within-species variation in our sample of Przewalski’s horses than in a single group of domestic horses is in accordance with what is observed at a wider level in equids [4] and other mammalian taxa [1, 3, 91]. This confirms that the strong artificial selection to which domestic horses have been subjected, aiming to produce various specific conformations and functional features, is a major evolutionary force driving shape diversification.

Cranial modularity and shape diversification

Stasis in patterns, changes in magnitude

Because selection may also alter the developmental and functional processes underlying the generation of variation, we examined the organization of horse cranial covariation. Our results support the idea that a conserved modular system does not limit cranial diversification in horses. The best-supported modular patterning of the skull is the same for all the studied groups and corresponds to that brought out in equids [4] and more generally in therian mammals [24]. Moreover, the vector angles between the main PLS axes reveal that covariation is characterized by similar patterns in all the groups. These results are concomitant with those obtained in various mammalian and avian taxa [1, 6, 27] and contribute to highlight the role of stabilizing selection on functional and developmental processes in maintaining the cranial integration structure [27].

Our results also reveal differences in the intensity of cranial integration among groups. The lowest magnitude of integration is obtained for racehorses whereas draft and Przewalski’s horses provide the highest values. This would appear to indicate that the strong degree of artificial selection to which the race breeds have been subjected may have decreased the magnitude of their cranial integration. This result is consistent with that obtained at the genus level which showed lower magnitudes of integration in domestic than in wild equids [4], although not observed more widely in mammalian domesticates [91]. Taken together, these results suggest that the diversification process related to artificial selection increases modular organization in horses. These relaxed covariation constraints may facilitate the shape diversification observed in the racehorses, and in domestic horses in general [4], without the need to disrupt modularity patterns and inter-modules relationships [26, 27, 92, 93].

Influence of morphological integration on shape variance

The modular structure of the skull may also foster independent variation of each module, resulting in a differential degree of shape variance and of intra and inter-module integration. In this study, the highest variance values were observed for the modules from the neurocranium (basicranium/CB and cranial vault/CV). The fact that these two modules are equally variable across draft horses, Mongolian horses, Przewalski’s horses, and racehorses confirms that this trend for the neurocranium applies across all studied groups. The lowest value of disparity was observed for the anterior oral-nasal (AON).

How the degree of integration influences morphological evolution, and thus impacts morphological disparity, is a lingering question [22, 26, 28]. In general, there are two opposing hypotheses: (1) that high magnitude of integration restricts the variation of individual traits, resulting in low morphological disparity and decreased evolutionary flexibility; (2) that high magnitude of integration promotes variation through coordinated morphological changes among traits, resulting in high morphological disparity and increased evolutionary flexibility. This study does not show a simple relationship between shape variance and magnitude of morphological integration, as it has been outlined in previous studies [22, 91]. Instead, the two most variable modules in our sample, which display respectively the lowest (CB) and highest (CV) eigenvalue dispersions, suggest that morphological integration would both constrain (CB) and facilitate (CV) morphological changes. This confirms that there is probably no simple rule relating morphological evolution and integration [22].

A different tendency entirely emerges in our sample of Przewalski’s horses, for which CB and CV are both the most variable and integrated modules, results which support the idea that strong integration facilitates variation in this taxon. In contrast, all three less variable modules (AON, ORB/orbital and MOL/molar) are the three least integrated in Przewalski’s horses. In a horse species which is currently not under artificial selection, our results are thus consistent with strong degrees of integration facilitating shape variation, in accordance with previous results obtained on dingoes and dogs [28], but differing from those obtained at a wider level in wild and domestic taxa of the genus Equus [4]. To sum up, while our overall results tend to more often support the facilitation hypothesis, discrepancies among our own results, as well as between our findings and those reported in other studies, leave few firm answers as to how the magnitude of integration impacts shape disparity.

Finally, our results also reveal that the three most variable modules (ZP/zygomatic-pterygoid, CB and CV), which are the three most integrated among wild horses, are also those which most poorly covary with other modules (excepting covariation between ZP and CV). This observation suggests that a strong level of modular independence could also be related to high values of shape variance.

The role of function and development in integration patterns

The six-module model coincides with functional groupings within the equine skull: the anterior oral-nasal (AON) and molar (MOL) parts are the main modules involved in the masticatory apparatus; the zygomatic-pterygoid (ZP) module is involved in mastication by comprising jaw muscle attachments; the orbit (ORB) module houses the visual structures; the cranial vault (CV) module provides protection to the brain; the cranial base (CB) module is both involved in supporting the braincase and constitutes in the attachment point between the skull and the axial skeleton. Consequently, examining the magnitude of within and inter-module integration enables us to investigate the patterns of functional relationships within the horse skull. Developmental processes, in particular the mode of ossification and tissue origin, also contribute to modular patterning by causing covariation among the structures derived from each of these origins [25]. AON, MOL and ORB are all derived from a same tissue origin (neural crest) and mode of ossification (dermal bones), although ORB may also include paraxial mesoderm-derived tissue and bone formed by endochondral ossification. The CB module only encompasses endochondral bones of paraxial mesoderm origin. Finally, the ZP and CV modules are composed of both neural crest and paraxial mesoderm derived bones and both dermal and endochondral bones [94, 95].

Antero-posterior patterning clearly emerges from the cranial shape variation and covariation in our analyzed sample of horses. The three most variable, integrated and independent modules in Przewalski’s horses are the three posterior ones (ZP, CV and CB). The relative independence of these posterior modules is further highlighted by our PLS analyses, which show no clear differentiation along the first PLS axes, but significant differences in covariation patterns (i.e. angular differences for ORB/CV and ZP/CB or parallel trajectories for CV/CB between the main PLS axes).

Anterior region

The common developmental origin of AON, MOL and to a lesser extent, ORB (i.e. neural crest-dermal bones), along with the fusion of facial prominences [25], explains strong inter-module relationships between these units, but functional reasons may also influence this pattern. For instance, strong covariation between AON and MOL is consistent with their shared functional role in the masticatory apparatus, as shared muscle attachments and mechanical activities (i.e. mastication) are known to produce covariation among structures [96]. Covariation between these anterior modules is accompanied by low intra-module integration and low shape variability, implying that inter-module integration could constrain the independent shape variation of the modules via developmental and/or functional constraints. This low shape variance observed in AON and MOL is in accordance with the weaker impact of domestication on the facial skeleton in horses than in many other domestic species [77, 97].

Mongolian horses seem to differ from other horse groups concerning the anterior modules with a specific pattern of variation and covariation in the MOL region. Indeed, within-module integration in MOL for Mongolian horses is higher than in CV, which contrasts with the results obtained for the other groups. To some extent, this difference may reflect lower CV integration in Mongolian horses, which, paired with their broader braincase, suggests that relaxed covariation constraints in CV would facilitate potential variations in brain volume. However, our results also reveal for Mongolian horses an unusual high variance value in MOL (paired with an absence of clear shape differentiation from draft and feral horses). In Carnivora, the MOL region has been identified as a strongly integrated and disparate module [22], a pattern linked to strong selective pressures applied on this area of crucial functional importance in a clade with high ecological and dietary diversity. Functional requirements related to mastication for varied diet is thus one potential explanation for the high variance and magnitude of integration in MOL in Mongolian horses.

Posterior region

Both developmental and functional processes also seem to be involved in the modular patterning of the posterior part of the horse skull. In this region, inter-module integration is generally low, which could be explained in part by the diversity of developmental origins across modules in this area of the skull. Strong covariation is, however, observed between ZP and CV, which could be due to their similar developmental pathway (these modules are the only ones to encompass bone deriving from both dermal and endochondral origin, and tissue from both neural crest and paraxial mesoderm).

Functional differentiation may also explain the low degree of inter-module integration in the posterior area of the skull. The brain being the largest organ in the skull of most mammals, its growth is an important driver for CV variation [98]. The evolutionary necessity of CV flexibility could thus explain the relative independence of the neurocranial modules. Domestication in horses is often hypothesized to have produced a reduction in brain size via selection for tameness, potentially modulated by changes in neural crest development [77, 97, 99]. Moreover, artificial selection could have deliberately influenced brain regions in some breeds, in relation to specific uses (e.g. for training or dressage). The high variance value we obtained for CV is thus consistent with variation in brain size linked to selective pressures in domestication. It should be noted that historical CV variance is probably underestimated in our analysis due to the absence of surviving undomesticated conspecifics of E. caballus. Strong selective pressures may also have impacted the Przewalski’s horse during its near-extinction in the twentieth century. Due to their time in captivity, Przewalski’s horses are considered as having been subject to a 14% decrease in cranial volume which would be similar to that of domestic horses [85, 100].

The CV module also emerges as the most highly integrated in our dataset. Whereas a strong degree of integration also characterizes the CV in the genus Equus [4], among mammals as a whole, this area is differentially integrated. The CV is highly integrated in carnivores, while a reduced magnitude of integration in primates [101, 102], may have allowed for the expansion of the brain [24]. Because of the hypothesized link between CV integration and encephalization, further research on a larger sample of feral horses and wild-raised Przewalski’s horses will help clarify the potential impact of free-roaming lifestyle and artificial selective pressures on cranial volume and CV integration.

Highly variable and poorly integrated across domestic horses, CB is the most integrated module in Przewalski’s horses. This finding is consistent with previous observations by Heck et al. [4], who showed that the CB module had the highest degree of integration in wild equids as compared to domesticated horses. A strongly integrated CB is also found more broadly across mammalian taxa [22, 24], but is usually associated with low shape variance [22]. Taken together, our results suggest that the CB module, generally considered as an evolutionarily conserved region [22, 103], has likely been subjected to significant selection pressure in domestic horses and that relaxed covariation constraints may have facilitated morphological diversification. This scenario would explain high CB variance in domestic horses, and the fact that the occipital appears as a main driver of morphological changes across our study groups. The CB module is derived from a single tissue origin (paraxial mesoderm) and is composed of bones using a single mode of ossification (endochondral). As a result, lower integration in this unit is probably caused by functional factors rather than developmental factors. The CB is implicated in two different functions (supporting the brain and connecting the skull to the axial skeleton), and differential selection could have favored functional dissociation of the within-modules traits. As a consequence of domestication and captivity, potential variations in posture [75] or brain size could have occurred, producing increased flexibility in the CB region with relaxed covariation constraints.


This study elucidated microevolutionary mechanisms underpinning phenotypic diversification of domestic horse breeds under artificial selection. We confirmed the ability of artificial selection to produce large amounts of shape diversity, in comparing domestic breeds to the extant wild form, the Przewalski’s horse. As already shown in several taxa, we also found that this drastic diversification did not rely upon changes in modularity patterns but rather upon variations in the magnitude of integration between morphological features. Our results also reveal that strong degrees of artificial selection are associated with lower intensity of integration, suggesting that increased independence of cranial modules facilitates rapid shape changes. A particularly high degree of autonomy was obtained for modules located in the posterior region of the skull, an area involved in holding brain and connecting the skull to the axial skeleton. Because of the potential variations in brain size and head posture associated with domestication, an enhanced need for flexibility in this anatomical region could explain this result. Further studies focused on reintroduced Przewalski’s horses will help identify evolutionary changes over generations under natural selection and disentangle these from plastic signals related to life conditions. Additional research involving a larger number of feral horse populations is also needed to evaluate how adaptive responses to natural selection impact cranial shape variation and covariations in domestic animals reverted into a wild state. Finally, by revealing morphological and microevolutionary responses of domestic horses to artificial selection, our findings may help better understand the domestication process in horses and other large mammals. Future research involving archaeological material from early domestication contexts could allow to track down morphological changes related to human impact on the horse over time.


  1. Curth S, Fischer MS, Kupczik K. Patterns of integration in the canine skull: an inside view into the relationship of the skull modules of domestic dogs and wolves. Zoology. 2017.

    Article  PubMed  Google Scholar 

  2. Darwin C. The variation of animals and plants under domestication. London: J. Murray; 1868.

    Google Scholar 

  3. Drake AG, Klingenberg CP. Large-scale diversification of skull shape in domestic dogs: disparity and modularity. Am Nat. 2010;175:289–301.

    Article  PubMed  Google Scholar 

  4. Heck L, Wilson LAB, Evin A, Stange M, Sánchez-Villagra MR. Shape variation and modularity of skull and teeth in domesticated horses and wild equids. Front Zool. 2018;15:14.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Trut L, Oskina I, Kharlamova A. Animal evolution during domestication: the domesticated fox as a model. BioEssays. 2009;31:349–60.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Young NM, Linde-Medina M, Fondon JW, Hallgrímsson B, Marcucio RS. Craniofacial diversification in the domestic pigeon and the evolution of the avian skull. Nat Ecol Evol. 2017;1:0095.

    Article  Google Scholar 

  7. Hanot P, Herrel A, Guintard C, Cornette R. The impact of artificial selection on morphological integration in the appendicular skeleton of domestic horses. J Anat. 2018;232:657–73.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Denis B. Les races de chevaux en France au XVIIIe siècle. Et les idées relatives à leur amélioration. Situ Rev Patrim. 2012.

    Article  Google Scholar 

  9. Fages A, Hanghøj K, Khan N, Gaunitz C, Seguin-Orlando A, Leonardi M, et al. Tracking five millennia of horse management with extensive ancient genome time series. Cell. 2019.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Weatherby J. An introduction to a general stud book. London: Weather Sons; 1791.

    Google Scholar 

  11. Gould SJ. The structure of evolutionary theory. Cambridge: Harvard University Press; 2002.

    Book  Google Scholar 

  12. Goswami A, Smaers JB, Soligo C, Polly PD. The macroevolutionary consequences of phenotypic integration: from development to deep time. Phil Trans R Soc B. 2014;369:20130254.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Hallgrímsson B, Willmore K, Hall BK. Canalization, developmental stability, and morphological integration in primate limbs. Am J Phys Anthropol. 2002;119:131–58.

    Article  Google Scholar 

  14. Wagner GP. Homologues, natural kinds and the evolution of modularity. Am Zool. 1996;36:36–43.

    Article  Google Scholar 

  15. Cheverud JM. Phenotypic, genetic, and environmental morphological integration in the cranium. Evolution. 1982;36:499–516.

    Article  PubMed  Google Scholar 

  16. Magwene PM. New tools for studying integration and modularity. Evolution. 2001;55:1734–45.

    CAS  PubMed  Google Scholar 

  17. Marroig G, Cheverud JM, Wainwright P. Size as a line of least evolutionary resistance: diet and adaptive morphological radiation in new world monkeys. Evolution. 2005;59:1128–42.

    PubMed  Google Scholar 

  18. Van Valen L. The study of morphological integration. Evolution. 1965;19:347–9.

    Article  Google Scholar 

  19. Wagner GP, Altenberg L. Perspective: complex adaptations and the evolution of evolvability. Evolution. 1996;50:967–76.

    Article  PubMed  Google Scholar 

  20. Olson EC, Miller RL. Morphological integration. Chicago: University of Chicago Press; 1958.

    Google Scholar 

  21. Calabretta R, Nolfi S, Parisi D, Wagner GP. Duplication of modules facilitates the evolution of functional specialization. Artif Life. 2000;6:69–84.

    Article  CAS  PubMed  Google Scholar 

  22. Goswami A, Polly PD. The influence of modularity on cranial morphological disparity in carnivora and primates (Mammalia). PLoS ONE. 2010;5:e9517.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. Martin DP, Van der Walt E, Posada D, Rybicki EP. The evolutionary value of recombination is constrained by genome modularity. PLoS Genet. 2005;1:e51.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Goswami A. Cranial modularity shifts during mammalian evolution. Am Nat. 2006;168:270–80.

    Article  PubMed  Google Scholar 

  25. Hallgrimsson B, Lieberman DE, Young NM, Parsons T, Wat S. Evolution of covariance in the mammalian skull. Novartis Found Symp. 2007;284:164–90.

  26. Marroig G, Shirai LT, Porto A, de Oliveira FB, De Conto V. The evolution of modularity in the mammalian skull II: evolutionary consequences. Evol Biol. 2009;36:136–48.

    Article  Google Scholar 

  27. Porto A, de Oliveira FB, Shirai LT, De Conto V, Marroig G. The evolution of modularity in the mammalian skull I: morphological integration patterns and magnitudes. Evol Biol. 2009;36:118–35.

    Article  Google Scholar 

  28. Parr WCH, Wilson LAB, Wroe S, Colman NJ, Crowther MS, Letnic M. Cranial shape and the modularity of hybridization in dingoes and dogs; hybridization does not spell the end for native morphology. Evol Biol. 2016;43:171–87.

    Article  Google Scholar 

  29. Boyd L, Houpt KA. Przewalski’s horse: the history and biology of an endangered species. Albany: SUNY Press; 1994.

    Google Scholar 

  30. International Union for Conservation of Nature (IUCN). IUCN red list of threatened species. Version 20111. 2012.

  31. Bowling AT, Zimmermann W, Ryder O, Penado C, Peto S, Chemnick L, et al. Genetic variation in Przewalski’s horses, with special focus on the last wild caught mare, 231 Orlitza III. Cytogenet Genome Res. 2003;102:226–34.

    Article  CAS  PubMed  Google Scholar 

  32. Goto H, Ryder OA, Fisher AR, Schultz B, Pond K, et al. A massively parallel sequencing approach uncovers ancient origins and high genetic variability of endangered Przewalski’s horses. Genome Biol Evol. 2011;3:1096–106.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Van Dierendonck MC, Wallis de Vries MF. Ungulate reintroductions: experiences with the takhi or Przewalski horse (Equus ferus przewalskii) in Mongolia. Conserv Biol. 1996;10:728–40.

    Article  Google Scholar 

  34. Der Sarkissian C, Ermini L, Schubert M, Yang MA, Librado P, Fumagalli M, et al. Evolutionary genomics and conservation of the endangered Przewalski’s horse. Curr Biol CB. 2015;25:2577–83.

    Article  CAS  Google Scholar 

  35. Gaunitz C, Fages A, Hanghøj K, Albrechtsen A, Khan N, Schubert M, et al. Ancient genomes revisit the ancestry of domestic and Przewalski’s horses. Science. 2018;360:111–4.

    Article  CAS  PubMed  Google Scholar 

  36. Sponenberg DP. Genetic resources and their conservation. In: Bowling AT, Ruvinsky A, editors. The genetics of the horse. Wallingford: CABI; 2000. p. 387–438.

    Chapter  Google Scholar 

  37. Petersen JL, Mickelson JR, Rendahl AK, Valberg SJ, Andersson LS, Axelsson J, et al. Genome-wide analysis reveals selection for important traits in domestic horse breeds. PLoS Genet. 2013;9:e1003211.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Hendricks BL. International encyclopedia of horse breeds. Norman: University of Oklahoma Press; 2007.

    Google Scholar 

  39. Sponenberg DP, Christman C. A conservation breeding handbook. American Livestock Breeds Conservancy; 1995. Accessed 21 Apr 2017.

  40. Birks J, Kitchener A. The distribution and status of the polecat Mustela putorius in Britain in the 1990s. Vincent Wildl Trust Lond. 1999;152.

  41. Lynch JM, Hayden TJ. Genetic influences on cranial form: variation among ranch and feral American mink Mustela vison (Mammalia: Mustelidae). Biol J Linn Soc. 1995;55:293–307.

    Article  Google Scholar 

  42. Pocock R. Ferrets and polecats. Scott Nat. 1932;196:97–108.

    Google Scholar 

  43. Sol D. Artificial selection, naturalization, and fitness: Darwin’s pigeons revisited. Biol J Linn Soc. 2008;93:657–65.

    Article  Google Scholar 

  44. Hanot P, Guintard C, Lepetz S, Cornette R. Identifying domestic horses, donkeys and hybrids from archaeological deposits: a 3D morphological investigation on skeletons. J Archaeol Sci. 2017;78:88–98.

    Article  Google Scholar 

  45. Wiley DF, Amenta N, Alcantara DA, Ghosh D, Kil YJ, Delson E, et al. Evolutionary morphing. In: Proceedings of IEEE visualization 2005. Minneapolis: IEEE; 2005. p. 431–8.

  46. Schlager S. Morpho: calculations and visualisations related to geometric morphometrics. 2016.

  47. Rohlf FJ, Slice D. Extensions of the procrustes method for the optimal superimposition of landmarks. Syst Biol. 1990;39:40–59.

    Google Scholar 

  48. Baylac M, Frieß M. Fourier descriptors, procrustes superimposition, and data dimensionality: an example of cranial shape analysis in modern human populations. In: Slice DE, editor. Modern morphometrics in physical anthropology. New York: Springer Science & Business Media; 2005. p. 145–65.

    Chapter  Google Scholar 

  49. Jolliffe IT. Principal component analysis. 2nd ed. New York: Springer; 2002.

    Google Scholar 

  50. Krzanowski WJ. Principles of multivariate analysis: a user’s perspective. New York: Oxford University Press, Inc.; 1988.

    Google Scholar 

  51. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57:289–300.

    Google Scholar 

  52. Monteiro LR. Multivariate regression models and geometric morphometrics: the search for causal factors in the analysis of shape. Syst Biol. 1999;48:192–9.

    Article  CAS  PubMed  Google Scholar 

  53. Zelditch ML, Swiderski DL, Sheets HD. Geometric morphometrics for biologists: a primer. Cambridge: Academic Press; 2012.

    Google Scholar 

  54. Bardua C, Wilkinson M, Gower DJ, Sherratt E, Goswami A. Morphological evolution and modularity of the caecilian skull. BMC Evol Biol. 2019;19:30.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Adams DC, Otárola-Castillo E. geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods Ecol Evol. 2013;4:393–9.

    Article  Google Scholar 

  56. Adams DC. Evaluating modularity in morphometric data: challenges with the RV coefficient and a new test measure. Methods Ecol Evol. 2016;7:565–72.

    Article  Google Scholar 

  57. Goswami A, Finarelli JA. EMMLi: a maximum likelihood approach to the analysis of modularity. Evolution. 2016;70:1622–37.

    Article  PubMed  Google Scholar 

  58. Adams DC, Collyer ML. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution. 2019;73:2352–67.

    Article  PubMed  Google Scholar 

  59. Goswami A, Watanabe A, Felice RN, Bardua C, Fabre A-C, Polly PD. High-density morphometric analysis of shape and integration: the good, the bad, and the not-really-a-problem. Integr Comp Biol. 2019;59:669–83.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Marshall AF, Bardua C, Gower DJ, Wilkinson M, Sherratt E, Goswami A. High-density three-dimensional morphometric analyses support conserved static (intraspecific) modularity in caecilian (Amphibia: Gymnophiona) crania. Biol J Linn Soc. 2019;126:721–42.

    Article  Google Scholar 

  61. Cardini A. Less tautology, more biology? A comment on “high-density” morphometrics. Zoomorphology. 2020;139:513–29.

    Article  Google Scholar 

  62. Bookstein FL. Morphometric tools for landmark data: geometry and biology. Cambridge: Cambridge University Press; 1997.

    Google Scholar 

  63. Bookstein FL, Gunz P, Mitterœcker P, Prossinger H, Schæfer K, Seidler H. Cranial integration in Homo: singular warps analysis of the midsagittal plane in ontogeny and evolution. J Hum Evol. 2003;44:167–87.

    Article  PubMed  Google Scholar 

  64. Rohlf FJ, Corti M. Use of two-block partial least-squares to study covariation in shape. Syst Biol. 2000;49:740–53.

    Article  CAS  PubMed  Google Scholar 

  65. Klingenberg CP. MorphoJ: an integrated software package for geometric morphometrics. Mol Ecol Resour. 2011;11:353–7.

    Article  PubMed  Google Scholar 

  66. Pavlicev M, Cheverud JM, Wagner GP. Measuring morphological integration using eigenvalue variance. Evol Biol. 2009;36:157–70.

    Article  Google Scholar 

  67. Wagner GP. On the eigenvalue distribution of genetic and phenotypic dispersion matrices: evidence for a nonrandom organization of quantitative character variation. J Math Biol. 1984;21:77–95.

    Article  Google Scholar 

  68. Machado FA, Hubbe A, Melo D, Porto A, Marroig G. Measuring the magnitude of morphological integration: the effect of differences in morphometric representations and the inclusion of size. Evolution. 2019;73:2518–28.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Melo D, Garcia G, Hubbe A, Assis AP, Marroig G. EvolQG—an R package for evolutionary quantitative genetics. F1000Research. 2016.

    Article  PubMed Central  Google Scholar 

  70. Wickham H. ggplot2: elegant graphics for data analysis. Berlin: Springer; 2016.

    Book  Google Scholar 

  71. Fox J, Weisberg S. An R companion to applied regression. Thousand Oaks: Sage Publications; 2018.

    Google Scholar 

  72. Gould SJ. The panda’s thumb: more reflections in natural history. New York: WW Norton & Company; 1992.

    Google Scholar 

  73. Heck L, Sanchez-Villagra MR, Stange M. Why the long face? Comparative shape analysis of miniature, pony, and other horse skulls reveals changes in ontogenetic growth. PeerJ. 2019;7:e7678.

    Article  PubMed  PubMed Central  Google Scholar 

  74. Wayne RK. Cranial morphology of domestic and wild canids: the influence of development on morphological change. Evolution. 1986;40:243–61.

    Article  PubMed  Google Scholar 

  75. Bendrey R. An analysis of factors affecting the development of an equid cranial enthesopathy. Vet Ir Zootech. 2008;41:25–31.

    Google Scholar 

  76. Taylor WTT, Bayarsaikhan J, Tuvshinjargal T. Equine cranial morphology and the identification of riding and chariotry in late Bronze Age Mongolia. Antiquity. 2015;89:854–71.

    Article  Google Scholar 

  77. Clutton-Brock J. A natural history of domesticated mammals. Cambridge: Cambridge University Press; 1999.

    Google Scholar 

  78. Petersen JL, Mickelson JR, Cothran EG, Andersson LS, Axelsson J, Bailey E, et al. Genetic diversity in the modern horse illustrated from genome-wide SNP data. PLoS ONE. 2013;8:e54997.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Philipps D. Wild horse country: the history, myth, and future of the mustang, America’s Horse. New York: WW Norton & Company; 2017.

    Google Scholar 

  80. Forbes C. Australia on horseback. Stuttgart: Macmillan Publishers Aus; 2014.

    Google Scholar 

  81. Law R. The horse in West African history: the role of the horse in the societies of pre-colonial West Africa. New York: Routledge; 1980.

    Google Scholar 

  82. Roe FG. The indian and the horse. Norman: University of Oklahoma Press; 1976.

  83. Trut LN. Early canid domestication: the farm-fox experiment: foxes bred for tamability in a 40-year experiment exhibit remarkable transformations that suggest an interplay between behavioral genetics and development. Am Sci. 1999;87:160–9.

    Article  Google Scholar 

  84. Brisbin I. The ecology of animal domestication: its relevance to man’s environmental crises–past, present and future. Assoc Southeast Biol Bull. 1974;21:3–8.

    Google Scholar 

  85. O’Regan HJ, Kitchener AC. The effects of captivity on the morphology of captive, domesticated and feral mammals. Mammal Rev. 2005;35(3–4):215–30.

    Article  Google Scholar 

  86. Price EO. Animal domestication and behavior. Wallingford: Cabi; 2002.

    Book  Google Scholar 

  87. Freedman AH, Gronau I, Schweizer RM, Vecchyo DO-D, Han E, Silva PM, et al. Genome sequencing highlights the dynamic early history of dogs. PLoS Genet. 2014;10:e1004016.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  88. Frantz LAF, Mullin VE, Pionnier-Capitan M, Lebrasseur O, Ollivier M, Perri A, et al. Genomic and archaeological evidence suggest a dual origin of domestic dogs. Science. 2016;352:1228–31.

    Article  CAS  PubMed  Google Scholar 

  89. Zhang S, Wang G-D, Ma P, Zhang L, Yin T-T, Liu Y, et al. Genomic regions under selection in the feralization of the dingoes. Nat Commun. 2020;11:671.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  90. Leroy G, Callède L, Verrier E, Mériaux J-C, Ricard A, Danchin-Burge C, et al. Genetic diversity of a large set of horse breeds raised in France assessed by microsatellite polymorphism. Genet Sel Evol. 2009;41:5.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  91. Wilson LA, Balcarcel A, Geiger M, Heck L, Sánchez-Villagra MR. Modularity patterns in mammalian domestication: assessing developmental hypotheses for diversification. Evol Lett. 2021;5:385–96.

    Article  PubMed  PubMed Central  Google Scholar 

  92. Brassard C, Merlin M, Guintard C, Monchatre-Leroy E, Barrat J, Callou C, et al. Interrelations between the cranium, the mandible and muscle architecture in modern domestic dogs. Evol Biol. 2020;47:308–24.

    Article  Google Scholar 

  93. Hansen TF, Houle D. Measuring and comparing evolvability and constraint in multivariate characters. J Evol Biol. 2008;21:1201–19.

    Article  CAS  PubMed  Google Scholar 

  94. Couly GF, Coltey PM, Le Douarin NM. The triple origin of skull in higher vertebrates: a study in quail-chick chimeras. Development. 1993;117:409–29.

    Article  CAS  PubMed  Google Scholar 

  95. Kuratani S. Craniofacial development and the evolution of the vertebrates: the old problems on a new background. Zool Sci. 2005;22:1–19.

    Article  Google Scholar 

  96. Zelditch ML, Mezey J, Sheets HD, Lundrigan BL, Garland T. Developmental regulation of skull morphology II: ontogenetic dynamics of covariance. Evol Dev. 2006;8:46–60.

    Article  PubMed  Google Scholar 

  97. Wilkins AS, Wrangham RW, Fitch WT. The “Domestication Syndrome” in mammals: a unified explanation based on neural crest cell behavior and genetics. Genetics. 2014;197:795–808.

    Article  PubMed  PubMed Central  Google Scholar 

  98. Jiang X, Iseki S, Maxson RE, Sucov HM, Morriss-Kay GM. Tissue origins and interactions in the mammalian skull vault. Dev Biol. 2002;241:106–16.

    Article  CAS  PubMed  Google Scholar 

  99. Librado P, Gamba C, Gaunitz C, Sarkissian CD, Pruvost M, Albrechtsen A, et al. Ancient genomic changes associated with domestication of the horse. Science. 2017;356:442–5.

    Article  CAS  PubMed  Google Scholar 

  100. Röhrs M, Ebinger P. Sind Zooprzewalskipferde Hauspferde. Berl Münch Tierärztl Wochenschr. 1998;111:273–80.

    PubMed  Google Scholar 

  101. Ackermann RR, Cheverud JM. Detecting genetic drift versus selection in human evolution. Proc Natl Acad Sci. 2004;101:17946–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  102. Cheverud JM. Developmental integration and the evolution of pleiotropy. Integr Comp Biol. 1996;36:44–50.

    Google Scholar 

  103. Wesley-Hunt GD, Flynn JJ. Phylogeny of the Carnivora: basal relationships among the carnivoramorphans, and assessment of the position of ‘Miacoidea’relative to Carnivora. J Syst Palaeontol. 2005;3:1–28.

    Article  Google Scholar 

Download references


This work was funded by a postdoctoral fellowship granted by the FYSSEN foundation. We are grateful to the researchers, curators and collection technicians who allowed us to access to the reference specimens: Olivier Pauwels (IRSNB-Bruxelles); Manuel Comte (ONIRIS-Nantes); Marion Bindé, David Cochard and Arnaud Lenoble (PACEA-Bordeaux); the TAKH Association (Hures-la-Parade); Joséphine Lesur, Christine Lefèvre, Aurélie Verguin and Céline Bens (MNHN-Paris); Joseph Cook and Jon. Dunnum (Museum of Southwestern Biology, UNM-Albuquerque), Scott Bender (Navajo Technical University), and Michael Stache (Central Natural Science Collections, University of Halle-Wittenberg). Finally, we would like to express our thanks to the two anonymous reviewers for their constructive comments on the manuscript.


Open Access funding enabled and organized by Projekt DEAL.

Author information

Authors and Affiliations



The study was designed by PH. PH, JB, CG, AH, EM, RS and WT collected the data. PH performed the analysis and wrote the main manuscript with inputs from WT, JB, CG, AH, EM and RS. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Pauline Hanot.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

List of the specimens.

Additional file 2.

Landmark definition.

Additional file 3.

Analysis of sexual dimorphism.

Additional file 4.

Study of the modularity structure.

Additional file 5.

Integration and modularity analyses on landmarks and curves only.

Additional file 6.

Pairwise comparisons of the effect sizes of PLS analyses for each group.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Hanot, P., Bayarsaikhan, J., Guintard, C. et al. Cranial shape diversification in horses: variation and covariation patterns under the impact of artificial selection. BMC Ecol Evo 21, 178 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: