Ecomorphology of the cervid calcaneus as a proxy for paleoenvironmental reconstruction

1 INTRODUCTION

Accurate paleoenvironmental reconstructions are important components of the testing and refinement of hypotheses regarding patterns in mammalian evolution, including that of our own species (Curran, 2012; DeGusta & Vrba, 2003; Potts, 1998; Reed, 1997). In order to construct high-resolution environmental frameworks, it is imperative to use a combination of different proxies in contrast with each other (Bishop et al., 2006). Although taphonomic factors dictate the availability of specific paleoenvironmental proxies at paleontological and archeological sites, the study of vertebrate remains often plays an important role in habitat reconstructions (Blain et al., 2014; Bobe et al., 2007; Geraards et al., 1986). This is often accomplished by looking at the presence or absence of certain taxa in the fossil record that are indicative of specific ecological conditions (Shipman & Harris, 1988; Vrba, 1975). A potential weakness of this approach, however, is that it assumes that the ecological preferences of extinct species were similar to those of their closest living relatives (Andrews, 1995). This is not necessarily a valid assumption (Dodd & Stanton, 1990).

Ecomorphology can circumvent this problem by examining the functional morphology of skeletal elements in a vertebrate group as it relates to ecological variables (Andrews & Hixson, 2014). Comparing morphological patterns in recent species with those observed in extinct forms allows us to make inferences about probable adaptations in fossils (Andrews & Hixson, 2014; DeGusta & Vrba, 2003). Although criticized (Klein et al., 2010), ecomorphology is often considered a taxon-free method, as it can be applied on mammalian remains that have not been identified beyond the family level (Kovarovic & Andrews, 2007). As a method it can therefore provide a more direct reconstruction of past habitats, as it considers how animals functioned in their ecosystem instead of simply which taxa were present at a site (Curran, 2009).

Paleoecologically focused ecomorphological studies have been conducted on a number of mammalian groups, including canids (Meloro & Louys, 2011), ursids (Figuerido et al., 2009), rodents (Fernandez & Campomanes, 2003), bats (Stimpson, 2010), equids (Schellhorn, 2009; Scott, 2004), suids (Bishop, 1994; Bishop et al., 2006), and—especially relevant to this study—ruminants (Barr, 2017; Kappelman, 1988; Kappelman et al., 1997; Plummer et al., 2008; Plummer & Bishop, 1994; Schellhorn, 2009; Scott, 2004). Within the latter group, the focus has been mostly on the Bovidae. Outside the work by Curran (2009, 2012, 2015, 2018), the Cervidae have received comparatively little attention, despite their ecological similarity and phylogenetic relatedness to the Bovidae (Geist, 1998; Janis, 2007). Although the capacity of deer to cope with more extreme environments is not of the same magnitude as that of bovids (Geist, 1998; Janis, 2007), this family is characterized by high species diversity and adapted to a range of environments (Putman & Flueck, 2011). As cervids are a common taxonomic group in European and Asian Pleistocene assemblages (Kurten, 1968), they have the potential to serve as valuable paleoecological proxies in the same way as bovids have for African sites (Forrest et al., 2018; Plummer & Bishop, 1994).

This article presents newly developed ecomorphological models for the cervid calcaneus. This element was not only selected because of its high chance of survival in the fossil record (Borerro, 1990), but also because it represents an important component of the locomotor apparatus, and because its morphology is thought to be primarily constrained by environmental factors (Curran, 2012; Köhler, 1993; Polly, 2007). Although the latter assumption is further explored in the next section, it suggests that the calcaneus can be used to predict habitat preference in extinct species. By adding fossil specimens of unknown ecological affinity to a training set of extant specimens, it should be possible to make inferences about the nature of Pleistocene cervid environments.

1.1 Theoretical basis and functional framework

The connection between cervid morphology and the environment finds its basis in the idea that predator evasion strategy is the main selective constraint on the locomotor system in artiodactyls (Curran, 2009; Geist, 1998). Ruminants especially have a range of methods to overcome encounters with predators (Geist, 1998; Ralls, 1974), but once detected, flight is the chief anti-predatory counter measure (Geist, 1998). Consequently, such animals have evolved several modes of locomotion, employed during flight. Despite limited data on how locomotion is precisely related to evasion tactics, there are five principal locomotor behaviors linked to predator evasion in ruminants: leaping (vertical jump), bounding (long horizontal jump), stotting (a bouncing gait), prancing (type of pronounced and exaggerated high step), and tacking or zig-zag running (sharp turn that suddenly changes course) (Caro, 1986, 1994; Caro et al., 2004). Prancing and stotting act more as signals to wrong-foot the predator than as true locomotor adaptations for flight (Caro, 1994; Caro et al., 2004), but tacking is the preferred flight mode among species living in open habitats (Caro et al., 2004). Leaping is more prevalent when fleeing animals have to deal with rugged topography or elaborate vegetation (Caro, 1994).

These observations support the idea that predator evasion and locomotor strategy are strongly connected to the characteristics of ruminant habitats (Barr, 2014a; Jenkins & Camazine, 1977; Kappelman, 1988; Scott, 2004). Because each locomotor strategy has certain biomechanic requirements, differences in habitat structure are predicted to be reflected in the morphology of the limb bones (Barr, 2018; DeGusta & Vrba, 2003, 2005a, 2005b; Kovarovic & Andrews, 2007; Schellhorn, 2009). Ruminants under predation pressure in dense vegetation have limb bones adapted for saltatorial evasion, allowing them to optimally jump over obstacles in the landscape, to maximize the distance between themselves and their pursuer (Curran, 2009; Kappelman, 1988; Köhler, 1993; Leinders, 1979). Ruminants that evade predators in open vegetation are thought to have limb bones adapted for cursoriality, allowing them to quickly generate speed and accommodate for the forces associated with quick tacking movements (Kappelman, 1988; Köhler, 1993; Leinders, 1979). Although most of this ecomorphological groundwork is based on bovids, it is thought that these principles also apply to cervids (Curran, 2009, 2012, 2015; Kovarovic, 2004).

There are multiple ways in which these adaptations are thought to be expressed in the morphology of the cervid calcaneus. The calcaneus, which together with the astragalus functions as a hinge point between the tibia and the metatarsus (the hock joint) (Barr, 2014b), acts as a lever for the triceps surae muscles that insert into the tuber calcanei via the Achilles tendon (Galvez-Lopez & Casinos, 2012). Being mostly restricted to movement in the sagittal plane (Schaeffer, 1947), the calcaneus pushes the limb against the ground and causes the animal to advance by acting as a moment arm of the triceps surae muscles (Alexander, 1983; Barr, 2018; Curran, 2012; Galvez-Lopez & Casinos, 2012). As part of this mechanism, the length of the lever arm of this muscle group (the tuber calcanei) to a large degree determines the force produced by the posterior limbs during locomotion. Increasing its length increases power but will make the joint move slower, a trait considered adaptive in saltatorial species (Curran, 2009, 2012) (Figure 1). A shorter tuber will lack power, but allows for quicker movement and acceleration, a trait optimal for cursoriality (Curran, 2009).

The rest position of the calcaneus is also thought to differ between cursorial and saltatorial forms. If the tuber calcanei is positioned more vertically relative to the astragalus and the cubonavicular, the distance covered by the lever arm is shorter, but less powerful (Curran, 2012; Polly, 2008). This adaptation is associated with cursorial animals of open environments that strongly increase their number of paces when accelerating (Curran, 2012; Gambardyan, 1974; Geist, 1998). When the tuber has a more horizontal position, more power is generated when the muscles contract, a trait advantageous in saltatorial species (Curran, 2009, 2012). Morphologically, this is thought to be expressed in a more oblique orientation of the tuber relative to the articular surfaces of the calcaneus in cursorial species (Curran, 2012). Additionally, there are indications that in cursorial forms the calcaneus is more tightly locked to the astragalus, due to a higher relief in the ridges and grooves on the articular surfaces between the calcaneus and astragalus and because the articular facet supporting the malleolus is thought to be larger and more developed (Polly, 2007; Scarborough et al., 2016). This way movement is minimized in the parasagittal plane, so as to provide optimal stabilization when running (Polly, 2007).

To translate this functional framework into an ecomorphological model, our study builds on earlier studies conducted by Curran (2009, 2012, 2015, 2018), who found correlations between the morphology of the calcaneus and habitat type, using a geometric morphometrics (GMM) model based on 10 landmarks. In her work, Curran (2012) saw evidence that cervids associated with open vegetation have calcanei adapted for more cursorial movement and vice versa. Although Curran (2009, 2012) observed some variation in mobility of the hock joint elements, with closed environment species having a wider sustentaculum than open adapted species, the primary morphological variation noted between open and closed environment specimens was a difference in the rest position of the calcaneus. More horizontally oriented calcanei were found in closed-adapted cervids and more vertically positioned calcanei in open environment forms (Curran, 2012). The models presented in our article are an extension of those pioneered by Curran (2009, 2012, 2015). Using an alternative methodology (i.e., a different data acquisition protocol, alternative habitat/functional group categorization, an altered set of landmarks, an extant dataset less skewed toward North and South American species, and different methods of analysis), this study aims to evaluate some of the ecomorphological correlations found by Curran (2009, 2012, 2015) and potentially establish new ones.

Based on the ecomorphological models and theoretical framework outlined above, we have summarized three functional hypotheses related to the calcaneus (Table 1). The viability of these hypotheses and the extent to which specific morphotypes can be correlated with locomotor behaviors and environmental parameters are addressed in this article.

TABLE 1. Functional hypotheses for the calcaneus

Hypothesis 1	The relative length of the tuber calcanei is predicted to be greater in species adapted to closed environments and shorter in species adapted to open environments.
Hypothesis 2	The orientation of the tuber and articular facets between calcaneus and astragalus are predicted to be more oblique in species adapted to more open environments and less oblique in species of closed environments.
Hypothesis 3	The articular surface supporting the malleolus is predicted to be larger and more developed in species adapted to open environments and shorter in species adapted to closed environments.

2 MATERIALS

Our models were mainly tested using extant cervid data collected by the authors at several European and American institutes (see Supplementary Data). In addition, 3D data from a small number of specimens were shared by the Max Planck Institute of Evolutionary Anthropology (Niven et al., 2009) and the Virtual Zooarchaeology of the Arctic project (Betts et al., 2011; Maschner et al., 2011). In total, 129 calcanei of extant individuals were studied, belonging to 26 species and 6 tribes (Table 2). The composition of the dataset was to an extent a function of the availability of species at the different institutes and sample sizes were in the same range as many similar studies (Curran, 2009; Schellhorn, 2009; Weinand, 2005). A maximum of 16 specimens were scanned per element for each species, to prevent the dataset from becoming skewed toward more common taxa. Although captivity is recognized to affect morphology in certain mammals (O'Regan & Kitchener, 2005), a number (n = 23) of known captive specimens were included to maximize diversity and sample size. The sample of captive individuals per species was too small for a robust analysis of the potentially confounding effects of this factor, but preliminary tests (not included here) conducted on six Axis axis specimens and six Rangifer tarandus specimens, suggested that at least for these species the effects are limited.

Both sexes were included in the sample and only adult, nonpathological specimens were analyzed. Skeletons were considered adult if all epiphyses were fused and all teeth erupted. When possible, the left calcaneus was selected, but in some cases scanned specimens had to be virtually mirrored using Meshlab 2.0 (Cignoni et al., 2008).

3 METHODS AND STATISTICAL ANALYSES

To emphasize the fact that different morphotypes vary as a result of functional differences related to locomotion (Barr, 2014a, 2014b), we used a combined functional/ecological category system, rather than a purely ecological one (Curran, 2012; DeGusta & Vrba, 2003; Kovarovic & Andrews, 2007). Although this difference in categorization was mainly theoretical, in cases where species had an atypical locomotor strategy to cope with certain environments (e.g., species of large body size lacking a capacity for saltatorial locomotion), this had a potential effect on how taxa were categorized. Our approach should be considered an attempt to take the complex interactions between habitat and locomotion into account, but what is ultimately needed for a more complete understanding of the relationship between locomotion and ecology is a series of experimental studies on living specimens (Wainwright, 1991). In the absence of such data, we assigned the cervid calcanei into five habitat/locomotor types based on different sources from the literature (see Table 2), a number chosen to avoid oversimplifying actual habitat variation, but maximize the ability of our statistical models to classify specimens into the correct category (Curran, 2009).

Type A animals (57 specimens) were characterized by a cursorial locomotor strategy, and preference for open environments. However, as cervids are rarely found in truly open, coverless environments (Geist, 1998), this category united taxa adapted to open woodland and the few species found in open plains. Open woodland was defined as an area of trees with an open canopy of 40% or less closure (Thomas & Packham, 2007). Type B species (35 specimens) were saltatorial forms, adapted to closed environments. This included animals of forest types, ranging from closed woodland (defined as an area of trees with more than 40% canopy closure; Thomas & Packham, 2007) to tropical rainforest. Type C species (eight specimens) were intermediate between Type A and Type B species and generalists adapted to open- or closed environments. Although some of the species in this group may be behaviorally closer to either the Type A or Type B forms, they lack the endurance of real cursorial species but are also less adept at leaping than saltatorial species. Because of this reason, and because of the limited number of specimens in this group, it was neither appropriate to add them to one of the other categories, nor to split them into multiple subcategories.

There were two additional types that had either a more cursorial or saltatorial strategy, but merited their own separate categories for being associated with very specific environmental conditions. The Type D group (eight specimens) included high altitude species from mountainous environments and may require a combination of saltatorial and cursorial traits (Geist, 1998). Type E cervids (21 specimens) were defined as species of open wetlands. They are generally considered cursorial forms, but the unusual structure of these habitat types may also favor certain saltatorial traits (Curran, 2009).

Although many authors have made use of linear measurements to quantify osteomorphological traits in ecomorphology (Bishop et al., 2006; Kappelman, 1988; Kovarovic & Andrews, 2007; Scott, 2004), recent studies (Brophy et al., 2014; Cucchi et al., 2009, 2011; Forrest et al., 2018) have demonstrated that GMM have advantages over linear measurement when applied on artiodactyl morphology. We therefore followed Curran (2009, 2012) and opted for a landmark-based 3D GMM approach to quantify shape. This allowed for an efficient control of isometric size effects (Viscosi & Cardini, 2011; Zelditch et al., 2004) through the application of generalized procrustes analysis (GPA) on the raw coordinate data (see Bookstein, 1991; Zelditch et al., 2004). It also permitted for a better preservation of object geometry in the measured data (Baab et al., 2012; Rohlf & Marcus, 1993; Slice, 2005), a distinct advantage when quantifying subtle shape differences inherent to the morphologically conservative cervids (Kurten, 1968; Perez et al., 2007; Strand Vidarsdottir et al., 2002). And thirdly, GMM allowed us to visualize shape differences using thin plate spline transformation grids that describe morphological variation as distortions in a grid (Bookstein, 1991; Zelditch et al., 2004).

We collected raw data on skeletal elements in the form of 3D surface scans generated with a NextEngine 2020i laser scanner. Using the associated software Scanstudio HD 1.3.2, individual scan divisions were cleaned and fused into 3D objects, after which landmarks were virtually placed on the objects with Landmark editor 3.0 (Wiley et al., 2005). Six landmarks were placed at discrete anatomical loci on the calcaneus (Figure 2). The positioning of the landmarks was inspired by earlier ecomorphological models (Curran, 2009, 2012, 2015; Kovarovic & Andrews, 2007), but also chosen because of their relevance to the functional hypotheses. The resulting coordinate data were then exported from Landmark editor and uploaded in Morphologika 2.5 (O'Higgins & Jones, 2006), where a GPA was conducted on the data. From this software, we generated an output-file with procrustes residuals and (log) centroid size as a size measure, for further statistical analysis in PAST 2.17 (Hammer et al., 2001).

Following earlier ecomorphological studies (Bignon et al., 2005; Figuerido et al., 2009; Forrest et al., 2018), we made use of principal components analysis (PCA) to simplify descriptions of variation between specimens and groups and as a primary method to explore morphological variation (Dryden & Mardia, 1998; Zelditch et al., 2004). This reduced the chances of overfitting the model, a problem that can arise when, for example, linear discriminant analyses are conducted on smaller datasets with unequally balanced samples (Kovarovic et al., 2011). Nevertheless, to emphasize between group differences, we made use of between groups PCA (bgPCA), where eigenvectors are derived from the variance–covariance matrix of the group means instead of all the data-points (Seetah et al., 2012). PCA was in a first instance also used to assess intraobserver error. To test the repeatability of the digitization procedure, we used an adjusted version of a protocol by Adriaens (2007), which consisted of randomly selecting and scanning five specimens and landmarking each specimen five separate times using our landmarking procedure. If the replicated specimens clustered tightly together on the first two axes of a PCA conducted on the dataset, the digitization error was considered low (Adriaens, 2007). The same procedure (but without rescanning the specimens) was repeated on another five specimens, to explore for potential error in repeating the landmarking protocol.

To assess statistical significance (p < .05) of cluster separations in the PCAs, an approach was taken that involved conducting a nonparametric MANOVA (npMANOVA) on the relevant PC scores (Hou et al., 2021; Marramà & Kriwet, 2017; Polly et al., 2013; Schutz et al., 2009). We opted for a permutational test, as assumptions required for parametric testing—including across-group homogeneity of variance–covariance matrices and normal distribution— are not necessarily met by highly dimensional data resulting from GMM (Cardini et al., 2015; Lopez-Lazaro et al., 2018). NpMANOVA tests were run on a data matrix including all relevant principal component scores, as indicated by a scree plot of the eigenvalue distribution (Jackson, 1993). When between-group differences were significant, the overall npMANOVA was followed by a post hoc test in the form of pairwise npMANOVAs between all pairs of groups, to asses which groups differed significantly. Bonferroni corrected probabilities were reported to adequately control for Type I errors (Dunn, 1961, 1964).

In addition, we accounted for a number of potentially confounding factors. To assess the effect of allometry, we regressed PC scores against the natural log of centroid size (Cucchi et al., 2011; Owen, 2013), using the ordinary least squares algorithm (following Kilmer & Rodriguez, 2016). If significant correlations between size and shape were found in the regressions, this was considered indicative of an allometric effect (Zelditch et al., 2004). A regression was also computed of the average centroid size per species against the natural log of average body mass for each species, to test if centroid size could be considered a good proxy for total body size. As no specific mass was known for the individual studied skeletons, we relied on species averages from the literature (see Table 2 for references).

To visually assess the confounding effects of a phylogenetic signal, mean shapes per species were calculated and projected on the PCA scatterplots derived from the variance–covariance matrix of the habitat/locomotion group means. Minimal spanning trees were calculated to estimate the minimal total lengths connecting all data points as an aid in grouping together taxa (Hammer et al., 2001). Furthermore, we made use of several quantitative approaches to evaluate the influence of phylogeny on patterns of morphological variation in the dataset. We first used Caper 1.0 in R (Orme et al., 2018) to conduct phylogenetic generalized least squares (PGLS) regressions (Martins & Hansen, 1997) and regressed the relevant PC scores on our functional/habitat groups (Barr, 2014a, 2014b; Curran, 2015; Meloro, 2007, 2008; Walmsley et al., 2012). In this weighed regression, phylogeny was incorporated as an error term during regression of the shape variables on ecological/locomotor categories transformed into dummy variables (Martins & Hansen, 1997; Walmsley et al., 2012). Using tree branch lengths to estimate phylogenetic covariance, the covariance for two given species was proportional to the sum of branch lengths from the root to the last common ancestor (Monteiro, 2013). It was assumed that cervid traits evolved according to a simple Brownian motion model (see Barr & Scott, 2014; Felsenstein, 1985; Monteiro, 2013). Habitat/locomotion groups were transformed into dummy variables and regressed against the mean shape coordinates for each species (Walmsley et al., 2012). Phylogentic tree distances were downloaded from the 10kTrees website (Arnold et al., 2010) and topologies checked for accuracy against Heckeberg's (2020) cervid phylogenetic framework. Two species present in the main analyses (Cervus mariannus and Cervus alfredi) were excluded from the PGLS regressions, as they were absent from the phylogenetic tree. Cervus elaphus and Cervus canadensis were considered conspecific in the phylogenetic tree and consequently lumped together in the PGLS regressions. Pagel's λ was used as a measure of phylogenetic dependence (Pagel, 1999). In this approach, a maximum likelihood estimate is used to find the value that best explains variation between species at the tips of a phylogeny (Edwards & Cavalli-Sforza, 1964; Kamilar & Cooper, 2013). More specifically, this measure was derived by multiplying all off-diagonal elements (or the covariances between species pairs in the phylogenetic variance–covariance matrix) by λ (Harmon, 2019; Pagel, 1999). Values close to 0 indicate a weak phylogenetic signal, and values close to 1 indicate that related species were morphologically more similar (Molina-Venegas & Rodriguez, 2017).

In addition, we used the phylogenetic tree from 10kTrees (Arnold et al., 2010) to calculate Blomberg's K-statistic (Blomberg et al., 2003) for the relevant principal components. This statistic gives the magnitude of phylogenetic signal as a ratio of the mean squared error in phylogenetic tip data measured from the phylogenetic corrected mean and the variance as expected under a Brownian motion model (Blomberg et al., 2003; Münkemüller et al., 2012). Values below 1 indicate that related species resemble each other less than expected under a Brownian motion model, while values above 1 suggest that related species have more similar traits (Diniz-Filho et al., 2012). This analysis was conducted in the Picante package for R (Kembel et al., 2010). Statistical significance (p < .05) of K was tested with permutation tests using 9,999 replications (Alvarez et al., 2011).

Although earlier work on cervids (Curran, 2009) already indicated that the effect of sexual dimorphism is small on the shape of the calcaneus, an exploratory analysis was performed on a subset of 31 calcanei of Dama dama (five females and six males), A. axis (five females and four males) and Capreolus capreolus (six females and five males). To examine if separations were confounded by sex differences, a PCA was conducted on the dataset, followed by an npMANOVA on the first four principal components.

4 RESULTS

The results of our repeatability tests revealed a close clustering together of replicates (Figure 3). This indicated that the error due to variation caused by the data acquisition and digitization protocols was minimal. Furthermore, a PCA on 31 D. dama, A. axis, and C. capreolus calcanei of known sex (Figure 4) revealed no clear visual separation between males and females, despite clear interspecific differences. This was confirmed by an npMANOVA on the first four principal components (F = 4.68, p < .001), with pairwise comparisons indicating significant differences between the three species (all p < .001), but nonsignificant differences between males and females (A. axis: p = .786, D. dama: p = .694, C. capreolus: p = .585).

The results of a PCA on all calcanei gave visual separation along the axes of the first four—and especially the first two—principal components (Figure 5). Based on a scree plot of the eigenvalues, PC1 to PC4 (summarizing 99.8% of the total variance) were retained for further analysis. The npMANOVA on the relevant PC scores indicated highly significant overall differences between the habitat/locomotor groups (F = 4.35, p = .0001). Pairwise comparisons (Table 3) showed that all but one relationship exhibited significant differences. The exception was the Type C group, that did not differ significantly from the type B group (p = .2276), suggesting the intermediate specimens were similar in morphology to the Type B specimens, as far as the shape differences described by the first four axes were concerned.

Visual assessment of the thin plate spline transformations associated with PC1—and representing the minimum and maximum scores along that axis—revealed two main shape changes (Figure 6; Table 4). A first shape change was a difference in angle of the tuber calcanei relative to the articular surfaces of the element. Specimens with a more positive score had tuber calcanei that were positioned more perpendicular relative to the anterior part of the bone with the articular surfaces. Specimens with a negative score had tuber calcanei positioned less perpendicular relative to the anterior part of the calcaneus. In addition, the first axis also appeared to summarize variation in the relative length of the tuber calcanei. Specimens with a high positive score tended to have calcanei with a relatively shorter tuber than specimens with a low score.

Shape differences noted along PC2 were mostly expressed in the articular surface that supports the malleolus, which tended to be longer in the anterioposterior direction in specimens with a high score (Figure 6; Table 4). In specimens with a low score, this particular surface was shorter in the anterioposterior direction. The variation in shape summarized by PC3 seemed to be mainly expressed in a difference in inferosuperior length of the posterior talar articular surface and in the distance between the articular surface of the malleolus and the posterior talar articular surface (Figure 6; Table 4). The main variation observed along PC4 was a difference in position of the anterior end of the articular surface supporting the cubonavicular relative to the landmarks on the posterior talar articular surface (Figure 6; Table 4). In specimens with a lower negative score, the anterior end of this facet was placed at a more superior height relative to the most inferolateral point of the posterior talar articular surface.

Although the thin plate spline visualizations suggested that a substantial shape variation was explained by the (first two) principal components, it remained to be evaluated how these patterns were reflected in specimen and group distributions in the scatterplots. Despite substantial overlap, Type A specimens gave a high score on the first axis, while Type B specimens gave a lower score (Figure 5). In light of the shape variation explained by this axis, this implied that the Type A and Type B groups had substantially different morphological characteristics. As expected, the intermediate Type C specimens gave intermediate scores on PC1. In other words, there was a continuum from saltatorial specimens of closed environments, over intermediate specimens, to cursorial specimens of open environments. This continuum was also observed along the second axis. It should, however, be mentioned that on both axes the intermediate Type C group showed more visual overlap with the Type B group than with the Type A group. This was in line with the results of the npMANOVA, which indicated that the differences between the Type C and the Type B group were nonsignificant.

The Type D specimens (high altitude) gave high scores on the first axis, most similar to the Type A group. On the second axis, the mountain group was less well separated from the other groups, but generally gave a higher positive score than the Type B specimens. This indicated that the calcaneus of the Type D species was at least morphologically different from the saltatorial Type B species of closed environment. The Type E group (cursorial/open wetland) was relatively well separated from the other specimens when PC1 was plotted against PC2. On PC1, this group produced predominantly negative scores, similar to the Type C and Type B groups. Somewhat unexpected, it suggested that the specimens were morphologically more similar to the intermediate- and saltatorial (closed habitat) species. On the second axis, however, this group gave a high score more similar to the cursorial Type A group. Specimen distributions along PC3 and PC4 did not show good visual separation between the groups (Figure 5). The morphological variation along these axes did not appear to be expressed in any discernable patters.

To test if the observed morphological variation was confounded by allometric size effects, we first regressed centroid size against log transformed average body mass per species. The results indicated that there was a highly significant correlation (R² = 0.875, p < .001) between the size of the calcaneus and the average body mass of the species. Assuming that cervid body mass is a good indicator of total body size (Curran, 2009), centroid size could be considered a good proxy for body size. Then, we assessed the results of the regressions of shape against log centroid size (Figure 7). When PC1 (R² = 0.0414, p = .0202) and PC2 (R² = 0.0306, p = .0462) were regressed against log centroid size a significant, but weak, correlation was found. The regression of the third component against log centroid size (R² = 0.005, p = .4199) indicated no significant correlation between size and shape. This implied that only a small amount of the shape variance explained by the first three components could be attributed to allometry. In PC4, the allometric effect was potentially somewhat stronger, with a—still relatively weak—correlation with centroid size (R² = 0.0975, p < .001).

A PGLS regression of habitat/locomotion dummy variables on PC1 to PC4 (R² = 0.283, overall p = .236) resulted in a high Pagel's λ (0.906), but one that was not significantly different from 1 (p(H₀:λ = 1) = .2902) or from 0 (p(H₀:λ = 0) = .3227). The λ value was high in PC1 (λ = 0.86, p(H₀:λ = 1) = .0106, p(H₀:λ = 0) = .0217) and PC4 (λ = 0.949, p(H₀:λ = 1) = .632, p(H₀:λ = 0) = .102), but lower in PC3 (λ = 0.683, p(H₀:λ = 1) = .034), p(H₀:λ = 0) = .0217), and especially in PC2 (λ = 0.33, p(H₀:λ = 1) < .001), p(H₀:λ = 0) = .348). These results indicated that there was a phylogenetic signal present in PC1, PC3 and PC4, but as p-values were not significantly different from the upper (1) and lower bound (0) for PC1, and both significantly different from the lower and upper bound for PC3 and PC4, it was implied that phylogeny was only partially responsible for the morphological differences summarized by the PCA. For PC2, the values showed that the phylogenetic signal was limited. Our calculations of Blomberg's K-values for the relevant principal components (PC1–PC4) confirmed this relatively limited phylogenetic signal. PC1 (K = 0.356, p = .032), PC2 (K = 0.407, p = .02), PC3 (K = 0.367, p = .014), and PC4 (K = 0.404, p = .015) gave significant results substantially below 1, indicating traits did not evolve as expected under a Brownian motion model (Blomberg et al., 2003).

To further explore the magnitude of this effect along the first two axes, giving the best separation (PC1 and PC2), we re-assessed the PCA scatterplot after calculating mean shapes per species (Figure 8). In several cases, some closely related taxa clustered together. Several members of the Cervini tribe (C. elaphus, C. canadensis, C. unicolor, Cervus timorensis, and A. axis) plotted relatively close together. On the other hand, Type E members of the Cervini tribe gave lower scores on PC1, similar to other, unrelated type E taxa, such as Hydropotes inermis.

The members of the genus Muntiacus also clustered together. When PC1 was plotted against PC2, Muntiacus muntjak was close to Muntiacus reevesi, but nevertheless closer to some other, unrelated taxa such as C. capreolus and C. mariannus. Moreover, the Elaphodus-members of the tribe gave different scores despite their relatedness to Muntiacus and plotted out with other species of similar habitat/locomotor strategy. A similar pattern was noted for the Rangiferini. Despite being most closely related to the Odocoileini, they were positioned near the main group of Type A species. The Alceini tribe, consisting of Alces alces only, was positioned near the somewhat related Rangiferini. It was, nevertheless, also close to the other cursorial Type A species. The Capreolini, with the genera Capreolus and Hydropotes, gave very different results on the first two axes of the PCA. Capreolus gave slightly negative scores on PC1, close to Axis kuhlii and M. muntjak. Hydropotes gave scores similar to the Type E Cervini.

5 DISCUSSION

Although our dataset did not allow us to fully test the effects of sexual dimorphism for all species, the analysis of a subsample of specimens of known sex suggested that shape-related differences between males and females were probably limited for the calcaneus and did not substantially obscure phylogenetic and functional patterns. Similarly, our regressions of the shape variables on centroid size indicated that allometric size effects played a limited role as a confounding factor.

Overall, it appeared that morphological variation along the first two axes of the PCA was mainly driven by function, and in the case of PC1, perhaps to some extent by phylogeny. Along the first axis, we noted a variation in angle of the tuber calcanei that appeared to be linked to hypothesis 2 that in cursorial forms of open environment the calcaneus is expected to be more vertically placed relative to the hock joint, while a more horizontally placed tuber calcanei was related to saltatorial species of closed environment (Figure 6). Furthermore, the variation in length of the tuber calcanei explained by PC1 could most likely be related to hypothesis 1, that the tuber is more elongated in saltatorial forms than in cursorial species. The third functional hypothesis proposed for the calcaneus appeared to be expressed in the shape differences summarized by PC2. Here, it was theorized that cursorial species have a more developed articular surface for the malleolus than saltatorial species. Overall, the Type A, B, and C groups were in line with a functional explanation. The fact that the group of mountain specimens (Type D) gave results similar to the cursorial Type A group was not unexpected in this light, as mountain cervids are probably less well adapted to such environments than mountain bovids (Geist, 1998)—and according to some authors (Flueck & Flueck, 2017) prefer flat terrain when found at high altitude—they may retain some cursorial traits. It should be noted that in bovid ecomorphological studies, mountain species were also often “misclassified” into different habitat categories and substantially overlapped with other groups (Kovarovic, 2004; Weinand, 2005, 2007). This may indicate that vegetation type, and not altitude, could be primarily driving the morphology of the calcaneus in ruminants.

The position of the Type E group was apparently somewhat contradictory as it was closer to the Type A group along the first axis, but closer to the Type B group on the second axis. This pattern may be related to the fact that the wetlands inhabited by Type E forms, are usually dominated by tall grasses and reeds that provide more cover than the open landscapes inhabited by more truly cursorial cervids (Curran, 2009). Possibly this group retains both cursorial traits as well as a significant amount of saltatorial traits. Overall, the patterns described by the habitat/locomotion groups along the first two axes appeared to indicate the presence of a functional signal.

Although it could not be excluded that some functional signal was contained in PC3 and PC4 (Figure 6), the weak separations between the pre-assigned groups and the lack of a clear connection between the observed shape variation and our hypotheses, would suggest that this signal was too obscured by other confounding effects, to act as a reliable predictor of locomotion or habitat preference.

Phylogeny may also have played a limited role in driving the shape differences along the axes of the PCA. Our high—albeit nonsignificant—Pagel's λ indicated that there was a phylogenetic signal present, perhaps primarily in the traits summarized by the first and fourth axis. Blomberg's K-values, nevertheless, indicated that the phylogenetic signal was low and that the traits did not depend heavily on heritability. The PCA on the species averages (Figure 8) showed that some closely related taxa had the tendency to cluster together. Nevertheless, it seemed clear that most of the major patterns were functionally driven. Many taxa of similar locomotor and ecological affinity clustered together, despite being distantly related. This was exemplified in a number of cases. The fact that Cervus unicolor and C. timorensis—despite their different ecological preferences—exhibited similar scores, might be due to phylogenetic relatedness. However, the observation that the type E members of the Cervini tribe gave negative scores on the first axis, similar to unrelated species of the same locomotor strategy and ecology (e.g., H. inermis) should probably be interpreted as functional. M. muntjak and M. reevesi were not only close to each other in PC space, but also to other unrelated taxa of similar locomotor/ecological affinity (C. capreolus and C. mariannus), which would at least in part imply a functional pattern. This was corroborated by Elaphodus cephalopus—also of the Muntiacini tribe—taking a different position, closer to other species of similar habitat/locomotor strategy. The position of the Rangiferini, being highly cursorial, and plotting close to the other Type A species, could also indicate a functional pattern, especially as the related Odocoileini gave different scores. The position of the Alceini might be partially confounded by a phylogenetic signal, as they are positioned close to the related Rangiferini. Its position was, nevertheless, also close to the Type A species of similar locomotion/ecological preference, suggesting function may also play a role. In the Capreolini tribe, the division between the genera Capreolus and Hydropotes may also be functional, as the members of the former plotted close to the unrelated, but ecologically more similar M. muntjak, and the latter to the ecologically similar Type E Cervini.

Taking these findings into consideration, we found support for the first two hypotheses proposed for the calcaneus (see Table 1). The shape variation along the first component reflected a gradient from saltatorial species of closed environment with long, horizontally positioned tuber calcanei to cursorial species of open environments, with short, vertically positioned tuber calcanei. Using these traits, the model effectively differentiated between the groups. Although the intermediate group had a transitional shape, overlap with the saltatorial specimens was more substantial, and the model did not discriminate as well from this group. Assuming that the Type D deer (mountain) prefer flat terrain at high altitudes (Flueck & Flueck, 2017; Geist, 1998) or are morphologically more driven by vegetation structure, and that Type E deer (wetland) retain some cursorial and saltatorial traits due to the specific vegetation structure of such open wetlands, these two groups corroborated the validity of hypotheses 1 and 2.

The third functional hypothesis (Table 1) was supported by specimen patterns along the second component in the PCA and by the associated deformation grids. Again, a gradient was observed from cursiorial, open habitat species (Type A) with long articular surfaces for the malleolus, to saltatorial, closed habitat species (Type B) with short articular surfaces (Figure 5). The model could effectively use this trait to differentiate between species of open environments and closed environments. Intermediate specimens (Type C) were morphologically more similar to the closed environment group. Similar to the length and orientation of the tuber, the articular surface supporting the malleolus in mountain species, was analogous to that of cursorial species of open environment, further confirming the interpretation that the morphology of the calcaneus in such species is driven primarily by vegetation structure. Unlike the shape of the tuber calcanei, the shape of the articular surface supporting the malleolus in Type E species was similar to that of dry, open environment species (Type A). The reason for this discrepancy was unclear, but it may again reflect the notion that cervids of open wetlands retain a combination of saltatorial and cursorial traits.

As a range of different methods have been reported in the artiodactyl ecomorphological literature (see, e.g., Curran, 2009, 2012, 2015; DeGusta & Vrba, 2003, 2005a, 2005b; Plummer et al., 2008; Schellhorn, 2009; Scott, 2004), it is difficult to compare the results of our models to those of other studies. This is also the case for the npMANOVA approach we used, that, despite its advantages, did not provide a measure of reclassification accuracy, as is sometimes seen in studies that make use of discriminant analyses (Curran, 2009, 2012; Kovarovic & Andrews, 2007). That said, the results were in line with those of earlier GMM-based models where the cervid calcaneus performed well as a habitat predictor (Curran, 2009, 2012, 2015). In contrast, studies of bovid ecomorphology have been less unanimous in their estimation of the calcaneus as a good habitat predictor (Barr, 2018; Kovarovic, 2004; Schellhorn, 2009; Schellhorn & Pfretzschner, 2015). Some researchers considered it a useful habitat predictor, but also warned for the potentially confounding effects of allometry and phylogeny in this element (Barr, 2018; Kovarovic, 2004). Despite these reservations, in our model allometric and phylogenetic effects did not seem to substantially obscure the functional signal. This could suggest a discrepancy between cervids and bovids, but could also be the consequence of methodological differences. In this context, it is worth noting that in Barr's (2018) bovid analyses many of the shape variations measured in the calcaneus were similar to those in the model presented here: variation in the length of the tuber calcanei and variation in the articular surface supporting the astragalus. Where the latter trait was considered mainly functional by Barr (2018), it was purported that the length of the tuber was more confounded by body size and phylogenetic relatedness (Barr, 2018). One possible explanation is that in bovids the larger species are driving the allometric signal. In bovids many larger forms exist (e.g., Syncerus, Bos, Taurotragus) that may require additional morphological accommodations to support their weight (Scott, 1979), and are probably too heavy to support the saltatorial or cursorial adaptations typically seen in smaller forms (Geist, 1998; Plummer et al., 2008). This could obscure the functional signal (Geist, 1998). Although some ecomorphological studies (Scott, 1979; Plummer et al., 2008) have a priori excluded species of very large body size from their analyses for precisely these reasons, most that have focused on the calcaneus, have not done so (Barr, 2018; Kovarovic, 2004; Kovarovic & Andrews, 2007; Schellhorn, 2009; Schellhorn & Pfretzschner, 2015).

Another explanation for the discrepancy between cervid and bovid studies of the calcaneus could be that the GMM methods, used in cervid studies so far (Curran, 2009, 2012, 2015, this study), more efficiently exclude size differences from the dataset than the linear size corrections often used in bovid studies (Barr, 2018; Kovarovic & Andrews, 2007). Although it should be noted that the GPA procedure, used here and in Curran's (2009, 2012, 2015) work, only accounts for size differences in the measured object itself, it does so in a more efficient way than traditional morphometrics (Viscosi & Cardini, 2011). If we accept that the size of the calcaneus is a good indicator for total body size, an assumption that appears to be confirmed by the correlation we found between calcaneus centroid size and total body mass, it is likely that our GMM approach more efficiently accounted for total body size differences as well.

The above explanations may account for the absence of an allometric effect in the cervid calcaneus, but they do perhaps not fully explain that only a limited phylogentic signal was found in our study, in comparison with certain bovid studies (Barr, 2018; Kovarovic, 2004; Kovarovic & Andrews, 2007). Our results suggested that many of the larger patterns observed in our analyses were functionally driven, and that there was only limited phylogenetic effect for the calcaneus in our model. These limited effects can possibly be seen in some related taxa clustering together in the scatterplot of PC1 and PC2 (Figure 8). This was not unexpected, because some cervid tribes have an evolutionary history of adaptation to certain habitats (e.g., muntjacs to closed habitats) (Geist, 1998).

DeGusta and Vrba (2003) argued that phylogenetic effects can be easily excluded from morphometric datasets, by selecting those anatomic features that co-vary with locomotion and habitat. While this is true to some extent, the underlying assumption is that a specific feature is either driven by phylogeny or by function. Although it should be remembered that our phylogenetic analyses were based on a limited number of species, what appeared from our results is that for most shape traits, this is not entirely possible. Features such as the length of the tuber calcanei are not driven by either phylogenetic relatedness or functional aspects, but most likely by a combination of both (Elton et al., 2016). That being said, when considering the major patterns observed between the large and diverse habitat/locomotion groups in our PCA, the effects of phylogeny were probably not strong enough to substantially confound the functional signal. This would suggest that our models are useful at predicting ecological and locomotor affinity in fossil samples.

Nevertheless, our method also had its limitations. First of all, it is important to keep in mind that, despite the significant overall differences observed between cervids of different habitat/locomotor strategy, there was still substantial overlap between the groups. For a robust analysis of a group of fossil specimens, it would therefore be advisable to test this on a sample of sufficient size. The results of such an analysis would also be limited in scope, in the sense that they are only informative about one aspect of the cervids' ecology: the vegetation structure of the habitat in which the animals evade predators. Depending on their behavioral repertoire, some deer can make use of different habitats for feeding, sleeping or reproduction (Geist, 1998). Looking at other parts of the skeleton as well—including the masticatory apparatus for feeding behavior—can lead to a more precise reconstruction of a species' ecology (Bishop et al., 2006). Ecomorphology is also but one approach in paleoecology. The results of ecomorphological studies are best considered in concert with other proxies for a robust paleoenvironmental reconstruction (DeGusta & Vrba, 2003).

Finally, it should be noted that the functional framework on which our predictive models were based relied on the assumption that similar selection pressures (i.e., predator evasion strategy) drove the morphology of the extinct deer species on which our models are ultimately intended to be applied. If we accept the notion that phylogeny played a limited role in driving the traits of the calcaneus we studied, this was a reasonably assumption, as the underlying biomechanics that drive functional morphology are constant (DeGusta & Vrba, 2003). It should, however, be taken into account that in certain insular ecosystems with an impoverished carnivore fauna, predation pressure can be reduced (Bouteaux, 2005; McClain et al., 2006). The morphology of deer living in such isolated conditions may therefore be driven by different factors than that of species from continental ecosystems.

6 CONCLUSIONS

In this article, we demonstrated that the morphology of the cervid calcaneus can be used to predict the habitat preferences of taxa of unknown ecological affinity. The shape of this bone was found to be primarily a good predictor of vegetation structure. In addition, we improved the way ecomorphological studies are conducted, by using GMM and by relying more heavily on the underlying functional aspects driving shape. Applying the methods presented here could contribute to a better understanding of Pleistocene environments.

AUTHOR CONTRIBUTIONS

Ben Gruwier: Conceptualization (lead); data curation (lead); formal analysis (lead); funding acquisition (equal); investigation (lead); methodology (lead); resources (equal); visualization (lead). Kris Kovarovic: Conceptualization (supporting); funding acquisition (equal); methodology (supporting); supervision (lead).