Causal explanations for the evolution of ‘low gear’ locomotion in insular ruminants

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Journal of Biogeography published by John Wiley & Sons Ltd 1German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Synthesis Centre for Biodiversity Sciences (sDiv), Leipzig, Germany 2Museum für Naturkunde, Leibniz-Institut für Evolutionsund Biodiversitätsforschung, Berlin, Germany 3ARAID Foundation, Instituto Universitario de Investigación en Ciencias Ambientales (IUCA)-Grupo Aragosaurus, Universidad de Zaragoza, Zaragoza, Spain 4Department of Rangeland, Wildlife, and Fisheries Management, Texas A&M University, College Station, TX, USA 5Don Sundquist Center of Excellence in Paleontology, East Tennessee State University, Johnson City, TN, USA

Bone fusions likely dampen instability of the joints and reduce lateral movement of the stride by restricting functionally undesirable lateral and planar movement (e.g. Bover & Fornós, 2005;Van der Geer et al., 2011). However, they are rare in insular ruminants, except for the fusion of the metatarsal with the naviculo-cuboid and cuneiforms in Myotragus and Hoplitomeryx (see e.g. Bover & Fornós, 2005;Van der Geer, 2014a). Myotragus balearicus is also characterized by shortened phalanges, that is, shorter feet, rare fusion of the distal fibula (os malleolare) and tibia (observed also in Hoplitomeryx), and occasional carpal fusions (scapholunar and pisiform-cuneiform fusions), making this species the most extreme case of 'low gear' locomotion (see e.g. Bover & Fornós, 2005;Sondaar, 1977;Van der Geer, 2014a).
A few studies presented an overview of taxa exhibiting morphological features associated with 'low gear' locomotion, discussed problems in assessing its generality and suggested that the main driver behind this evolutionary pathway would be the lack of mammalian predators on the focal islands Sondaar, 1977;Van der Geer, 2005, Van der Geer, 2014Van der Geer et al., 2011). Nevertheless, a quantitative investigation of the causal forces influencing the acquisition of this peculiar type of locomotion in insular ruminants has yet to be conducted. Bovidae is a highly diverse clade of large mammals (e.g. Bibi, 2013;Bibi et al., 2009), which provides an excellent opportunity to study how locomotion on different terrains affects loading regimes experienced by distal limb elements and, consequently, their morphology (Higgins, 2014). Here we focus on insular bovids, which encompass a variety of species exhibiting different degrees of metapodial shortening and thickening (see also  and which inhabited or are still living on islands located in different regions and characterized by varied palaeogeographic histories and evolutionary pressures. We expand the dataset for insular bovids as given in Rozzi and Palombo (2014) and investigate the selective factors hypothesized to influence the evolution of 'low gear' locomotion extensively for the first time. We concentrate on the two main morphological traits associated with this peculiar type of locomotion (i.e. short and robust metapodials) and we predict that they should be most pronounced for: 1. Smaller species and/or species that are strongly reduced in size relative to their mainland ancestors or closest mainland relatives.
2. Species that inhabited or are still living on islands with the fewest predators ( Figure 1).
3. Species that inhabited or are still living on islands characterized by homogenously mountainous landscapes ( Figure 1).

| Dataset
The dataset includes the majority of extant (6) and Plio-Holocene (15) insular bovids (Table S1). Linear measurements (recorded to the tenth of a millimetre) of distal limb elements were used to estimate shortening and degree of robustness (see below and Tables S1 and S2). Additional specimens were used for morphological comparisons. The sources for these data and a complete list of specimens (780), stored in 13 European and American museums, are reported in Table S3. In three cases (Myotragus balearicus, Ovibos moschatus wardi, and Nesogoral cenisae), estimates were based on additional data given in publications (Table S1). In addition to these measurements, we recorded 10 predictor variables (see below) describing topographic and ecological conditions of the focal taxa or islands (see Table S4, for characterization of 11 focal islands).

| Description of main variables associated with 'low gear' locomotion
We calculated mean robustness indices (RI; transverse diameter at midshaft/total length) of metacarpals (RI Mc) and metatarsals (RI Mt) (Tables S1 and S2, and a few examples in Figure 2). Because stout metapodials are not necessarily short and vice versa, we also defined species average shortening index (SI) of metacarpals (SI Mc; length of metacarpal/length of radius) and metatarsals (SI Mt; length of metatarsal/length of tibia) (see Tables S1 and S2, and a few examples in Figure 3). Evaluating changes in relative proportions of limb elements of fossil insular bovids is extremely difficult because remains belonging to a single individual are rare. Accordingly, we measured bones with matching epiphyseal articulation of comparable size recovered from the same site. Furthermore, because living insular bovids do not exhibit significant sexual size dimorphism (see e.g. Jass & Mead, 2004;Martin & Barboza, 2020;Rozzi, 2017Rozzi, , 2018 and references in those papers), robustness and shortening indices based on specimens of male and female individuals were averaged together to obtain reference mean values.
We performed a correlation analysis of these morphometric variables in R ver. 3.5.0 (R Core Team, 2013) using the function ggpairs in the package 'GGally' (Schloerke et al., 2018). Hindlimb and forelimb shortening and robustness indices are positively correlated (r = 0.633 and 0.953; see Figure S1). Therefore, we selected as our response variables the two indices associated with the highest numbers of observations in the dataset (SI Mt and RI Mc).

| Description of predictor variables: body size and S i , ecology, topography
To test whether the most common morphological features associated with 'low gear' locomotion (i.e. shortening and thickening of metapodials) are dependent on body size and body size change (prediction 1), we took body mass values from the literature (Rozzi, 2018). We also included estimates of S i (= insular body size divergence, sensu Lomolino, 1985Lomolino, , 2005, that is the mass of the insular taxon divided by that of its putative mainland relative, the latter based on geographic proximity and taxonomic designation (see Table S1).
Variables describing ecological conditions of the focal insular communities (see prediction 2) included species richness of large mammals and number of predators or competitors likely to directly interact with the focal insular taxa. Values were taken from Rozzi (2018) and were estimated by first developing a list of all other mammals co-occurring on each island, and then consulting references on the diet and habitats of those species to assess which ones were likely to be predators or competitors of the focal taxa (see Materials and Methods and Supporting Information in Rozzi, 2018 for details). To examine the role that habitat selection had on each island either in driving the acquisition of 'low gear' locomotion or in maintaining/increasing a cursorial aptitude (prediction 3), we obtained island (palaeo) area from different published sources (Table S4) and we calculated different topographic indices (see e.g. Yu et al., 2015).

F I G U R E 2 Selection of metatarsals of extant and fossil insular bovids arranged in order of increasing
We took maximum elevation values for each island (= range of elevation sensu Yu et al., 2015) from the UN Island database (Table S4) (Hijmans, 2018). To obtain the perimeter of each island in the present we used a global land-mask shapefile downloaded from https:// www.natur alear thdata.com/downl oads/10m-physi cal-vecto rs/, and selected our islands with qgis (QGIS Development Team, 2020, ver. 2.18.20). To quantify the irregularity of the terrain, we used a roughness index that measures the degree of irregularity of the surface by estimating the largest inter-cell difference of the central pixel and its eight surrounding cells using qgis. We then reclassified the roughness map into a binary map identifying mountainous areas, that is pixels with a roughness index larger than 100, indicating a difference of over 100 m in a circle of a ~2 km ratio (at the equator).

| Data analysis
We fitted multiple linear regression models to investigate the relationship between response variables (SI Mt and RI Mc) and predictor variables describing life history, ecological, and topographic conditions of the focal taxa and islands. We scaled selected predictor variables with a Z-transformation and we verified model assumptions by inspecting diagnostic plots of the residuals. For both dependent variables, we tested the effects of each predictor separately, and we fitted a 'maximal' model with all predictors, three models including variables associated with the focal predictions (see above; body size and S i , ecology, topography), and a null model (intercept-only).
To further assess the most important predictors influencing the evolution of 'low gear' locomotion in the focal taxa, we ran a regression tree analysis on the RI Mc dataset (which encompasses the highest number of observations; n = 19). Although regression trees-as all other machine learning methods-do not assume data independence (Davidson, Hamilton, Boyer, Brown, & Ceballos, 2009;Lomolino et al., 2012;Melo, Rangel, & Diniz-Filho, 2009;Rozzi, 2018;Westoby, Leishman, & Lord, 1995), we decided to include subfamily of the focal taxa as a predictor variable in the analysis to further control for the effects of phylogenetic affinities. The main product of regression tree analysis is a branching tree that describes the contextual relationships between the response variable (here RI Mc) and a subset of predictors (life history, ecological and topographic variables) (see e.g. Durst & Roth, 2012Lomolino et al., 2012;Lyons et al., 2016;Rozzi, 2018;Van der Geer, Lomolino, & Lyras, 2018). We adopted a holdout validation approach (split ratio = 75:25) and built a regression tree using MSE (mean squared error) splitting criterion. To prevent overfitting we pruned the tree by setting a maximum depth of two. In other words, given the small sample size, we focused only on the first split to interpret our results. Furthermore, we assessed the importance of predictors across many different possible trees by performing a random forest analysis. Random forests fit a number of decision trees on various subsamples of the dataset and use averaging to control overfitting and improve predictive accuracy (Prasad, Iverson, & Liaw, 2006). Accordingly, they do not overfit as easily as decision trees and deal well with high-dimensional feature spaces and small sample sizes (Biau & Scornet, 2016;Breiman, 2001). We created a random forest with 100 trees and used MSE (mean squared error) splitting criterion to build the trees. To visualize the most important features of the dataset, we reported feature importance scores for both the regression tree and random forest models. Furthermore, we used the test data to calculate R 2 and MSE and validate the models' performance. All machine learning analyses were performed using Python 3.7.3 and the packages 'scikit-learn' (Pedregosa et al., 2011), 'pandas' (McKinney, 2011, 'seaborn' (Waskom et al., 2020), 'matplotlib' (Caswell et al., 2020), 'numpy' (Oliphant, 2006) and 'graphviz' (Bank, 2019).  (Table S7). A more complex model, with all ecological independent variables, although statistically significant (see Table 1) and including non-collinear predictors (vifs < 4; Hair, Black, Babin, Anderson, & Tatham, 1998), received a higher (i.e. worse) AICc score (AICc ecology = −33.73; Table S7). Results of regression tree analyses and random forest analyses confirm the importance of large mammal richness and number of competitors, but they also suggest that body size and S i might have played a role in influencing the evolution of stout metapodials in insular bovids. In fact, richness of large mammals, body size, number of competitors and  Table 1) and does not vary across subfamilies (intercepts and slopes are not significantly different between Antilopinae and Bovinae; Table S6). However, the effect of competitors within the model is statistically significant per se (p = 0.02). No significant relationships between SI Mt and other predictors were found and none of our models-including the 'maximal' model-received a better AICc score than the null model (Table S7). Nevertheless, model selection based on ΔAICc suggested that number of competitors (AICcWt = 0.18) contributes, albeit marginally, to influence the evolution of short limbs in island bovids (Table S7).

TA B L E 1
Results of regression analyses of the relationships among response and predictor variables in the study. P-values < 0.05 are in bold. For additional statistical details, see Tables S5, S6, S7

| Allometry and phylogenetic conservatism of bauplan
Our assessment of causal explanations influencing the acquisition of 'low gear' locomotion suggests that the evolution of short and robust metapodials in insular ruminants does not parallel the direction and magnitude of their island rule pattern, that is, it cannot be explained by a simple allometric downscaling of the animals (prediction 1; see also . In fact, no significant relationships were found between our response variables-RI Mc and SI Mt-and body size of the island species or S i (Table 1), and multiple models that included these two predictors received worse AICc scores than null models (Table S7) Tables 1, S5, S6 and S7). (d) Report showing the importance of features included in regression tree (on the left) and random forest (on the right) models and how much they contribute to their predictive accuracy. All importance scores are rescaled to have values between 0 and 1 and only predictors with scores higher than 0.05 are included body size among the features contributing to the predictive accuracy of regression tree and random forest models (Figure 4e) might also reflect the significant effect of subfamily in influencing the evolution of robust metapodials in the focal taxa.
Assessing whether insular bovids exhibit a shortening of limb length and an increase of limb bone robustness relative to their ancestors or closest mainland relatives is challenging because of the dearth of available fossil and phylogenetic information (see e.g. Bover et al., 2019;Rozzi & Palombo, 2013b, 2014. However, our results suggest that short and stout metapodials are traits of convergent evolution, in that 'low gear' has evolved multiple times on different islands across phylogenetic clades (see also Bover et al., 2019;Bover et al., 2010;. In fact, although the effect of subfamily on some of the patterns that we observed is significant, a variety of trait combinations characterizes metapodials of insular bovids, with taxa exhibiting different degrees of robustness and shortening within the same clade (see e.g. Figures 2 and 3). Accordingly, unveiling the role of habitat selection and interspecific interactions in driving the evolution of those traits is essential.

| Ecology
The most extreme cases of 'low gear' locomotion evolved in the absence or scarcity of predators (e.g. Myotragus balearicus, Bubalus mindorensis, Bubalus cebuensis; see also ). This has led many authors to hypothesize that predator release on islands would be one of the main factors driving the acquisition of this peculiar suite of traits (e.g. Bover et al., 2010;Caloi & Palombo, 1994;Rozzi & Palombo, 2013b;Sondaar, 1977; Van der Geer et al., 2011, cf. prediction 2 and Figure 1). Results obtained here challenge this view, as we did not find any significant relationships between number of predators and our response variables (Tables 1, S5,  since the Pliocene and evolved a highly specialized locomotion (see above and Bover, 2004;Bover & Fornós, 2005;Bover et al., 2019;Bover et al., 2010;.

| Topography
Topographic complexity is a well-known measure of habitat heterogeneity and available niches and has been linked to species richness in many ecological communities (see e.g. Hutchinson, 1959;Owen, 1990;Yu et al., 2015 and references in those papers). Thus, we expected taxa that inhabited or are still living on islands characterized by relatively complex physiographies to display long and slender metapodials and maintain a less specialized gait (cf. prediction 3 and Figure 1). This hypothesis garnered some indirect support from our assessment of the importance of ecological predictors in the evolution of 'low gear' locomotion (large mammalian richness and the number of competitors in the focal communities likely reflect habitat heterogeneity of the respective islands; see above).
Nevertheless, we did not find significant relationships between RI Mc and SI Mt and proxies for topographic complexity (mean roughness and accumulated roughness).
Animals from rugged-mountainous areas-including humans and ruminants-tend to exhibit relatively robust and short lower limb bones (see Figure 1 and e.g. Higgins, 2014;Rozzi & Palombo, 2013a,b;Scott, 1985 and references in those papers).
However, much variation exists, even within the same family. For instance, a few mountainous Antilopinae, such as the Japanese serow, Capricornis crispus and the Appennine chamois, Rupicapra pyrenaica ornata, exhibit long and slender metapodials (see also Rozzi & Palombo, 2013b, 2014. Moreover, some Bovinae reported as inhabiting flat terrains (Higgins, 2014)-i.e. the muskox, Ovibos moschatus-are actually characterized by short and robust metapodials (Bover et al., 2019;.   Figure 2 and Table S1). In this regard, the evolution of 'low gear' locomotion and the loss of adaptation to mountain climbing in Myotragus would reflect its great body mass in relation to body height (Bover & Fornós, 2005)

| CON CLUS IONS
The evolution of 'low gear' locomotion in insular ruminants does not simply result from phyletic dwarfing and from the absence or scantiness of predators in the focal communities. Competitive release on species-poor islands plays an essential role in prompting adaptations for this peculiar type of gait, such as shortening and thickening of metapodials. Given that predator diversity is the main factor influencing body size evolution of insular bovids (Rozzi, 2018), different traits associated with the island syndrome may evolve in response to varied drivers, even within the same taxonomic group (see also ; Van der Geer, 2014a). Furthermore, island topography is not as relevant as interspecific dynamics in influencing the evolution of the focal morphological traits, although the amount of mountainous terrain occurring on each island seems to significantly affect the evolution of robust metapodials in insular bovids.
Short and stout metapodials are traits of convergent evolution for insular ruminants, in that they evolved multiple times on different islands across Antilopinae and Bovinae. The marked variation in the responses to ecological and topographic predictors that we observed within Bovidae supports the idea that the evolution of 'low gear' locomotion would be the product of a complex interplay of biotic and abiotic factors, and calls for caution in drawing conclusions on this phenomenon on the basis of isolated, albeit significant cases (e.g. Myotragus). All in all, investigations on other large herbivores, especially cervids and graviportal taxa (e.g. elephants, hippopotamuses), might be crucial to provide further insights into the mechanisms triggering the acquisition of such a peculiar locomotion or promoting the maintenance of a more generalistic gait.

ACK N OWLED G EM ENTS
We thank the following people for their fruitful advice and discussion on the evolution of insular ruminants: Mark V. Lomolino

DATA AVA I L A B I L I T Y S TAT E M E N T
The authors declare that the data supporting the findings of this study are available within the paper and its Supporting Information file.