Anatomical Correlates of Cursoriality are Compromised by Body Size and Propensity to Burrow in a Group of Small Mammals (Lagomorpha)

Highly cursorial animals are specialised for fast, sustained running via specific morphological adaptations, notably including changes in limb segment length and mechanical advantage. Members of the order Lagomorpha (hares, rabbits and pikas) vary in cursorial ability; hares are generally highly cursorial, rabbits more frequently saltate, and pikas predominantly trot. Previous investigations of lagomorphs have identified anatomical trends correlated with this ‘cursoriality gradient’, however, the phylogenetic sampling of such investigations has been limited to three American species, namely the American pika (Ochotona princeps), brush rabbit (Sylvilagus bachmani), and black-tailed jackrabbit (Lepus californicus). Here, we expand the phylogenetic sample and body size range by including novel data from Australian samples of the European rabbit (Oryctolagus cuniculus) and European hare (L. europaeus), alongside unpublished data on the Eastern cottontail (S. floridanus). X-ray Computed Tomography and digital landmarking were used to capture proportions within the appendicular skeleton of ~ 40 specimens of each European species. In doubling the number of species studied, we find the previously-identified morphological gradients associated with cursorial behaviour are complicated when evaluated in the larger sample. The relative length and joint velocity of limbs was found to be lower than predicted in European rabbits and hares. Furthermore, we present a novel assessment of morphological integration in the lagomorph appendicular skeleton, finding between-limb covariation patterns that are generally similar to those of other mammals. Broadly, these results suggest cursoriality is only one of many selective forces driving lagomorph skeletal evolution, with variations in body size and fossoriality potentially having measurable impacts.


Introduction
Cursorial is a term widely applied to any animal that is anatomically and physiologically specialised for running relatively fast or for long periods (Álvarez et al., 2013;Carrano, 1999;Stein & Casinos, 1997;Young et al., 2014). Typically, this locomotor style is associated with a consistent set of morphological adaptations that increase the efficiency of high-speed locomotion; for example, elongate limbs that have proximally concentrated mass and hinge-like joints that articulate parallel to the axial skeleton (Coombs, 1978). Such morphology optimises running performance through increasing either stride length or stride rate (Bramble & Carrier, 1983), often countering specific biomechanical and metabolic constraints via similar anatomical adaptations (Bertram & Biewener, 1990;Biewener, 1989Biewener, , 1990Christiansen, 2002;Hildebrand & Hurley, 1985). Therefore, evaluating skeletal morphology is an established way to infer running ability (Howell, 1944;Gambaryan, 1974;Coombs, 1978;Garland & Janis, 1993;Young et al., 2014).

Aims and Hypotheses
Despite the phylogenetic diversity of Lagomorpha, diversity in their locomotion and associated postcranial adaptations have been previously investigated using only three North American species (Camp & Borell, 1937;Young et al., 2014). Thus, the predicted anatomical cursoriality gradient among lagomorph morphotypes has been estimated from a statistically limited sample that may not be indicative of the wider order. Here we use a larger phylogenetic sample to test whether the expected gradated relationship between cursoriality and postcranial morphology holds. Our sample comprises the American pika (Ochotona princeps), brush rabbit (Sylvilagus bachmani), and black-tailed jackrabbit (Lepus californicus)-data collected by Young et al. (2014)-and a second Sylvilagus species (S. floridanus), the European rabbit (Oryctolagus cuniculus) and European hare (L. europaeus). The latter two species were introduced to Australia during European colonisation, where there has been a need for population control following their subsequent rise as pertinent pest species (Stott, 2015;Thompson & King, 1994). Bone length ratios were calculated following Young et al. (2014) to examine the presence of a cursoriality gradient in two functionally-relevant variables; distal limb segment elongation and joint architecture. Three proxies for cursoriality were used, where the six species were ranked in order of increasing running ability according to (1) maximum reported running speed, (2) degree of craniofacial tilt (per Kraatz & Sherratt, 2016), and (3) propensity to burrow. We then examined the degree of morphological integration (sensu Olson & Miller 1958) between and within the limb bones of lagomorphs to understand the role that integration plays in adaptation to cursoriality. As such, two of Young et al.'s (2014) hypotheses were re-examined alongside a third, novel hypothesis: Hypothesis 1 Distal limb segment elongation will be present to a higher degree in hares than rabbits, both of which will exceed the pika.
The phenomenon of distal limb segment elongation as an adaptation to cursorial locomotion is well-reported (Christiansen, 2002;Coombs, 1978;Gregory, 1912), particularly through comparison of metatarsal-to-femur ratios (Gambaryan 1974;Garland & Janis, 1993;Elissamburu & Vizcaíno, 2004;Fostowicz-Frelik, 2007;Young et al., 2014). This adaptation occurs in both limb pairs in quadrupeds, enabling the extension of stride length. By keeping muscle mass concentrated proximally to decrease rotational inertia, a high stride rate is maintained (Gregory, 1912;Williams 2007b), thereby allowing cursorial animals to move greater distances in less time whilst expending less energy (Hildebrand & Hurley, 1985;Marsh et al., 2006). As such, this adaptation is well-suited to the lifestyle of herbivores that use extreme speed and stamina to outrun predators.

Hypothesis 2
The lowest limb joint velocity (measured via anatomical mechanical advantage) will be observed in the pika and hares will tend towards the highest.
Limbs and their joints can be characterised using the same principles that govern the motion of levers and fulcrums (McHenry, 2012;Smith & Savage, 1956), including mechanical advantage (see Supplementary Materials for comprehensive description). Fundamentally, mechanical advantage provides an approximation of whether a particular taxon favours joint force or joint velocity (Arnold et al., 2011;Bergmann & Hare-Drubka, 2015;Stern, 1974). Extending an out-lever relative to its in-lever grants greater velocity, as the far end is displaced further by the same inforce, but the increased distance from the fulcrum and the dynamics of muscle resistance consequently decreases outforce (Biewener, 1989;Camp & Borell, 1937;Christiansen, 2002). This information can be used to infer the functional 1 3 specialisations of the limb. Joints with low out-force but high velocity are speed-oriented (i.e., low mechanical advantage) (Bergmann & Hare-Drubka, 2015) and are associated with the rapid, sustained joint motion of cursoriality as opposed to the high force production necessary for tasks such as digging (Elissamburu & Vizcaíno, 2004). Due to the oppositional nature of this joint adaptation, cursoriality and burrowing are typically viewed as mutually exclusive functions (a biological "trade-off").
Hypothesis 3 More specialised cursors will display less integration between the forelimb and hindlimb, as well as greater integration of functionally analogous bones compared to serially homologous bones.
One of the key anatomical features enabling mammals to optimize cursoriality was the evolution of the tri-segmented, dorsoventrally oriented limb that allowed for the parasagittal (non-sprawling) stride observed in therians (Erickson et al., 2002;Fischer & Blickhan, 2006;Fischer et al., 2002). This limb structure, alongside vertebral flexion, facilitates a functional differentiation of the limbs; the hindlimbs act as the primary propulsion, whereas the forelimbs provide stabilisation and deceleration (Demes et al., 1994;Fostowicz-Frelik, 2007;Williams et al., 2007b). This differentiation forms the basis of asymmetrical gait patterns, which are predominantly found in mammals, common in cursors, and occur when the pairs of feet contact the ground with uneven spacing (Álvarez et al., 2013). The impacts of this differentiation can be examined through morphological integration (Olson & Miller 1958;Young & Hallgrímsson, 2005), which shows the degree of evolutionary covariation between established 'modules' of a biological system (Klingenberg, 2008). Previous studies have identified the advantages and constraints of certain locomotor modules in varying contexts (Bennett & Goswami, 2011;Garland et al., 2017;Martin et al., 2019;Young & Hallgrímsson, 2005), however, the effect of the asymmetrical half-bound on limb integration in placental mammals remains relatively understudied (Schmidt & Fischer, 2009). The half-bound that predominates in lagomorphs involves greater amounts of skeletal strain than other gaits (Gambaryan 1974), thus, it was hypothesised that cursorial lagomorphs may exhibit decreased integration between the forelimbs and hindlimbs compared to generalists and cursors that use other types of gaits (Martín-Serra et al., 2015).
This functional differentiation is facilitated by a limb segment restructure that involves the mobilisation of the scapula and extension of the pes. This causes a functional reorganisation of the limbs where serially homologous bones no longer perform corresponding locomotor functions, leading instead to functionally analogous segments (Schmidt & Fischer, 2009). Serially homologous elements are derived from the duplication of a module and as such are assumed to be tightly integrated ancestrally (Young & Hallgrímsson, 2005), however, we hypothesise that highly cursorial lagomorphs have evolved to favour integration between functionally analogous bones over serial homologs.

Cursoriality Gradients
Two hare species, three rabbit species, and one pika species were included in this study. We build on Young et al.'s (2014) original sample, which included one species of each lagomorph morphotype (Lepus californicus, Sylvilagus bachmani, and Ochotona princeps), adding unpublished data on S. floridanus and new data on O. cuniculus and L. europaeus. Table 1 summarises several basic aspects of the ecology of the six species, with climate zone information corroborated using the climate type maps of Peel et al. (2007). As specimens of O. cuniculus and L. europaeus in this study are from Australia and not their native Europe, it must be noted that an assumption is made that these specimens are representative of their species at a wider level, when it is possible any deviation from the predicted trends shown by these specimens may be an artefact of their isolation from the environment in which they originated. However, the range of L. europaeus in Australia is restricted to a region with a high climactic match for its natural habitat, potentially limiting the impact of any major factors that may affect the morphology in this species.
With three species, Young et al. (2014) presented a simple cursorial gradient (hare > rabbit > pika) based upon gross behavioural differences. Increasing the sample size requires a higher-resolution hypothetical gradient to rank cursorial ability among six species. Three possible gradients were derived (Table 2) based upon maximum running speed, degree of facial tilt and propensity to burrow. Maximum running speed values were obtained from the literature ( Table 2), but should be considered with caution. The values provided by Garland (1983) are drawn primarily from sources of unknown methodology, a fact acknowledged by the author. However, the pacing of Howell (1944) supports that L. californicus is unlikely to reach or exceed the reported speed of L. europaeus, conferring credibility to the placement of both hares on the gradient in this study. The angle of craniofacial tilt in leporids has been shown to correlate with cursoriality (Kraatz & Sherratt, 2016), whereby a more ventrally tilted face may provide greater visual acuity at speed. Species average angles were taken from Kraatz and Sherratt (2016) for L. californicus, S. bachmani, O. cuniculus. For Oc. princeps and S. bachmani, facial tilt angle was measured by one author (ES) from lateral diagrams of the 1 3 skull presented in Smith and Weston (1990) and Chapman (1974) respectively. For L. europaeus, a species average tilt angle was measured from the CT scans by ES. Propensity to burrow is included chiefly for comparison with the anatomical mechanical advantage measurements. The force vs velocity trade-off associated with mechanical advantage (wherein adaptation for running and burrowing are viewed as oppositional and mutually exclusive) means that a high propensity to burrow is likely to be negatively associated with cursoriality, and vice versa (Elissamburu & Vizcaíno, 2004;Smith & Savage, 1956). Propensity to burrow was inferred through qualitative descriptions in the literature, and although the differences between rabbits and hares overall are well-documented, there is a possibility the differences within the Sylvilagus and Lepus genera may be relatively negligible.

Samples, Imaging and Landmarking Protocol
We collated postcranial skeletal data from 280 lagomorph specimens (Table 3), with adults identified via the presence of epiphyseal and pelvic fusion. Specimens of L. europaeus and O. cuniculus were measured for this study from the preserved carcasses of 90 individuals acquired via population control culls during October/November 2020, from Wallangarra, Dalveen, Toowoomba, Oakey and Gatton, Queensland, Australia. Specimens of the remaining four species were obtained as skeletonised samples in natural history museum collections previously measured in Young et al. (2014). Specimens of S. floridanus were measured by one author (CDF) using methods identical to those used for all other species reported in Young et al. (2014).
Of the 90 Australian specimens, 16 exhibited damage to the skeleton that prevented certain measurements being taken, however this occurred in no more than two bones per specimen in each case. In the Young et al. (2014) dataset, some skeletons were incomplete, particularly missing the 3 rd metatarsal. A summary of the sample sizes by species for each bone is given in Table S2.

Skeletal Variables
From these landmark data, 14 linear variables were calculated across the bones, and nine variables were derived from combinations of the measurements following Young et al. (2014). All analyses were performed in R Statistical Environment (R Core Team, 2021), using the geomorph v. 4.0.0 Baken et al., 2021) and corpcor v. 1.6.10 (Schafer et al., 2021) packages.
We calculated three measures of distal limb segment elongation, referred to as relative length indices (RLI; Table 4). In the forearm, the brachial index provides a measure of how elongate the radius is relative to the humerus. Likewise, in Table 1 A description of the average weight and several aspects of ecology/behaviour of the six study species 1 Bock (2020), 2 Chapman & Flux, 1990, 3 Stott 2015, 4 Best 1996, 5 Melo-Ferreira et al. 2012, 6 Costa et al. 1976, 7 Chapman 1974, 8 Walker 1975, 9 Chapman et al. 1980, 10 Thompson & King 1994 Cowan & Bell 1986, 12 Smith & Weston 1990, 13 Huntly et al. 1986 Flee to cover 9 Flee to cover 11 Flee short distances to cover 8 1 3 the hindlimb, the crural index provides a measure of elongation in the tibio-fibula relative to the femur. Finally, the metatarsal-to-femur index provides an indication of elongation in the foot. The six measures of mechanical advantage used in this study are identical to those used by Young et al. (2014) and are named primarily for the muscle involved in the respective joint action. Each measure consists of a ratio of the length of the skeletal element that acts as a biological in-lever to the element acting as an out-lever (Table 5). This metric is referred to as anatomical mechanical advantage (AMA), and provides an indication of the impact of the bone architecture on joint action in the absence of muscle data (Carrier, 1983;Young, 2005). Comparisons with muscular mechanical advantage data gathered by Williams et al., (2007aWilliams et al., ( , 2007b, as well as previous studies (Young, 2009), show that such AMA proxies provide values that correlate well with more comprehensive musculoskeletal assessments of mechanical advantage. AMA was calculated for six joint actions in this investigation; extension of the foreleg/elbow (triceps brachii), stifle (quadriceps femoris), and hock (triceps surae), as well as three different actions at the hip (iliosoas, hamstring and lesser gluteal) (Table 5). A more comprehensive explanation of the principle of mechanical advantage is provided in the Supplementary Materials.
We evaluated the strength of correlation between morphological variables and the order of species predicted by each cursorial gradient using a Spearman's rank-order correlation (Kendall, 1949). Significance of the correlation coefficients was obtained through a permutation analysis, wherein the observed correlation coefficients for each gradient were compared to a distribution obtained from 64 random gradients (species ranked in a randomised order).

Morphological Integration
Modules are commonly identified using patterns of statistical covariation, where measurements within a module will exceed the covariation of variables between modules  1 3 (Armbruster et al., 2014). We assessed functional integration at two levels: between the forelimb and hindlimb, and among individual bones (Fig. 2). Two statistical approaches were used to identify patterns of covariation for functional integration: two-block partial least squares (2B-PLS) (Rohlf & Corti, 2000) for modules comprised two or more variables, and partial correlation (Pearson, 1915) analysis for modules comprised of a single variable. These were implemented with the two.b.pls function in geomorph and the cor2pcor function in corpcor. Statistical significance in the 2B-PLS analysis was identified to the 5% level by permutation approach (1000 random iterations), and the partial correlation coefficients using edge exclusion deviance (EED) scores exceeding 3.84 (as per Magwene, 2001). Correction for body size differences among species was performed to prevent differences in overall size (scale) from obscuring meaningful differences in proportions. This was performed using the log-shape ratio method (Mosimann 1979), which was found by Jungers et al. (1995) to be the most effective geometric scaling technique for linear morphological data. A geometric mean is calculated for each specimen by calculating the product of all variables and dividing by the number of variables. Each measurement variable is then standardised through division by the geometric mean, and then logarithmically transformed. In this case the geometric mean, representing size, was derived from the humerus, radius, tibia and femur lengths due to cases of incomplete data in other variables.

Relative Limb Indices
Variation in the three relative limb length indices (RLI) among species are shown in Fig. 3, with species arranged according to the maximum speed hypothetical gradient. The hares (Lepus spp.) display higher mean relative limb index (RLI) scores than the rabbits (Sylvilagus spp. and Oryctolagus) in all but two instances (L. europaeus is exceeded by S. bachmani in the crural index and by S. floridanus in metatarsal/femur index), and the pika (Oc. princeps) shows the lowest mean RLI in all three variables.
Contrary to the predictions of the maximum speed hypothetical gradient, relative distal limb segment lengths do not increase in a stepwise manner from the slower to faster species. A better alignment is often provided by the hypothetical gradient based on craniofacial tilt angle, as L. californicus consistently displays higher RLI than L. europaeus. However, results are more mixed in the three rabbit species. Contrary to the predictions of this gradient, O. cuniculus does not display higher RLI than the two Sylvilagus species, and S. bachmani does not display lower RLI than S. floridanus.
The low values of Oc. princeps align with both gradients. Correlation coefficients for each variable and the three hypothetical gradients (Table 6) indicate that the brachial and crural indices best align with the craniofacial tilt hypothesis (ρ = 0.860 and 0.488), and metatarsal-femur index the burrowing propensity hypothesis (ρ = 0.508).

Anatomical Mechanical Advantage
The distribution of the anatomical mechanical advantage (AMA) values are shown in Fig. 4. The highest AMA at the foreleg joint (triceps brachii) is displayed by the pika, and the lowest by the two species of hare (Lepus spp.). The rabbit O. cuniculus displays higher AMA values than predicted by the maximum speed and craniofacial tilt hypotheses, particularly in the hamstring. Similarly, L. europaeus consistently displays AMA values higher than L. californicus, and either equivalent to or higher than both Sylvilagus species in several cases. Most notably, the AMA of L. europaeus in the lesser gluteal is higher even than in the pika, Oc. princeps, which (based on two of the three hypothetical gradients) is predicted to, and generally does, display the highest AMA values of all sampled species.
We find the mean AMA scores of most variables do not correspond explicitly with the patterns predicted by any of the hypothetical cursoriality gradients. The relatively low correlation coefficients given in Table 7 demonstrate low support for the three hypothetical gradients amongst most of the AMA variables. The triceps brachii and iliopsoas AMA variables best align with the maximum speed gradient, the quadriceps femoris and lesser gluteal with the craniofacial tilt gradient, and the hamstring and triceps surae with the burrowing propensity gradient. However, the permutation analysis of these results revealed a large number of randomised hypotheses were better correlated to the variables than any of the three hypothetical gradients, particularly in the case of the triceps surae and lesser gluteal AMA measures (see Figure S1).

Between Fore-and Hind-Limbs
Scatterplot visualisations of the first PLS axis of the hindlimb against the first PLS axis of the forelimb (Fig. 5) demonstrate variation among species in the patterns of morphological integration between limbs. The slope of the PLS fit indicates the proportional relationship between the two limbs, whilst the distribution of points around the line of best fit provides the strength of integration (denoted by r-PLS). All species except Oc. princeps returned statistically significant degrees of integration between the fore-and 1 3 hindlimb. Integration between the limbs was greatest in the highly cursorial L. californicus, and lower in the less cursorial rabbits. Oc. princeps shows the lowest r-PLS value, however, this correlation is rejected at the 5% significance level (P-value = 0.243). The second specialised cursor, L. europaeus, exhibits a lower degree of integration between limbs than L. californicus, however, the value displayed by the rabbit S. bachmani is lower still. The rabbit S. floridanus displays a higher degree of integration than the two other rabbits. Patterns of integration between the fore-and hindlimbs are similar between S. bachmani and Oc. princeps, between L. europaeus and O. cuniculus, and between S. floridanus and L. californicus, as demonstrated by the angles of the slopes (Fig. 5). The similar distribution observed in adult and juvenile specimens indicates between-limb integration patterns are generally consistent throughout ontogeny.

Among Limb Bones
The results of the partial correlation analysis (Fig. 6) shows that in each species there are significant correlations between all serially homologous limb bones, except the metapodials in L. californicus, O. cuniculus and S. bachmani. Strength of the correlation between each pair of serial homologues varies between species; strongest between the humerus and femur (stylopods) in O. cuniculus and L. europaeus, between the tibio-fibula and radius (zeugopods) in S. bachmani, between the scapula and pelvis (girdles) in L. californicus, and between the metapodials in S. floridanus and Oc. princeps. Conversely, only three instances of correlation between functionally analogous bones reached statistical significance: the humerus/tibio-fibula as well as the radius/ metatarsal in L. californicus and the humerus/tibio-fibula in Oc. princeps. The patterns of integration within the limbs also differ notably between species. In L. californicus, a high degree of integration is present in the hindlimb, whereas no significant pattern was detected in the second hare, L. europaeus. Similarly, three instances of within-limb integration are detected in the two Sylvilagus rabbits, but only one is detected in Oryctolagus.

Discussion
Expanding the species sampling of Young et al. (2014) to include six lagomorph species reveals a more complicated association between cursoriality and appendicular skeletal adaptations than previously considered. In terms of distal limb elongation, the predicted stepwise increases from pika to rabbit to hare were present in the forelimb RLI values, but levels of elongation in the hindlimb showed a weaker directional trend. The predicted pattern was also present in the foreleg joint velocity (triceps brachii AMA), but inconsistent in all hindlimb joints. This indicates that mechanical advantage is likely affected in a complex, multifactorial manner and may not be suitable as a direct correlate for cursorial ability. Morphological integration varies among the study species, and the hypothesis of decreased integration between limb sets in cursors was generally unsupported. The hypothesis that increasing cursoriality would displace integration between serially homologous bones in favour of functional analogous bones was also shown to have low support, indicating integration patterns of lagomorphs are overall congruent with other therian taxa.

Distal Limb Segment Elongation will be Present to a Higher Degree in Hares than Rabbits, Both of Which Will Exceed the Pika
Our study identified the brachial index as the strongest correlate with cursorial behaviour based on our hypothesised cursoriality gradients. The crural and metatarsal-to-femur indices deviated from the predicted trends, with the metatarsal-to-femur index presenting the lowest correlation coefficients with each hypothesised gradient. This is contrary to both Young et al. (2014), wherein the metatarsal-to-femur index was the best fit for the three-species gradient, and to its wide application in the literature (Garland & Janis, 1993). One likely factor in the disparity of the hindlimb indices comes from the RLI values of L. europaeus and O. cuniculus; based on both the maximum speed and craniofacial tilt cursoriality hypotheses, these species are predicted to exhibit greater RLI values than the other sample species of their type.
For O. cuniculus, the low RLI values are potentially a result of burrowing propensity; the length of the femur is derived from a landmark placed on the femoral head (see Fig. 1 & Table 4), which is elongated in O. cuniculus to enable greater abduction of the hip for the rabbit to shovel soil between its hindlimbs more effectively (Fostowicz-Frelik, 2007). It is likely an increased femur length, rather than a decrease in elongation of distal elements, has caused the lower RLI values in O. cuniculus. This has implications for the mechanical advantage variables discussed in the next section.
Lepus europaeus also exhibits lower RLI values in its hindlimb morphology than expected given its speed and behaviour (Garland, 1983), however, this disparity is unrelated to fossorial ability as this species does not burrow (Bock, 2020). An alternative explanation may be body size; maximum stride length is limited by body size, which drives a disproportionate increase in limb length amongst smaller cursors (Coombs, 1978). Furthermore, evidence exists to support that increasing body mass necessitates changes 1 3 in limb proportion to decrease breakage risk, particularly in cursors (Gambaryan 1974;Biewener, 1989;Steudel & Beattie, 1993;Schmidt & Fischer, 2009). Having the largest size in this sample-and amongst most wider lagomorphs (Chapman & Flux, 1990)-suggests L. europaeus approaches the threshold at which size impacts this anatomical relationship. A similar phenomenon in distal elongation was noted by Gambaryan (1974), wherein the measurements  Table 5 The joint levers from which anatomical mechanical advantage (AMA) is calculated in this investigation, as well as the actions performed by these joints and the bone lengths from which they are derived # denotes landmarks in Fig. 1 Joint Action In-lever Out-lever

Triceps surae
Plantar flexion of the hindlimb at the hock Calcaneal tuber (#24-#25) + 3rd metatarsal (#22-#23) 3rd metatarsal (#22-#23) 1 3 of relative foot and tibial length of L. europaeus were less than in the mountain hare (L. timidus), despite the former being the more cursorial animal. Gambaryan (1974) attributed this discrepancy to the necessity of L. timidus to leap at a higher angle through thickly forested habitats, however, the proportions given for the Tolai hare (L. tolai) also exceeded L. europaeus, despite it favouring open regions (Smith et al., 2010). The average weight range of L. europaeus (3.5-5 kg) exceeds that of both L. tolai (1.6-2.6 kg) and L. timidus (2.9-4.7 kg in the largest subspecies), presenting the possibility that increased body mass affects the degree of distal limb elongation in leporids (Angerbjörn & Flux, 1995;Chapman & Flux, 1990;Smith et al., 2010). A mass-based explanation for the lower RLI values observed in L. europaeus also has implications for mechanical advantage; the lengths of the bones used in calculation of the hindlimb RLI act as the biological out-levers in several AMA variables. As such, any factors working to increase the AMA of the joints by limiting elongation in these species will cause the RLI values to decrease, as elongation occurs preferentially in the distal segments. Body mass is also strongly correlated with mechanical advantage itself (Gambaryan 1974;Biewener, 2005), the implications of which will be examined in the following section.

Hypothesis 2: Anatomical Mechanical Advantage
The Highest Anatomical Mechanical Advantage Values will be Observed in the Pika, and Hares will Tend Towards the Lowest.
The triceps brachii AMA (Fig. 4) supported the predictions of Hypothesis 2, where the maximum speed and craniofacial tilt hypothetical cursoriality gradients correlated negatively with this variable. However, the remaining five AMA variables were neither congruent with Hypothesis 2, nor strongly correlated to the finer-resolution predictions of the hypothetical cursoriality gradients. This is particularly evident in the lesser gluteal, wherein two correlation coefficients revealed a trend opposite to the prediction, and the permutation analysis ( Figure S1) revealed these observed values were no different to any random ranked order.
Relatively stronger support for the burrowing propensity hypothetical gradient supplied by the triceps surae and hamstring AMA measurements is likely to be influenced by the high AMA values in O. cuniculus. The burrowing style of the rabbit requires strengthened abduction of the hip, which depends in part on the action of the biceps femoris that originates on the ischial tuberosity of the pelvis (Fostowicz-Frelik, 2007;Williams et al., 2007b). This factor could also contribute to the remarkably high hamstring AMA value observed in O. cuniculus, but does not explain all other disparity in our results. Of note, an increase in the triceps brachii AMA is not observed in O. cuniculus. An elongated olecranon process provides greater force during digging (Winkler et al., 2016) and would be expected to increase triceps brachii AMA if present.
As stated in the previous section, disparity in AMA should be considered with body mass. Both Fig. 4 (AMA) and Fig. 3 (RLI) show that O. cuniculus and L. europaeus display higher/lower values than would be expected given their respective levels of cursoriality. As decreases in distal limb elongation will be reflected by increases in AMA, this may imply a factor that works to increase AMA is operating on both of these species-this is potentially body mass. Based on mean body mass, both European lagomorphs are larger than others of their morphotypes (Table 1). The larger an animal is, the more stress is placed on its bones by the muscle force required for locomotion. Without adaptations to counter this, bones become prone to fracture. Biewener (1989) identified an increase in effective mechanical advantage in the fore-and hindlimb of galloping animals as body size increases. Increasing effective mechanical advantage in lieu of muscle force decreases the stress placed on the bones without the joint losing the torque necessary for locomotion. In large animals, this leads to a transition from a crouched to an upright posture. As smaller animals, lagomorphs 1 3 maintain a relatively crouched posture, however Schmidt and Fischer (2009) mentioned 3-4 kg as a point where limb kinematics begin to differ in mammals-a range which the larger lagomorphs occupy (see also; Reilly et al., 2007). Lagomorphs also fall well below the weight threshold suggested by Bertram and Biewener (1990) after which skeletal proportions are strongly affected by allometric scaling, however, this 80-100 kg threshold is based on analysis of carnivores, bovids, and ceratomorphs. This raises questions regarding the limb scaling patterns of lagomorphs, given their rather unique form and function. Changes in the skeletal arrangement that affect AMA are clearly present in these results, which provides some grounds for further inquiries.
Our study supports the observation of Young et al. (2014) that AMA at the three hip joint actions, particularly in the hip extensors, was lower in Lepus than predicted given their cursorial capabilities. Gambaryan (1974) stated that the hip extensors bear the heaviest load during the hind support phase of the half-bound, from which forward propulsion is most strongly generated. The reliance on hindlimb propulsion in the asymmetrical gait used by cursorial lagomorphs may cause adaptations that reduce the stress of muscle force to be more pronounced in the hindlimb than the forelimb, particularly as this gait offers less force stability than the gallop (Gambaryan 1974;Álvarez et al., 2013). This potentially explains why the triceps brachii AMA of the forelimb-which produces less propulsion-is the sole AMA variable that is strongly correlated to the hypothetical gradients proposed in this study. Furthermore, calculations of gear ratio for the gastrocnemius muscle action (i.e., triceps surae AMA) by Fostowicz-Frelik (2007) also found that L. europaeus displayed a comparatively low gear ratio (i.e., high AMA) in a species sample that included O. cuniculus and S. floridanus, which the authors attributed to increased body mass. Young et al. (2014) observed decreases in bone robusticity as body size increased across their sample, a trend that is converse to the general expectation in mammalian species (Bertram & Biewener, 1990). This adaptation is proposed to decrease the weight of the limb and thusly increase speed, but places the bones at an increased risk of fracture. If these robusticity trends are true to other leporids, it is possible that the fracture risk may be countered by an increased mechanical advantage at the joints where muscle force is highest. This provides another potential explanation for the high AMA values of the hares in this investigation, as well as how lagomorphs can possess bone robusticity trends Fig. 3 Boxplot series showing the changes in distribution of the three ratios of distal limb elongation (RLI) among species. Species are arranged from left to right as per the order of the hypothetical cursoriality gradient based on maximum speed (see Table 2 that are converse to wider mammalian trends. Combined, these factors support the hypothesis that increased body mass may impact hindlimb AMA in lagomorphs, despite similar sized-based phenomena being rare amongst other orders (Day & Jayne, 2007;Schmidt & Fischer, 2009).

More Specialised Cursors will Display Less Integration Between the Forelimb and Hindlimb, as well as Greater Integration of Functionally Analogous Bones Compared to Serially Homologous Bones.
Our analyses of morphological integration do not support the hypothesis that more specialised cursors will display less integration between the forelimb and hindlimb, as evidenced by the results of the 2B-PLS analysis. Instead, lagomorphs appear to present a similar pattern as observed carnivores, where increased between-limb integration in limb geometry is found in cursors (Martín-Serra et al., 2015). This pattern purportedly arose due to the tendency for both limb pairs to be placed under similar biomechanical requirements, as their functions are both cursorial in nature. The differentiation of hindlimbs and forelimbs as propulsors and stabilisers respectively occurs in carnivorans when running (Álvarez et al., 2013), meaning the degree to which this gait pattern is used may have less effect on limb integration than Hypothesis 3 predicts. Martín-Serra et al. (2015) (and Young & Hallgrímsson, 2005 also suggest integration is lower in non-cursorial carnivorans as their forelimbs often perform non-locomotory tasks such as grappling. This presents the possibility that burrowing propensity in O. cuniculus would cause a decrease in integration between limb sets. Curiously, Fig. 5 shows O. cuniculus displays a level of integration similar to L. europaeus and higher than the non-burrower S. bachmani. This may indicate that cursorial bone length proportions convey a greater advantage to O. cuniculus than developing stronger fossorial ability, which is congruent with its lower than predicted triceps brachii AMA values. However, this investigation cannot satisfactorily conclude this point. The patterns observed between the fore-and hind-limb provide some support for the proposition of Schmidt and Fischer (2009) that postcranial integration in therians has evolved to optimise consistent proportions in the hindlimb over the forelimb due to their reliance on hindlimb propulsion. The forelimb demonstrates greater within-species variation than the hindlimb, implying that the proportion of the bones in the hindlimb is conserved more strongly than in the forelimb. Interestingly, the greatest variation in PLS1 (Fig. 5) of the hindlimb is shown in the cursorial L. californicus, which could indicate that strongly cursorial lagomorphs are less reliant on hindlimb proportions alone for locomotion. A greater incorporation of both limbs into the gait of more cursorial animals may also explain why between-limb integration is the highest in this species.
Regarding the differing levels of between-limb integration within lagomorphs exclusively, the 2B-PLS analysis shows similar slope patterns exist between disparate pairs of the sample species (S. bachmani/Oc. princeps, L. europaeus/O. cuniculus, and S. floridanus/L. californicus). If higher between-limb integration and greater forelimb variation are indeed indicative of greater cursoriality in lagomorphs, these results are congruent with the other variables measured in this study with regards to the lower amount of cursorial adaptation displayed by L. europaeus and O. cuniculus relative to the other species in their respective morphotypes. This may also explain the similarity of S. floridanus to L. californicus, as S. floridanus displayed marginally higher degrees of cursorial adaptation than the other rabbit species in several variables evaluated in this investigation.
Our results reject the hypothesis of greater integration among functionally analogous bones compared to serially homologous bones. We found only three instances of significant covariation between functional analogues, specifically in two modules of the highly cursorial L. californicus, but also in the generalist Oc. princeps (Fig. 6). Furthermore, the correlation coefficients were low compared to those between the serial homologues. Martín-Serra et al. (2015) also identified functional analogues were not consistently more strongly integrated than serial homologues in carnivorans; they found covariation between the humerus and tibia, but less so between the scapula and femur (attributed to their disparate developmental origins). This is congruent with the results of this investigation, wherein two in three cases of significant integration between functional analogues exist between the humerus and tibio-fibula, and none between the scapula and femur. However, in the serial homologues, Martín-Serra et al. (2015) observed stronger covariation between the zeugopods than the stylopods. This is converse to the pattern observed in lagomorphs, wherein the stylopods are more integrated than the zeugopods in all species except S. bachmani (Fig. 6). Similarly, low integration existed between the girdles and other bones in  carnivores, whereas several significant and strong instances of integration between the girdles and other bones exist in lagomorphs. These discrepancies indicate that although overall between-limb integration patterns may be similar for lagomorph and carnivoran cursors, they likely arise from differing fine-scale patterns. However, comparable geometric data in lagomorphs are needed to ascertain this finding. Previous studies of morphological integration in the limbs of mammals have found integration to be strongest between both bones with shared developmental origin, and bones that form functional units (Fabre et al., 2014). Our study suggests lagomorphs generally favour the presence of modules comprised of the humerus and tibio-fibula in the more cursorial species (Lepus spp. and Oryctolagus), and between the zeugopods and metapodials in the less cursorial species (Sylvilagus spp. and Ochotona). Of note, the scapula also shows greater integration with the rest of the forelimb bones in the less cursorial species. These observations may corroborate the discussion within Schmidt and Fischer (2009), which notes the zeugopods as being the least variable in proportion amongst mammals, whereas there is a fluctuation of the proximal and distal elements across species. The functional unit comprised of the proximal elements may therefore be of greater relevance in cursors, whereas generalists may favour consistent proportions in the distal elements.
A caveat exists in the sampling sizes of the morphological integration analysis in this study. The metapodial lengths were excluded from the between-limb integration analysis, as the large numbers of missing measurements (see Table S2 in the Supplementary Materials) caused the sample size of three species to drop below the 20-specimen accuracy threshold identified in the rarefaction analyses of Garland et al. (2017). Similarly, the partial correlation values shown in Fig. 6 relating to the metapodials come from reduced sampling of all North American lagomorphs. S. floridanus in particular only meets the 20-species threshold even in the wider partial correlation and 2B-PLS analysis, meaning the high between-limb integration of this species may be a false positive. Nevertheless, we have attempted the first investigation of morphological integration in postcranial skeleton of lagomorphs and provided some testable hypotheses for future studies. To date, between-limb variation in bone proportion has been shown to vary remarkably across therian taxa (Schmidt & Fischer, 2009), and even subtle changes in function within closely related species have purportedly led to distinct arrangements of limb integration (Martín-Serra et al., 2015). This complexity, reflected in the results of this investigation, shall hopefully encourage wider taxonomic studies to be done to better understand mammalian locomotor diversity.

Conclusion
It is evident that a cursorial gradient derived from observations of three lagomorphs by Young et al. (2014) is complicated in the context of a wider phylogenetic sampling, potentially as a result of the greater range of body sizes and burrowing habits introduced by the species in this study. This finding challenges the applicability of directional 'sliding scale'-type gradients to specialised morphological adaptation, particularly in the presence of differing body sizes. Such gradients can imply the existence of monotonic, directional changes in morphology across phylogeny, which-although possible (see Sherratt et al., 2016 and references therein)-can be hard to detect through univariate means and are often simply not present. The morphology of an animal is impacted at a holistic level by all features of its ecology such that adaptation in one aspect of its form will likely impact other, less related aspects in the way that morphological integration theory predicts. It is likely this consideration applies to many morphological correlates of cursoriality outside of this investigation. There are several avenues for the expansion and verification of this study; primarily, expanding the sample size even further will provide an indication of whether the inconsistency of the RLI and AMA trends is present across the existing diversity of lagomorph sizes and ecologies, or if the anomalous values of O. cuniculus and L. europaeus are unique to these two species. It may also be worthwhile to assess native European samples of these species to ascertain whether the anomalous results in this investigation are a result of adaptation to Australian ecosystems. Certainly, there is much to be learned from lagomorphs regarding the evolution of specialised locomotor adaptation.