Geography of artiodactyl locomotor morphology as an environmental predictor

We investigate locomotor function in artiodactyls, represented by calcaneal gear ratio, as it relates to multiple environments. Using an ecometric approach, we develop a trait–environment model to investigate ecosystem‐level changes through time and to reconstruct past environments. We apply the trait–environment model to a case study of six sites in Kenya to evaluate changes over the past 100 years.


| INTRODUC TI ON
Functional traits are measurable features that influence an organism's interaction with its environment (McGill et al., 2006;Violle et al., 2007), and, because of this relationship, the morphological composition of a community will change as environment changes (Enquist et al., 2015;Violle et al., 2007). Study of these functional trait-environment relationships over spatial and temporal scales and at the community level has come to be called "ecometrics" (Eronen, Polly, et al., 2010;Polly et al., 2011;Vermillion et al., 2018). Ecometric indices have been used to reconstruct palaeoenvironment (Fortelius et al., 2002, evaluate extinction risk (Polly & Sarwar, 2014), understand the impact of non-ecological processes on patterns of biodiversity  and estimate community vulnerability to environmental change (Barnosky et al., 2017;Polly & Head, 2015).
In ecometric studies, the focus on morphological change makes it possible to examine the diversity of associations between climates and biotas through time, and the application of ecometric methods has been used to expand conservation palaeobiology (Eronen, Polly, et al., 2010;Lawing et al., 2012;Polly & Head, 2015).
Communities dominated by fauna with high-crowned teeth occur in drier, more open, grassland-like environments and communities dominated by fauna with low-crowned teeth occur in wetter, more closed, forest-like environments (Eronen, Puolamäki, et al., 2010;Fortelius et al., 2002;Janis et al., 2000). Higher-crowned teeth enable animals to tolerate high levels of environmental grit associated with arid environments (Damuth & Janis, 2011;Jardine et al., 2012;Lucas et al., 2014;Semprebon et al., 2019) as well as increased roughage in a diet of primarily grass (Erickson, 2014;Merceron et al., 2016;Strömberg, 2002Strömberg, , 2006. In ecometric analyses of post-cranial skeletal elements, some traits are functionally related to an animal's locomotion strategy, which is directly linked to the landscape in which the animal moves and, thus, the environment in which they live. Mammalian hind limbs are the primary contributors to forward propulsion during quadrupedal locomotion (Dutto et al., 2006;Kilbourne & Hoffman, 2013), and the calcaneum, or heel bone, is the primary lever and insertion point for the gastrocnemius muscle (Radinsky, 1987). Calcanea are excellent elements to use for palaeobiological studies because they have strong trait-environment relationships and are likely to be preserved in the fossil record because they are relatively dense elements without much soft tissue (Hill, 1996). Previous research has correlated calcaneal morphology and locomotor strategy, which, when averaged across a community, can be indicative of habitat (Panciroli et al., 2017;Polly & MacLeod, 2008). The calcaneal gear ratio (i.e. a measure of the overall calcaneal length relative to the length of the calcaneal tuber) is related to ecoregion province, vegetation cover and mean annual temperature in carnivoran communities (Polly, 2010), but this work has not yet been applied to other taxonomic orders.
Whereas carnivorans are the primary carnivores in most ecosystems, artiodactyls are the primary large herbivores. Artiodactyls are typically cursorial with limbs that are restricted to parasagittal movement by hinge joints (Foss & Prothero, 2007;Hildebrand & Goslow, 2001), and, although they are cursorial, there is still a reasonable degree of variation in their locomotor strategies, producing functional trait variation in their tarsals and metatarsals that is related to vegetation cover (Barr, 2017(Barr, , 2020. These circumstances render an investigation of artiodactyl calcanea an important next step in building ecometric models to relate morphology to environment. Previous research on the ecomorphology of artiodactyl postcrania has focused primarily within the Bovidae and at the species level. Trait-environment relationships have been found between habitat and morphology of the bovid astragalus (Barr, 2014(Barr, , 2017Plummer et al., 2008), metatarsal (Barr, 2017), phalanges (DeGusta & Vrba, 2005), calcaneum (Barr, 2020) and femur (Kappelman, 1988), as well as between precipitation and morphology of the astragalus and metatarsal (Barr, 2017). Similarly, Curran (2012) connected habitat type to the morphology of hindlimb elements in Cervidae. Recently, Dunn and Avery (2020) expanded this work to include five families of artiodactyls in their analysis of 3D calcaneal morphology and found similarities across taxonomic families supporting the use of calcanea in studies of ecomorphology. However, artiodactyl post-crania have not been evaluated at the community level in these studies.
Here, we test whether the trait-environment relationship between calcaneal gear ratio and ecoregion that was documented in Carnivora by Polly (2010) is also present in Artiodactyla, and we explore other environmental variables that may also produce good ecometric models. We contribute and analyse newly collected data of functional trait measurements on the calcanea of artiodactyls.
Functional traits are analysed for relationships to richness, taxonomic family and environment, including precipitation, temperature, elevation, vegetation cover and ecoregion division. All five environmental variables are expected to have a relationship with artiodactyl calcaneal gear ratio, because the environmental variables are closely related to vegetation, substrate and topography. We expect that ecoregion division will have the strongest relationship with the gear ratios of artiodactyl communities because ecoregion divisions are designated by climatic conditions, physiography and vegetation physiognomy, which all contribute to the environmental conditions in which an animal locomotes (Bailey, 2005). Furthermore, as a case study, we demonstrate the application of community-level artiodactyl calcaneal gear ratios at six sites in Kenya documented in Tóth et al. (2014a).

| ME THODS
To investigate calcaneal gear ratio as an ecometric trait in artiodactyls, we collected gear ratio measurement from skeletal specimens preserved in museum collections and cross-referenced the spatial distributions of species to determine community composition. The variation in gear ratio, summarized at the community level, was evaluated for fit with five environmental variables. All analyses were performed in the R statistical computing environment (R Core Team, 2016), and maps were made in ArcGIS 10.2.2 for Desktop (ESRI, 2017).

| Study system
With 240 species in 89 genera in 10 families (Wilson & Reeder, 2005), Artiodactyla has a nearly global distribution in almost all ecosystems, and species frequently overlap geographically creating a myriad of unique communities. Calcanea of artiodactyl individuals (n = 572) from 157 (65%) extant species were examined from natural history collections (see Acknowledgements; Data S1.1). We resolved discrepancies in taxonomy and followed Wilson and Reeder (2005).
Only adults were measured and, when possible, both male and female specimens were included to account for sexual dimorphism. Domestic species were removed from analysis following Gentry et al. (2004).
Species were retained in the analysis if the species are extant and native or reintroduced. Zoo specimens were kept in the dataset because, in artiodactyls, the captive environment has been shown to not affect post-cranial morphology (Bormet & Lawing, 2012;Bormet & Polly, 2015). Number of individuals per species ranges from one to 13, and the proportion of males and females varies with many recorded as "unknown." Measurements of individuals were used to calculate mean and standard deviation values for each species. Mean values were evaluated for differences at the family level.
Calcaneal gear ratio is a measurement of the overall length of a calcaneum divided by the length of its tuber, or in-lever (Polly, 2010; Figure 1). These two linear measurements were measured with callipers from the medial tubercle of the calcaneal tuber to the furthest extent of the cuboid facet and the furthest extent of the dorsal edge of the astragalar facet, respectively. The differing orientation of the sustentacular process means the methodology used by Polly (2010) to measure the length of the calcaneal tuber in carnivorans must be modified for use on artiodactyls. The ratio of these measurements indicates the position of the sustentaculum, which articulates with the astragalus to form the ankle joint. In carnivorans, a low gear ratio indicates a long in-lever and a more plantigrade stance whereas a high gear ratio indicates a short in-lever and a more unguligrade stance (Polly, 2010). Here, we explore this relationship, which has not yet been determined in artiodactyls.

| Sampling strategy and environment
Species' geographic range maps were downloaded from IUCN (2018).
Communities were sampled from overlapping range maps using a global 50 km equidistant points (Polly, 2010 available at https://polly lab.india na.edu/data/index.html). Evenly spaced points at this scale are representative of geographic mixing and the derived patterns are appropriate for comparing with samples from the fossil record (Polly, 2010). Where Hurlbert and Jetz (2007) suggest to only use cells of at least 200 km per side when evaluating global diversity trends, Lawing et al. (2017) showed that, in mammals, ecometric relationships do not change when evaluating records at 50, 100 or 250 km distance. Thus, we chose to use 50 km points to be consistent with previous ecometric studies. In total, there were 54,090 terrestrial global points after removing Antarctica from the dataset, with 47,420 being occupied by at least one artiodactyl species. The remaining points were dropped from the dataset because they are not within the geographic distribution of artiodactyls. For each community with three or more species with trait data available, we calculated mean and standard deviation of calcaneal gear ratio from species' means (n = 21,827). cultivation was excluded from this analysis (Matthews, 1983(Matthews, , 1984(Matthews, , 1999; Figure 3a). Bailey's ecoregions define large-scale areas with similar ecosystem features, such as climate, physiography and vegetation (Bailey, 1989(Bailey, , 1998(Bailey, , 2005Bailey & Hogg, 1986;Figure 3b

| Ecometric summaries and analyses
Community means of calcaneal gear ratio were evaluated for relationships with species richness, and each of the five environmental variables. Continuous variables were transformed for normality as follows: cube root of richness and elevation, cube of mean annual temperature and log of annual precipitation. Continuous variables were analysed with a Pearson's correlation, and categorical variables were analysed with an ANOVA-derived R 2 . We conducted a sensitivity analyses to evaluate whether community means of calcaneal gear ratio change when only a subset of species are present in a community. We performed 10 iterations of randomly sampling species within communities to 95%, 90%, 80%, 70%, 60% and 50% and calculated the average community gear ratio for each sample.
We compared the difference between the gear ratio values of each sample with gear ratio values for each original community.
For each environmental variable that explained a large amount of gear ratio variation, an individual ecometric space was constructed following Lawing et al. (2012) and Vermillion et al. (2018). Gear ratio values were binned into 25 × 25 cells by mean and standard deviation within an ecometric space. We used maximum likelihood to estimate the most likely environmental conditions for each trait bin.
We fitted a Gaussian distribution to each bin individually to create a likelihood surface. Then, we extracted the maximum point on the likelihood surface to determine the most likely environmental value.
We used the ecometric models to estimate environments for all communities in the dataset. For model verification, we repeated five iterations of randomly partitioning the dataset into training and testing data. We used 80% of the communities to build the ecometric spaces and estimated the environment for the remaining 20% of the communities (Table S1.1). With these environmental estimates, the models can be evaluated by comparing the observed environment to the estimated environment. Differences between the observed and estimated values were used to produce anomaly maps (Vermillion et al., 2018). Continuous variables were evaluated using numeric anomaly values. Categorical variables do not allow for numeric anomalies, so these were evaluated using the percentage of estimates that were "correct" (i.e. matching the observed) or "incorrect" (i.e. not matching the observed).

| Conservation palaeobiology
We use six sites in Kenya that are protected areas with historical  and modern ) data on mammal composition (Tóth et al., 2014a). Sites are from a variety of habitats-Kakamega Forest Reserve (forest), Maasai Mara National Reserve (grassland), F I G U R E 2 Environmental variables used to test for trait-environment relationships. (a) Mean annual temperature displayed in 10 colour bins using Jenks Natural Breaks. Dark blue is the lowest temperature (−26.9 to −17.2°C), and dark red is the highest temperature (23.9-31.4°C); (b) Annual precipitation displayed in 10 colour bins using Jenks Natural Breaks. Bright yellow is the lowest precipitation (0-220 mm), and dark green is the highest precipitation (4,388-9,916 mm); (c), Elevation displayed in 10 colour bins using Jenks Natural Breaks. Dark green is the lowest elevation (−1,216-214 m) and white is the highest elevation (4,245-6,231 m)

| RE SULTS
Species-level calcaneal gear ratios range from 1.40 (Ammodorcas   lower mean. There is more taxonomic diversity in the high gear ratios with the five highest values from five different families, including Tragulidae, Suidae, Giraffidae, Camelidae and Hippopotamidae, that have many fewer, specialized species with much more restricted geographic ranges and higher gear ratios than the bovids and cervids ( Figure 4b). We did not find a relationship between calcaneal gear ratio and log body mass (R = .077, R 2 = −.00076, p = .3482).
Mean gear ratios are highest in the tropical regions and lowest in the mid-latitudes (Figure 5a). Although community mean gear ratios are

| Trait-environment relationships
All five environmental variables are related to community mean calcaneal gear ratios (Table 1), although mean annual temperature (36.4%) and elevation (21.9%) explain little relative variation.
Ecoregion division is the most strongly correlated with gear ratio and explains 68.6% of the trait variance. Vegetation cover and annual precipitation are also strongly correlated with gear ratio and explain 63.5% and 60.7% of the variance, respectively. Because of their strong relationships with community mean calcaneal gear ratio, these three environmental variables are investigated further to produce ecometric spaces. When using the models to estimate environment, a sensitivity analysis determined there are only minimal differences between using all of the communities and dividing the communities into training and testing datasets (Table S1.1). Results for all of the communities are presented here.
In the ecometric space of ecoregions, humid tropical savannas make up the largest division that is nearly centred in the ecometric space (Figure 6a) The ecometric space for vegetation cover is dominated by evergreen and grassland classes (Figure 6c). Both evergreen and grassland classes range from low mean and low standard deviation to high mean and high standard deviation, but communities in evergreen classes generally have higher gear ratio means than those in grassland classes. Communities in desert classes follow the same lowlow to high-high pattern that evergreen and grassland communities do, but desert communities are consistently at a lower mean than most of the communities in evergreen or grassland communities.
Communities in deciduous classes have low-to-moderate means and standard deviations. The only trait bin of arctic communities has a moderate mean and low standard deviation. In general, as density of vegetation cover increases, the mean gear ratio also increases.
When the ecometric model is used to estimate vegetation cover based on community calcaneal gear ratio, the model is accurate for only 49.1% of the communities (Figure 6d;

| Conservation palaeobiology
Artiodactyl trait shifts support the homogenization of communities in Kenya (Figure 7). Nairobi National Park did not shift trait bins from historical to modern, but Maasai Mara, Lake Naivasha, and Tsavo East and West converged towards the trait bin in which Nairobi occurred.
At Maasai Mara, the mean calcaneal gear ratio increased, and at Lake Naivasha and Tsavo East and West, the mean and standard de-

Mara had the largest increase in precipitation while Tsavo East and
West had the largest decrease in precipitation (Figure 7c; Table S1.4).
Because trait composition at Nairobi National Park did not change, historical and modern environmental conditions did not shift either.

| D ISCUSS I ON
Here, we demonstrate the ecometric value of artiodactyl calcanea by quantifying the trait-environment relationship between mean gear ratio and five environmental variables, three of which make reasonably good ecometric models. Functionally, the calcaneum is directly related to an animal's locomotion strategy and, consequently, the environment in which the animal lives. In artiodactyls, there are strong relationships between calcaneal morphology and taxonomic family ( Figure 4b); species richness ( Figure 5b); and ecoregion division, vegetation cover and precipitation ( Figure 6, Table 1).

| Community trait composition
Communities with high species richness have moderate gear ratio values whereas communities with low species richness have gear ratios across the range ( Figure 5). This is particularly evident in eastern and south-eastern Africa, which is where all communities with species richness greater than 18 occur. Although the richness is high in these communities, the mean gear ratio when richness is 18 or higher is 1.50, equivalent to the global average. As an example of the differences in functional diversity in these communities with high richness, stotting bovids have longer calcanea than non-stotting bovids in the same habitats, possibly to facilitate upward propulsion and predator avoidance (Barr, 2020). As a result, the stotting bovids contribute lower gear ratios than non-stotting bovids to the community trait value. In a sensitivity analysis, we demonstrated that, even with the diversity that comes with richness, the community trait value is not sensitive to changes in richness ( Figure S1.1).
Conversely, in areas of low richness, such as North and South America, mean gear ratio occurs across the range of values, depending on the environments in which they occur. With low richness, species that have a more specialized locomotor strategy related to the environment may shift the mean gear ratio higher or lower than the mean. Community trait composition is expected to be related to taxonomic designations because of the constrained morphology that is typical within clades Polly, 2010;Polly et al., 2017). Often ecometric analyses are conducted for communities of species within a particular loosely defined guild (e.g. carnivorans or artiodactyls), because these species function in a similar way in terms of usage of resources and exploitation of similar environments and therefore may have a similar morphological response. However, this relationship does not typically mask the overall trait-environment signal at the community level. In an analysis of North American mammal communities, accounting for phylogeny did not significantly change the trait-environment relationships .
Notably, carnivorans and artiodactyls have nearly the same range of variation in calcaneal gear ratios but are functionally reversed during locomotion. Globally, carnivoran gear ratios range from 1.08 in Melursus ursinus (sloth bear; plantigrade) to 1.46 in Ichneumia albicauda (white-tailed mongoose; digitigrade; Polly, 2010). Artiodactyls are unguligrade, and, in the data presented here, species-level mean gear ratios range from 1.40 in Ammodorcas clarkei (dibatag) to 1.74 in Hexaprotodon liberiensis (pygmy hippopotamus; Figure 4a). The highest values of carnivoran gear ratios overlap the lowest values of artiodactyl gear ratios and this overlap represents the highest efficiency in terms of speed in both groups (Figure 8). In carnivorans, a longer calcaneal tuber (i.e. smaller gear ratio) provides more power and less speed, whereas a shorter calcaneal tuber (i.e. larger gear ratio) provides less power and more speed (Carrano, 1997). This is reversed in artiodactyls, possibly due to specialized adaptations, such as a double pulley astragalus and the reduction and lengthening of distal skeletal elements. As vegetation cover becomes denser and precipitation increases, we see community-level carnivoran gear ratios decrease and artiodactyl gear ratios increase. In carnivorans, calcaneal morphology is related to predation style (Panciroli et al., 2017), so it is notable that artiodactyls have maximized the same outcome (i.e. more power/less speed or less power/more speed) in the same habitats but with different morphology.
Community-level calcaneal gear ratios of artiodactyls and carnivores respond to different suites of environmental variables.
Ecoregion division (68.6%), vegetation cover (63.5%) and precipitation (60.7%) explained the most variance in artiodactyl calcaneal gear ratio (Table 1). Similarly, Polly (2010) found that ecoregion province (70%) and vegetation cover (49%) are strongly related to calcaneal gear ratio of carnivores. But, rather than precipitation as with artiodactyls, temperature (48%) explained a large amount variance in North American carnivoran gear ratios. Because ecoregion and vegetation cover explain the most trait variance in both groups, it is likely that temperature and precipitation are not directly contributing to trait variance but are part of an interaction that produces the associated vegetation cover and ecoregion designation (Polly, 2010  richness and decreasing beta-diversity in the mammal communities. These authors suggest the homogenization pattern was driven by increased presence of common species across the region. We note functional trait composition, captured in the ecometric spaces ( Figure 7), was also homogenized. Over the same time, these sites saw a homogenization of habitat because of increased human land use (Tóth et al., 2014b) as well increased annual precipitation in the west and decreased annual precipitation in the east (Gebrechorkos et al., 2019). Trait values indicate a homogenization of vegetation cover towards grassland-based categories, which is consistent with human land use patterns (Tóth et al., 2014b). Furthermore, the trait-based predictions of precipitation change are also consistent with regional patterns documented across Kenya (Gebrechorkos et al., 2019). In the west, Kakamega indicates an increase in precipitation (i.e. increase in mean) and, in the east, Tsavo East and West indicate a decrease in precipitation (i.e. decrease in mean). showed that only 25% of species within a community are needed to predict accurate ecometric relationships (Polly & Sarwar, 2014).

| Limitations
As geographic ranges and trait databases become more commonplace, increased availability of data will increase the fit of traitenvironment relationships.
No-analog communities include species that are extant today, yet they do not occur in the same combinations today as they did in the past (Williams & Jackson, 2007). No-analog communities are expected to increase as species reshuffling occurs in response to climatic and environmental changes (Hobbs et al., 2018;Williams & Jackson, 2007). Ecometrics provides a method for comparing noanalog communities with trait compositions within the modern trait bins of the ecometric spaces (Vermillion et al., 2018). Thus, maximum likelihood predictions are possible for no-analog communities, but this method is limited for communities with trait compositions that do not exist in the trait bins of the constructed ecometric spaces.
Maximum likelihood also produces an ecometric model that underpredicts precipitation in wet areas and overpredicts precipitation in dry areas. This tendency towards the mean is expected because maximum likelihood uses the most frequent precipitation of the communities in a trait bin to predict precipitation for all communities in that trait bin. However, in an analysis of different methods used for estimating ecometric trait-environment relationships, maximum likelihood consistently produced the most accurate estimates . Additionally, the severely over-and underestimated areas may also be highlighting areas of environmental change where the faunal community has not yet responded. Many of the areas that are overestimated are also drylands that are or will be vulnerable to increased aridity by 2100 (Berdugo et al., 2020). In these areas, the fauna may be lagging aridification, and thus, the morphology of native, extant communities conveys past environment.
Community trait composition is likely to be impacted by domestics and other non-natives filling new ecological roles or replacing native fauna (Smith et al., 2016). Many domesticated and livestock animals are artiodactyls, and we have removed these species from the dataset used here following Gentry et al. (2004). We expect domestic species that are heavily managed may weaken the ecometric relationships because the species are not naturally interacting with their environment but are dependent on human activities. However, non-native species that are well suited to their new environment and are not dependent on humans for survival may strengthen the ecometric relationships. As scientists consider agricultural land as a distinct system in the middle of a continuum ranging from unmanaged wilderness to human-dominated urban areas (Swinton et al., 2007), the impact of these species on community composition should be further investigated.

| CON CLUS IONS
Here, we present a new application of ecometrics by demonstrating a trait-environment relationship between calcaneal gear ratio and ecoregion division (68.6%), vegetation cover (63.5%) and precipitation (60.7%). These methods have been used before but have not been applied to artiodactyl post-crania or applied at a global scale.
Artiodactyls make up a large proportion of the world's megabiota, which is disproportionately threatened, yet critical for ecosystem health and resiliency (Enquist et al., 2020). Our ecometric models presented here have the potential to be tools in understanding biodiversity responses to environmental change and could be applied to the palaeontological and historical record. For instance, in the Kenyan sites, detailed historical environmental data are not available, but with the trait-environment relationships, it becomes possible to track environmental changes using calcaneal gear ratios.
As habitats changed, the communities of the protected areas were changed by the movement of common species between the protected areas (Tóth et al., 2014b). With ecometric models, it will be possible to detect community responses to environmental changes in a way that informs conservation strategies.

ACK N OWLED G EM ENTS
We thank E. Panciroli and an anonymous reviewer for their comments, P.D. Polly and P. Barboza for comments on a previous version of this manuscript and S. Holte, C. Janis and J. Martin for discussions that improved this manuscript. J. Martin helped with data collection. K. Brown worked as an undergraduate intern on this project and contributed to data curation. A.-C. Fabre graciously hosted the lead author's visit to London. We are grateful to the museums and collections managers who supported our data collection for this