Ontogeny of effective mechanical advantage in eastern cottontail rabbits (Sylvilagus floridanus)

Sub-adulthood can be a risky time of life for many animals. Smaller body sizes, immature musculoskeletal tissues and general ecological naïveté all compromise locomotor performance in juvenile animals (Carrier, 1996). Despite these limitations, juveniles must often compete and survive in the same environment as adults. Demographic studies of age-based mortality in natural populations have shown that ecological stresses on sub-adult prey taxa are particularly acute, with the rate of death from predation being highest in the first few weeks of life (Case, 1978; Promislow and Harvey, 1990; Adams et al., 1995). Therefore, strong selection for mechanisms that allow juvenile prey to survive past this demographic bottleneck should be expected, permitting them to survive to reproductive adulthood.

The ratio r/R quantifies the effective mechanical advantage (EMA) at a given limb joint. Maintaining all else equal, increases in EMA will necessarily lead to a relatively high output force for a given muscle force (Smith and Savage, 1956; Gray, 1968; Biewener, 1989, 1990, 1991; Roberts et al., 1998).

Increasing extensor muscle EMA, either through growth-related changes in r or though postural adjustments to shorten R, represents one hypothetical strategy that juvenile animals could use to overcome absolutely weaker limb muscles and other growth-related limits on locomotor performance (Carrier, 1983, 1996; Young, 2005). Here, we test three specific predictions related to this hypothesis. Firstly, we predicted juvenile cottontail rabbits would have greater hindlimb joint EMA than adult conspecifics, either through negative ontogenetic allometry of r (i.e. a greater EMA numerator in young animals), positive ontogenetic allometry of R (i.e. a smaller EMA denominator in young animals) or a combination of these two scaling trends. Secondly, we predicted that ontogenetic variation in R will be associated with variation in joint posture, GRF orientation, limb length distal to the joint or a combination of these factors. Finally, as a functional consequence of greater hindlimb muscle EMA, we predicted that juvenile S. floridanus would require less active muscle volume to support their body weight and power locomotion (proportional to the metabolic cost; Roberts et al., 1998; Biewener et al., 2004; Pontzer et al., 2008).

We obtained approval from the Northeast Ohio Medical University (NEOMED) Institutional Animal Care and Use Committee for all procedures (protocols 10-032 and 13-026, to J.W.Y.) and secured Scientific Collecting permits from the Ohio Division of Natural Resources Division of Wildlife (permits 14-310, 15-173 and 16-128, to G.A.S.) prior to initiating this research.

Eastern cottontail rabbits (S. floridanus) were trapped at several sites in Portage, Summit and Columbiana Counties, Ohio, USA, where we secured prior permission for trapping from site owners or controlling authorities. All sites were public parks, nature preserves or private properties where hunting was forbidden. We used a combination of wooden rabbit traps (Wildlife Control Supplies, East Granby, CT, USA) and medium-sized wire traps (Tomahawk Model 606SS, Tomahawk Live Trap, Hazelhurst, WI, USA). Traps were baited with a variety of lures, including oats, peanut butter, apples and rabbit urine, set in the early evening and checked the subsequent morning.

We transported captured animals to the Portage Park District Breakneck Creek Field Station, where we recorded body mass to the nearest gram, shaved and marked major hindlimb joints with small circles of retro-reflective tape for later video tracking, and measured hindlimb segment lengths to the millimeter using digital calipers. A list of joints marked, landmarks used for identification and hindlimb segment definitions are presented in Table 1.

We trapped a total of 61 rabbits over three field seasons from 2013 to 2015. Owing to variation in animal condition and motivation, and occasional equipment malfunction, we were only able to collect locomotor performance data in 42 rabbits. Of the animals contributing locomotor performance data, we euthanized 14 animals for detailed anatomical study, including measurements of r for several hindlimb extensor muscles and fiber architecture (see below; Butcher et al., 2019). The sample for the present study consists of a subset of 13 individuals for whom we obtained sufficient kinematic and morphological data to calculate EMA. Overall, the ontogenetic S. floridanus dataset used here extended over an order of magnitude in body mass (i.e. 0.106–1.277 kg). Individual S. floridanus with a body mass <1 kg were designated as ‘juveniles’; animals with a body mass >1 kg were designated as ‘adults’. This 1 kg body mass cut-off corresponds to the age at which body mass growth begins to reach asymptotic values and long bone epiphyses fuse (Hale, 1949; Lord, 1963). Based on a published body mass growth curve (Lord, 1963), we estimated that rabbits in our sample ranged from juveniles approximately 22 days of age to fully grown adults. Following these criteria, the animal sample for this study included eight juvenile and five adult rabbits.

Following joint marking, rabbits were placed in a dark transport box that was secured to the end of an enclosed runway (4.0×0.2×1.0 m, length×width×height). Rabbits were coaxed to move from the transport box to a corresponding box at the other end of the runway using mildly aversive stimuli (e.g. puffs of air from a compressed air canister) to encourage maximum effort. Two HE6×6-16 small animal force plates (AMTI, Watertown, MA, USA; 15×15 cm in area) were embedded in the floor of the runway to measure GRF during bounding strides of burst acceleration. Forces were sampled at 500 Hz using proprietary software (NetForce, AMTI). The force plates were covered with stairway tread tape (Safety-Walk Slip Resistant Tread, 3M Corp., Minneapolis, MN, USA). Locomotion was filmed with two high-speed cameras (Fastec TS3 100-L, Fastec Imaging, San Diego, CA, USA) also using proprietary software (HiSpec, Fastec Imaging). The cameras were placed approximately 45 deg to the direction of travel and positioned approximately 90 deg relative to one another. Force platform and camera recordings were synchronized by means of a common trigger.

Video and force platform data were processed in MATLAB R2018a (MathWorks, Natick, MA, USA). We used the open source motion-tracking software DLTdv5 (Hedrick, 2008) to track the position of all landmarks throughout each video clip. The 3D position of each anatomical landmark was reconstructed using direct linear transformation (DLT) (Hedrick, 2008). DLT coefficients were generated using DLTcal5, a companion program to DLTdv5, by digitizing a control object of known geometry that had been filmed immediately prior to each locomotor data collection. The orientation axes of the calibration object were made coincident with those of the force platforms. We also used the DLTdv5 software to demarcate the beginning and end of the hindlimb push-off phase of the stride (i.e. from hindlimb touchdown to lift-off). Only full- or half-bounding strides, where the separation between leading and trailing hindlimb contacts were separated by ≤10% of hindlimb contact duration (Hildebrand, 1977), and where the rabbit was fully supported by one or both force platforms for the duration of hindlimb push-off, were included in the dataset.

Synchronized force and kinematic data were imported in a custom-written MATLAB routine for additional processing. Joint coordinates were fit to a quintic smoothing spline function (tolerance of 0.75 mm²), allowing us to mitigate digitizing error and interpolate the position of a feature for frames where the marker was not visible (Walker, 1998). Force traces from each channel of the force plates were smoothed using a zero-lag fourth-order Butterworth low-pass filter. Filter cut-off frequencies were selected using an automated procedure developed by Winter (2005), designed to balance noise reduction against signal attenuation [mean selected frequency±95% confidence interval (CI)=51±2.9 Hz]. Force plate baseline drift during the experiment was corrected by sampling the average values of unloaded periods immediately prior to and following contact and subtracting these values from the force traces. Further analyses only consider movements in the sagittal plane (i.e. we ignored mediolateral forces and displacements). Overall, average mediolateral force magnitudes were about half of fore–aft force magnitudes (mean±95% CI=54.8±25.4%). Moreover, given the narrow width of the testing runway, mediolateral impulse could not contribute to animal displacement in the net direction of travel, and should therefore be considered ‘wasted effort’.

where M_x is the moment about the x-axis (i.e. mediolateral axis) of the force platform, F_y and F_z are the forces in the fore–aft and vertical directions, respectively, and a_z0 is the thickness of the vertical offset of the force platform's cover plate relative to the position of the sensors. Instantaneous values of M_x, F_y and F_z were determined using the NetForce software, whereas a_z0 was provided by the force plate manufacturer. Given that, by definition, the half-bounding gaits used by the rabbits involve simultaneous contact by the right and left hindlimbs – which may or may not be spatially aligned – the raw CoP position indicated by the force plate calculation alone did not always sensibly align with the kinematic coordinates (which were only digitized on the left side of the animal). The CoP position was therefore set to the midpoint of the distance between the metatarsal head and the distal tip of the third pedal digit at the beginning of the support phase, and then allowed to subsequently translate, as indicated by the frame-to-frame changes in CoP position calculated from the force plate data (with the constraint that the CoP could never move beyond the tip of the third digit or proximal to the metatarsal head). We also assumed a common fore–aft CoP position across both limbs in subsequent calculations of joint moments (McGowan, 2005).

where ρ represents a standard value of density (1.06 g cm⁻³) for mammalian skeletal muscle (Méndez and Keys, 1960).

Average muscle lever arm lengths at the hip, knee, and ankle joints (i.e. r_h, r_k and r_a) were measured as the mean of all individual extensor lever arm lengths at the joint, weighted by the PSCA of the corresponding muscle. This method assumes that muscle stress (i.e. force per unit area) is distributed equally across the muscles at a joint. Additionally, this method assumes that lever arm lengths remain constant through changes in joint position. Data on instantaneous changes in muscle lever arm lengths in European hares (Lepus europeus, a closely related leporine taxon) indicate that for the muscles studied here, values of r remain fairly constant over the ranges of joint excursion observed in our S. floridanus locomotor dataset, particularly relative to changes in GRF load arm lengths (Williams et al., 2007).

where t₁ and t₂ represent the beginning and end of the support phase interval. This quotient indicates how much force, on average, was exerted on the substrate per unit muscle force. This method has the benefit of characterizing mean EMA throughout the support phase, rather than at finite instances or intervals, and allows for the variation in R that can result from changing joint postures throughout the support phase.

where muscle abbreviations follow those listed in Table 2. This set of equations was solved for each instance (i.e. video frame) of the support phase, first solving Eqn 6 and subsequently solving Eqns 7 and 8 as a set of simultaneous linear equations (using MATLAB's ‘linsolve’ function). We then calculated total muscle force impulse at each joint as the finite integral of instantaneous muscle forces throughout the support phase, and then ultimately mean EMA over the support phase as the ratio of GRF impulse to muscle impulse. Total EMA across the hindlimb was calculated as the mean of individual joint EMA, weighted by the total extensor muscle PCSA at that joint.

As the goal of this study was to investigate how developmental changes in EMA impact locomotor performance (i.e. external work against the environment), we did not include calculations of the muscular effort required to resist gravitational moments or perform the internal work required to accelerate limb segments relative to the CoM (Wells, 1981; Biewener and Full, 1992; Winter, 2005). Though such components of total limb work could impact total energy expenditure, studies of non-human mammals have shown that inertial and gravitational moments forces are minor relative to the external moments engendered by GRF, and therefore have less influence on estimates of muscle torque and leverage (Witte et al., 2002; Dutto, 2004; McGowan, 2005; Williams et al., 2009).

Following Roberts et al. (1998), we calculated L_mean as the mean of individual muscle fascicle lengths weighted by their PCSAs. We then calculated V_musc as the quotient of L_mean and EMA.

The locomotor data for this study were explicitly collected to measure peak acceleration performance in S. floridanus and were therefore not collected at steady-state velocities. Given that acceleration has the potential to impact EMA through an influence on limb posture and GRF orientation (Roberts and Scales, 2002; Biewener et al., 2004; Williams et al., 2009), we included average acceleration (calculated as average fore–aft force throughout hindlimb support phase divided by body mass) as a potential predictor of EMA.

Average R across the support phase was calculated as the quotient of r and EMA (algebraically equivalent to dividing muscle moment by GRF moment; Eqns 1 and 2; Fig. 1).

Joint angles were calculated as the two-dimensional vector angle between the relevant limb segments, with increasing values indicating greater limb joint extension. We calculated mean hip, knee and ankle joint angles throughout support phase to evaluate how changes in joint posture affected GRF load arm lengths (R).

where F_y and F_z refer to the GRF components in the fore–aft and vertical directions, respectively.

Analyses of the full dataset (i.e. where each individual was represented by multiple trials) were carried out using mixed-effects models (Pinheiro and Bates, 2000), where individual rabbit was included as a random factor. To improve mixed-effect model fits, raw variates were first Box–Cox transformed to improve normality (Box and Cox, 1964; Sokal and Rohlf, 1995), and then scaled and centered (i.e. converted into z-scores). Conversion to z-scores also permitted direct comparisons of the resulting partial regression coefficients (i.e. β weights) from a multiple regression to evaluate which predictors best explained variance in the dependent variable. Because the residuals of such models were often characterized by heteroscedasticity, the variance of the error term was allowed to vary as an exponential function of the independent (i.e. predictor) variable. Coefficients of determination (R²) for these models were calculated following Johnson (2014). We used standard log–log Model II (reduced major axis) regressions on body mass to test for ontogenetic allometry of muscle lever arm lengths, average GRF load arm lengths and limb lengths. These analyses were conducted on reduced datasets, where only a single value existed for each individual. Significance for all tests was accepted at P≤0.05. All statistical analyses were conducted in the R statistical platform (version 3.5.2, ‘Feather Spray’) (https://www.r-project.org/), including the add-on packages dplyr (https://CRAN.R-project.org/package=dplyr), emmeans (https://CRAN.R-project.org/package=emmeans), MASS (Venables and Ripley, 2002), MuMIn (https://CRAN.R-project.org/package=MuMIn), nlme (http://CRAN.R-project.org/package=nlme) and rptR (Stoffel et al., 2017).

EMA at the hip and ankle declined with significant negative allometry in S. floridanus, whereas knee EMA remained unchanged (i.e. was not significantly associated with body mass) (Table 3, Fig. 2A–C,E–G). Given that EMA declined at two of the three joints examined, total hindlimb EMA also significantly declined with increasing body size, the reverse of scaling patterns previously observed in interspecific studies of adult mammals of varying body size (Fig. 2D,H). Controlling for variation in body mass, increasing acceleration decreased hip EMA, increased knee EMA and had no significant effect on EMA at the ankle or total hindlimb EMA (Table 3).

As detailed above (see Fig. 1), EMA can be increased via longer r or shorter R. Therefore, an ontogenetic decrease in EMA implied negative allometry of r (i.e. greater values in younger animals), positive allometry of R (greater values in older animals) or a combination of these scaling trends. Values of r scaled with positive allometry at hip and negative allometry at the knee and ankle joints (Table 4, Fig. 3A–C), though the 95% confidence intervals about these slopes included the isometric expectation of 0.33 in each case. Inversely, R scaled with significant positive allometry at each joint (i.e. the 95% confidence intervals did not include isometry; Table 4, Fig. 4D–F). We also examined the ontogenetic scaling of distal limb length (i.e. thigh length at the hip, leg length at the knee and foot length at the ankle; see Table 1 for segment definitions) as a possible determinant of ontogenetic changes in R (see below). Thigh and leg length scaled with positive allometry, whereas foot length scaled with slight negative allometry (and the 95% confidence intervals included isometry; Table 4, Fig. 3G–I).

Because R is defined as the perpendicular distance between the GRF vector and the joint center of rotation, its magnitude can vary as a function of joint angle, GRF orientation or limb length distal to the joint (Biewener, 1989; Polk, 2002; Young, 2009). Multiple mixed-effects regression models indicated that R was positively associated with distal limb length, though the relationship was only significant at the hip and ankle joints (Table 5). Neither average joint angle nor θ_GRF were significantly associated with variation in R at any of the joints analyzed.

Mean hindlimb extensor muscle fascicle length (i.e. L_mean) of the hindlimb extensors scaled with positive ontogenetic allometry in S. floridanus, though the 95% confidence interval on the estimate included isometry [slope (95% CI)=0.48 (0.303, 0.749), R²=0.50, P=0.007]. Active muscle volume (V_musc) increased with body mass and, when controlling for body mass, decreased with increasing total hindlimb EMA (Table 6). Mean acceleration was not independently associated with V_musc.

Our hypotheses were generally supported by our analyses. EMA at the hip and ankle joints, as well as in the hindlimb overall, significantly declined during S. floridanus ontogeny. Lack of significant change in knee EMA is likely due to the position of the knee joint relative to the GRF vector. Given that our locomotor dataset explicitly sampled acceleration, the average GRF vector had a net cranial orientation (mean±95% CI=67±6.9 deg), indicative of the net propulsive function of the hindlimbs during the support phase. As a result, the GRF vector typically passed very near the knee joint, falling slightly cranial or slightly caudal to the joint throughout the support phase. The variable position of the GRF with respect to the knee resulted in highly variable GRF load arm lengths at the knee joint (see Fig. 3E), highly variable EMA ratios, and therefore no evidence of consistent ontogenetic change in knee joint EMA. The variable position of the GRF vector relative to the hip and knee joints also accounts for observed effects of acceleration on EMA at these joints. At greater accelerations, the average GRF vector tilted cranially, increasing R at the hip joint but decreasing R at the knee joint, resulting in a negative correlation between acceleration and hip joint EMA and a positive correlation between acceleration and knee joint EMA.

Overall, our results are consistent with previous research documenting negative allometry of anatomical mechanical advantage (i.e. bony muscle lever arm length scaled to skeletal limb length below the joint) in other mammalian taxa (e.g. jackrabbits, Carrier, 1983; capuchin monkeys, Young, 2005), with attendant benefits for juvenile locomotor performance (Carrier, 1995). Given the combination of morphometric, kinematic and kinetic data required, EMA has only been measured in a small sample of adult animals (Biewener, 1989; Carrier et al., 1998; Blob and Biewener, 2001; Polk, 2002; Biewener et al., 2004), and we are aware of only one other study (Smith and Wilson, 2013) that has attempted to evaluate EMA across ontogeny. Smith and Wilson (2013) measured overall hindlimb EMA in growing ostriches (Struthio camelus) over the first 10 months of life (during which time body mass increased from 5 to 75 kg), finding no consistent changes in mean hindlimb EMA with increasing size (although the lowest EMA was observed in largest ostriches; their fig. 4). Discrepancies in the ontogenetic trajectory of EMA between cottontail rabbits and ostriches could be due to several factors, principally (1) differences between avian and mammalian limb anatomy and posture (Gatesy and Pollard, 2011) and (2) differences in overall body size. As the largest extant birds, adult ostriches would require high EMA to overcome the geometric limits on muscle force and bone bending strength faced by all terrestrial organisms (Biewener, 1989, 1990, 1991; Gatesy, 1991; Reilly et al., 2007). It could also be that the apparent lack of size-related change observed in the Smith and Wilson (2013) study arose from assumptions made during their analyses. Specifically, rather than calculate r at each age, they assumed that lever arms scaled proportionally to segment lengths, using allometric equations from their previous studies of ostrich growth and musculoskeletal anatomy (Smith et al., 2007, 2010). Because R also closely tracked segment length during ostrich growth (Smith and Wilson, 2013), it may be that constant EMA values during ontogeny arise from calculating a ratio of two values that fundamentally track the same dimension.

Declining EMA in the hindlimb of S. floridanus during growth resulted from isometric ontogenetic scaling of muscle lever arm lengths combined with positive ontogenetic allometry of GRF moment arm lengths. In other words, the denominator of EMA (see Eqn 2) increased with body size, resulting in increasingly smaller ratios as animals grew. This is the reverse of the pattern observed in interspecific studies of adult mammals of varying body size, where size-related increases in EMA result from the negative allometry of R (in this case, principally due to the tendency of larger animals to stand and move with more extended limb postures) (Biewener, 1989, 1990, 1991; Gatesy, 1991; Reilly et al., 2007; but see Ren et al., 2010). Comparison of observed values with predictions using published scaling equations of adult mammals (Biewener, 1989, 1990; Biewener et al., 2004) specifically show that young S. floridanus have much greater hindlimb EMA than predicted for a generalized mammal of their size, whereas EMA in adult S. floridanus more closely aligns with predicted values (Fig. 2H). As such, EMA may represent an example of size-related allometric changes during ontogeny inversely mirroring those observed in interspecific studies. The ontogeny of long bone cross-sectional robusticity represents another example of this trend. Young animals typically have geometrically more robust bones than their adult counterparts, resulting in negative ontogenetic allometry of second moments of area, section moduli and other geometric correlates of bone strength (e.g. Carrier, 1983; Currey, 1984; Brear et al., 1990; Carrier and Leon, 1990; Heinrich et al., 1999; Currey, 2001; Lammers and German, 2002; Ruff, 2003; Main and Biewener, 2004; Main and Biewener, 2007; Young et al., 2010; Russo and Young, 2011; Kilbourne and Makovicky, 2012; Patel et al., 2017), whereas geometric measures of bone strength increase with positive allometry across the size range of extant mammals (Alexander et al., 1979).

Future work should seek to characterize ontogenetic changes in EMA in other terrestrial animals, particularly given the differing findings of this study and Smith and Wilson (2013). It may be possible to make inferences about the ontogenetic scaling of EMA in other mammals by studying the ontogenetic allometry of limb bone lengths. In our study, we found that distal limb length was the best predictor of R across the ontogenetic S. floridanus dataset, with positive allometry of limb lengths generally leading to positive allometry of GRF load lengths. Young (2009) and Smith and Wilson (2013) made similar observations with respect to growing squirrel monkeys (Saimiri boliviensis) and ostriches, respectively. In fact, an examination of scaling exponents from allometric studies of limb growth in 26 mammal species spanning six orders of magnitude in body mass supports the prediction that positive allometry of limb bone lengths – and, by extension, GRF load arm lengths – is characteristic of mammalian ontogeny (Fig. 4). Specifically, proximal and middle limb segments (i.e. brachium/thigh and antebrachium/crus) generally scale with positive allometry (proximal: 81–87% of taxa; middle: 72–79% of taxa; Fig. 4B–E), whereas distal segments (i.e. manus/pes) generally scale with negative allometry (75–77% of taxa; Fig. 4F,G). Scaling exponents for the proximal segments tend to decrease with increasing body mass, whereas those for distal segments tend to increase with body mass, though most of these trends were not statistically significant, particularly when analyzed in a comparative framework using phylogenetic generalized least squares regression. Overall, these data suggest that, in general, young mammals are characterized by relatively short limbs for their size, though distal segments tend to be relatively long early in life (see also Jungers and Fleagle, 1980; Young and Heard-Booth, 2016; Druelle et al., 2018; Young et al., 2019). Inasmuch as distal limb length determines R, ontogenetic declines in muscle leverage may also be characteristic of mammals in general.

What might be the functional consequences of ontogenetic declines in EMA? We found that greater hindlimb extensor EMA significantly decreased the estimated muscle volume required to be activated during locomotion in young S. floridanus. Because active muscle volume is directly proportional to metabolic cost (Roberts et al., 1998; Biewener et al., 2004; Pontzer et al., 2008), this finding suggests that burst locomotion may have been energetically more efficient in young cottontail rabbits. Mitigating the energetic cost of locomotion would allow young rabbits to shunt more energy to tissue growth and differentiation, accelerating growth and ultimately limiting the time that must be spent in the ecologically risky state inherent to being a small-bodied juvenile (Williams, 1966; Werner and Gillam, 1984; Janson and van Schaik, 1993; Carrier, 1996; Arendt and Reznick, 2005; Herrel and Gibb, 2006). Additionally, greater EMA necessarily reduces the muscular forces that the skeleton must withstand during locomotion, perhaps increasing limb bone safety factors in young rabbits (Biewener, 1989; Main and Biewener, 2004). Similarly, by increasing the output forces resulting from a given muscle contraction, greater extensor muscle EMA may facilitate heightened acceleration capacity in juvenile S. floridanus, a fitness-critical performance metric during predator evasion (Carrier, 1983, 1995; Walker et al., 2005). Overall, the findings of this study strongly suggest that allometric growth trajectories, and resulting performance capacities, are indicative of adaptive processes acting on juveniles and not just adults (Carrier, 1996).

The authors wish to thank Christine Craycroft and the Portage Park District for permission to use the Breakneck Creek Field Station for the duration of this study. We also thank Summit Metro Parks, the University of Akron, and multiple private landowners in Portage, Summit and Columbiana counties, Ohio, USA, for permission to trap rabbits on their properties. Many high school, undergraduate, graduate and medical students contributed to data collection and processing, including: E. Barkette, B. Barnette, K. Burns, B. Crawford, F. Galloway, Z. Glenn, E. Grubb, S. Hamrick, B. Herbert, E. Kunzler, A. Metz, E. Moreland, K. Reardon, J. Rose, J. Rupert, M. Slavens, N. Tatomirovich, A. Thakore and D. Thomas. Connie Fellmann provided valuable discussion and assistance during the initial planning phases of this work. We thank the two external reviewers and the editor for their valuable comments on an earlier version of this manuscript.