Limb bone scaling in hopping macropods and quadrupedal artiodactyls.

Bone adaptation is modulated by the timing, direction, rate and magnitude of mechanical loads. To investigate whether frequent slow, or infrequent fast, gaits could dominate bone adaptation to load, we compared scaling of the limb bones from two mammalian herbivore clades that use radically different high-speed gaits, bipedal hopping (suborder Macropodiformes; kangaroos and kin) and quadrupedal galloping (order Artiodactyla; goats, deer and kin). Forelimb and hindlimb bones were collected from 20 artiodactyl and 15 macropod species (body mass M 1.05-1536 kg) and scanned in computed tomography or X-ray microtomography. Second moment of area (I max) and bone length (l) were measured. Scaling relations (y = axb ) were calculated for l versus M for each bone and for I max versus M and I max versus l for every 5% of length. I max versus M scaling relationships were broadly similar between clades despite the macropod forelimb being nearly unloaded, and the hindlimb highly loaded, during bipedal hopping. I max versus l and l versus M scaling were related to locomotor and behavioural specializations. Low-intensity loads may be sufficient to maintain bone mass across a wide range of species. Occasional high-intensity gaits might not break through the load sensitivity saturation engendered by frequent low-intensity gaits.

Introduction

During daily rest and activity in development, growth and adulthood, bones experience a range of mechanical loading conditions that relate to each behaviour's physical intensity. Kangaroos, wallabies and their macropodiform kin are famed for their hopping hindlimb gait which they use for bursts of efficient high-speed locomotion [1-3]. They are less well known for their slower pentapedal gait, wherein their powerful tail acts as the third point of a tripod with the forelimbs during hindlimb protraction [4] (figure 1). The pentapedal gait is used during grazing and, along with other slow-speed activities, dominates macropods' locomotor behaviour. Tammar wallabies (Macropus eugenii) spend up to twice as much time per day in pentapedal walking than in bipedal hopping (6% versus 3–5%), and both gaits are eclipsed in the locomotor time budget by bipedal standing (50–70%), quadrupedal crouching (15–30%) and bipedal rearing (3–12%) [3,5,6]. During hopping, the forelimbs are held away from ground contact for the entire stride cycle and thus are relatively unloaded [2], while hindlimb tissues experience near-ultimate stresses from ground reaction forces and muscle-tendon action, especially in larger Macropodiformes [7]. The tail's role in pentapedal locomotion during slow-speed locomotion might enable reduced forelimb mass, potentially assisting more efficient bipedal hopping [4]. In extinct sthenurine macropods, the thoracic limb displays features of a browsing adaptation with elongated manus, reduced lateral digits, slender radius, ulna and humerus, and a ‘human-like’ scapula, which may have enabled these animals to forage above their heads [8]. Hopping is probably not possible at body mass over approximately 160 kg, at which the distal tendons' safety factor (ratio of actual to ultimate stress) drops below 1, meaning that extinct ‘giant kangaroos’ would have used slower gaits [7-10].

Figure 1.

Bennett's wallaby (Macropus rufogriseus) in the hindlimb suspension phase of the pentapedal gait (a) and blackbuck (Antelope cervicapra) in a lateral sequence walk (b), indicating the limb bones measured in the study. These two species have femora of similar length (199 mm and 186 mm, respectively) and are presented here approximately to scale. Drawing by Manuela Bertoni may be reused under the CC-BY licence.

In contrast to macropods, artiodactyl mammals (even-toed ungulates in the eutherian lineage; deer, sheep, camels and kin) have limited manual dexterity, and quadrupedal gaits in which the loads are spread more evenly among fore- and hindlimbs during both slow and fast gaits, reflected in similarity of forelimb and hindlimb bones' cross-sectional properties [11]. Artiodactyls and macropods spend a large proportion of their time grazing or resting as they are foregut fermenter herbivores [12] and may be considered ecological equivalents [13]. Extinct giant macropods probably could not hop due to tissue strength being exceeded by scaling of muscle and tendon stress [9], yet even the largest artiodactyls retain high-speed gaits. Bison, buffalo and giraffe are capable of galloping [14,15], while hippopotami achieve high land speeds by a fast walk or trot [16]. Scaling of limb bones in artiodactyls is relatively well characterized, exhibiting isometric or modestly allometric patterns across their three orders of magnitude body mass range [17-20].

The distributions of occasional maximal loads and habitual moderate loads vary within the skeleton and depend on locomotor activity, which should appear as a morphological signal in clades that adopt very different characteristic gaits [21]. Although direct bone strain gauge data do not exist, positive allometry of hindlimb muscle physiological cross-sectional area, reduced duty factor with increasing speed and constant effective mechanical advantage of hindlimb joints probably lead to relatively increased muscle force, and increased stress and reduced safety factors in larger macropods’ hind limb bones [10,22,23]. Bennett [24] pointed out that kangaroos' tibial cross-sections (section modulus Z and second moment of area I, which relate to fracture strength and resistance to bending, respectively) scale more strongly than other quadrupeds [24], whereas McGowan et al. [22] found that the macropod femur is more robust in larger animals lending support to the concept that intense hopping could relate to increased hindlimb robustness. Musculotendinous forces generated during hopping could incur relatively larger loads on tendon insertion sites around the metaphyses compared to artiodactyls. Those larger loads in macropods may manifest as stronger scaling of cross-sectional parameters in macropods’ metaphyses, evidenced as higher scaling exponents. Conversely, if the more frequently used, slower gaits' loading environment drives bone shape, then we should expect to see similar scaling between macropods’ fore- and hindlimbs, and between equivalent bones in macropods and artiodactyls, because the low-speed pentapedal gait and quadrupedal walking, respectively, dominate these clades' locomotor repertoires.

Using artiodactyls as a quadrupedal comparator clade that has a wide body mass range, we ask whether macropod limb bones exhibit structural scaling features that relate to their pentapedal and hopping locomotor specializations. In particular, we predict that the forelimb bones of the macropods, which are used for grasping and low-speed locomotion (and are essentially unloaded during hopping), should become relatively more gracile with increases in body size, and consequently have lower scaling exponents than artiodactyl forelimb and macropod hindlimb bones. We hypothesize that scaling exponents should be more similar between fore- and hindlimb bones in artiodactyls than in macropods due to artiodactyls’ more even distribution of stresses between fore- and hindlimbs during high-speed locomotion.

Material and methods

We selected the humerus, radius, ulna and metacarpal bone (III in macropods and fused III–IV in artiodactyls), along with the femur, tibia and metatarsal bone (IV in macropods and fused III–IV in artiodactyls) from 15 macropod and 20 artiodactyl species (table 1). All specimens were skeletally mature as determined by fusion of the growth plates and free from obvious skeletal disease. We imaged the bones in clinical computed tomographic (CT) scanners (LightSpeed 16, Ultra, or Pro 16, GE Medical Systems, Pollards Wood, UK) or for the smallest specimens, in an X-ray microtomographic scanner (XMT; X-Tek HMX ST 225, Nikon Metrology, Tring, UK) with the bone's long axis positioned parallel to the image's z-axis, and applied a similar image processing technique used elsewhere [25,26]. Scans where the long axis of the bone was oblique to the z-axis of the scanner were aligned with BoneJ's Moments plugin, so that the bone's principal axes of inertia were parallel with the scan's x-, y- and z-axes. Scans with large numbers of image slices were downsampled without interpolation to contain 100–200 slices, providing 5–10 values for averaging in each 5% length bin. Fat in the marrow cavity and other bony or metal elements were manually replaced with a pixel value corresponding to air. Where nearby or fused bones could not be excluded by a rectangular region of interest (ROI), they were manually removed by replacing them with pixels of an air-equivalent value. Bones containing lesions or severe post-mortem deterioration were excluded from the study. Image analysis was performed with BoneJ v. 1.4.2 [27,28] for ImageJ v. 1.51c [29].

Table 1.

Complete list of specimens, their body masses and lengths of bones used for scaling calculations. Some bones from some specimens were not available to study and are indicated as a dash. NHM, Natural History Museum (London); UMZC, University Museum of Zoology, Cambridge; RVC, The Royal Veterinary College (authors' collections); mc, metacarpal; mt, metatarsal.

						bone length (mm)
source	accession number	clade	binomial	familiar name	M (kg)	humerus	radius	ulna	mc	femur	tibia	mt
UMZC	H15052	Artiodactyla	Tragulus kanchil	lesser mouse-deer	1.62	68.7	57.0	67.7	39.3	78.3	91.0	58.4
UMZC	H14975	Artiodactyla	Tragulus napu	greater mouse-deer	5.27	84.7	—	—	49.8	101.1	117.1	70.6
UMZC	H.15532	Artiodactyla	Muntiacus reevesi	Reeves’ muntjac	13.50	93.3	80.2	99.3	59.6	114.7	123.4	83.9
NHM	ZD 1863.12.29.1	Artiodactyla	Cephalophus dorsalis	Bay duiker	20.00	124.2	106.9	132.3	85.8	152.8	151.7	97.8
NHM	ZD 1974.414	Artiodactyla	Antilope cervicapra	blackbuck	36.30	151.4	171.9	207.0	185.4	185.7	227.8	192.6
RVC	alpaca1	Artiodactyla	Vicugna vicugna	alpaca	45.00	—	—	—	—	265.8	—	—
RVC	goat1	Artiodactyla	Capra hircus	goat	48.00	—	—	—	—	216.5	—	—
RVC	sheep2	Artiodactyla	Ovis aries	domestic sheep	57.00	—	—	—	—	192.8	—	—
NHM	ZD 1961.8.9.80	Artiodactyla	Cephalophus sylvicultor	yellow-backed duiker	61.30	193.2	176.7	219.6	155.4	232.3	242.6	170.2
UMZC	H.16232	Artiodactyla	Cervus davidianus	Pere David's deer	182.00	255.2	261.5	314.0	244.6	327.8	339.4	263.9
NHM	ZD 1963.10.21.1	Artiodactyla	Oryx beisa	beisa	200.58	243.2	273.2	336.5	224.1	309.4	322.6	232.9
UMZC	H.20302	Artiodactyla	Okapia johnstoni	okapi	230.00	317.5	350.2	418.0	309.4	328.4	353.1	324.7
UMZC	H.16634	Artiodactyla	Cervus elaphus	red deer	240.87	—	266.9	314.2	231.6	278.0	329.0	265.1
NHM	47	Artiodactyla	Bos indicus	zebu	275.00	178.3	190.5	229.1	—	221.8	228.8	—
UMZC	H.17691	Artiodactyla	Alces alces	Eurasian elk	461.90	387.8	401.1	491.5	343.0	445.9	492.7	406.4
RVC	WindfallLgGiraffe	Artiodactyla	Giraffa camelopardalis	giraffe	482.00	410.3	628.0	416.3	587.7	457.9	525.8	580.5
RVC	cow2	Artiodactyla	Bos taurus	cattle	500.00	—	—	—	—	415.9	—	—
UMZC	H.17535	Artiodactyla	Megaloceros giganteus	Irish elk	585.00	405.1	402.1	501.3	334.9	—	477.6	370.5
NHM	ZD 1874.11.2.4	Artiodactyla	Syncerus caffer	African buffalo	592.67	347.8	316.9	404.0	187.3	411.0	387.5	216.0
UMZC	H.10707-H.10715	Artiodactyla	Hippopotamus amphibius	hippopotamus	1536.31	453.7	276.4	384.9	—	474.4	339.0	—
NHM	ZD 1851.4.24.2	Macropodiformes	Potorous tridactylus	long-nosed potoroo	1.05	39.3	45.1	55.6	12.0	75.2	83.4	30.3
NHM	ZD 1858.5.26.23	Macropodiformes	Bettongia penicillata	woylie	1.18	33.6	44.1	52.7	10.1	80.1	104.4	44.4
NHM	ZD 277.p	Macropodiformes	Bettongia lesueur	boodie	1.45	22.8	28.5	35.6	6.9	56.8	75.6	38.9
UMZC	A12.79/1	Macropodiformes	Bettongia gaimardi	eastern bettong	1.67	33.3	38.2	47.6	—	84.4	99.5	46.4
NHM	ZD 1984.1002-1003	Macropodiformes	Macropus parma	parma wallaby	4.16	70.1	85.2	98.4	14.1	134.2	168.5	57.9
UMZC	A12.59/2	Macropodiformes	Onychogalea fraenata	bridled nail-tail wallaby	4.95	71.9	79.4	92.6	—	119.5	185.0	56.3
NHM	ZD 1962.5.22.1	Macropodiformes	Macropus eugenii	tammar wallaby	5.28	56.6	74.8	86.4	12.0	123.3	167.4	60.1
UMZC	A12.50/1	Macropodiformes	Thylogale billardierii	Tasmanian pademelon	5.87	78.3	91.5	107.5	—	134.7	165.1	52.1
UMZC	A12.52/2	Macropodiformes	Petrogale xanthopus	yellow-footed rock-wallaby	8.50	82.6	91.9	110.4	—	157.7	202.7	69.5
NHM	ZD 1970.2189	Macropodiformes	Macropus agilis	agile wallaby	11.82	85.3	108.4	126.0	20.5	169.8	239.4	94.5
NHM	ZD 1976.184	Macropodiformes	Macropus rufogriseus	Bennett's wallaby	16.85	126.2	138.4	157.2	19.4	199.3	285.1	98.6
NHM	ZD 1961.12.11.1	Macropodiformes	Macropus fuliginosus	western grey kangaroo	25.56	82.3	109.0	128.0	17.3	189.3	290.4	101.4
UMZC	A12.19/1	Macropodiformes	Macropus robustus	wallaroo	25.98	97.8	135.7	159.5	—	207.5	326.4	113.8
UMZC	A12.17/1	Macropodiformes	Macropus giganteus	eastern grey kangaroo	60.00	175.1	235.0	266.4	—	267.5	473.9	166.3
UMZC	A12.17/4	Macropodiformes	Macropus giganteus	eastern grey kangaroo	60.00	216.5	280.2	326.6	—	299.7	543.2	183.8
NHM	ZD 2010.8	Macropodiformes	Macropus rufus	red kangaroo	72.50	212.3	310.8	348.2	36.2	300.5	572.7	179.6

The maximum second moment of area (Imax) was measured on every slice of each scanned specimen with Slice Geometry in BoneJ. Other parameters including Imin, cross-sectional area (CSA) and section modulus were also measured and are available in the associated datasets [30], but are not reported here due to their close mathematical relationships: I is calculated by multiplying area by distance from the principal axis squared, and section modulus is calculated by dividing I by chord length. Imax can reflect the stiffness of a member in bending, although we use Imax here as a geometric parameter of cross-sectional size and disposition that includes non-bending-related features such as the tibial crest, which is loaded mainly in tension by the patellar ligament. To calculate bending stiffness, I is calculated around the neutral plane of bending, which moves during the stride cycle and which usually does not coincide with the principal axis used to calculate Imax [31]. Because the ratios between specimen size, image resolution and pixel spacing were not constant, we applied a correction for partial filling of pixels which maintains comparable cross-sectional area measurements when image resolution, pixel spacing and resolution vary with respect to each other (figure 2). Partial filling correction was set by excluding pixels less than −800 HU to eliminate artefacts with values close to air (−1000 HU) and scaling linearly between −1000 HU (0% bone, 100% air) to 2300 HU (100% bone). Pixel values over 2300 HU were considered 100% bone. XMT images lacking HU calibration were set by taking a histogram of an ROI positioned in the background and using its mean for the 100% air scaling value and its maximum as the minimum cut-off value. Another histogram was made in a thick region of cortical bone and its mean used as the 100% bone scaling value. The partial volume correction approach was validated using synthetic images and an exemplar CT image, and resulted in a high degree of stability compared to global thresholding; test scripts and data are available online [30]. Bone length (l) was measured using the image data, which we validated against physical measurement of the bones. Body mass (M) was unknown for most of the specimens so was estimated from literature values [32-35]. The red and eastern grey kangaroo specimens were male, so we used body masses near the high end of the estimate to account for the sexual dimorphism in these species.

Figure 2.

Interaction between specimen size, image resolution and pixel spacing. As pixel spacing increases and resolution decreases relative to specimen size, a greater proportion of pixels represent the edge of the specimens compared with the mid-substance. (a) Progressive downsampling of a well-sampled image of a bone cross-section (top left) increases pixel spacing (vertical axis), and Gaussian blurring with increasing radius simulates lower instrument resolution (horizontal axis). High-resolution images from X-ray microtomography (b) and lower resolution clinical CT images (c) relate to different pixel spacing/image resolution combinations within this scheme. We corrected for imaging condition and specimen size variation using a weighted pixel sum approach in BoneJ's Slice Geometry plugin.

We analysed scaling of bone dimensions using the general equation y = ax [36], where y is the bone parameter, x is a measure of size (body mass M or bone length l), a is the scaling coefficient and b is the scaling exponent. The exponent b expresses the rate of change in y as a function of body size, while a is the magnitude of y when x = 1. Scaling analysis relies on linear fitting to the log-transformed variables, log(y) = log(a) + blog(x), where b becomes the slope of the line and log(a) the y intercept or ‘elevation’. The purpose of a scaling analysis is to determine whether relative proportions of the variables under study (here, l, Imax and M) vary with size. Isometry occurs when the larger animal is a ‘to-scale’ version of the smaller animal and the dimensions retain the same proportions, with an isometric scaling exponent b. Deviations from isometry occur when b is greater or less than b, and the corresponding variation in proportions with size is termed allometry. Positive allometry occurs when there is a disproportionately greater increase in y with increasing x and b > b, while negative allometry occurs when there is a disproportionately lesser increase in y with increasing x, and b < b. Differences in scaling coefficient between groups may be calculated when exponents are equal, and represent a fixed ratio with changing size. Scaling estimates were calculated using smatr v. 3.4-3 [37] for R [38], using the standardized major axis (SMA, also known as RMA), which accounts for error in x as well as in y and is suited to scaling analyses, which typically deal with data that have a large amount of unidentifiable error or ‘natural variation’ [39]. Exclusion of the exponent values for isometry (l ∼ M, b = ; Imax ∼ M, b = ; Imax ∼ l, b = 4) by 95% confidence intervals returned by smatr was used to test the null hypothesis of isometry, and R2 and p from an F-test were used to test the null hypothesis of scale invariance, in which there is no correlation between variables. Cross-sectional parameters were averaged within each 5% increment of length, and scaling exponents and elevations were calculated for each 5% bin across all the individuals in each clade, for each bone in the study. Normalized cross-sectional parameters were calculated by dividing the nth root of the parameter by length. Second moment of area has units of mm4, so it was normalized by taking the fourth root and dividing by bone length in mm. Normalized parameters are unitless and a size-independent measure of shape.

To control for non-independence of samples due to their phylogeny, a phylogenetically informed version of RMA (phyl.RMA) from the R package phytools [40] was called from a custom script to calculate scaling relationships versus M for bone length and Imax at mid-shaft (50% of bone length). Calibrated phylogenetic trees, one for macropods and one for artiodactyls, were used for this analysis and were constructed based on divergence time estimates from a previous publication [41]; values from the two eastern grey kangaroo specimens were averaged for the phylogenetic analysis (figure 3). We report the results of the phylogenetic RMA with five parameters: the estimate of the scaling exponent (b); the squared correlation coefficient (R2); the p-value for deviation from isometry comparing exponents using an F-test (p); Pagel's λ, a measure of how well the phylogeny explains the data (λ = 1 suggests that the evolution of the observed traits follows a Brownian motion model with variance proportional to divergence times, whereas λ = 0 suggests phylogenetic independence) and the log-likelihood (L) [42] in table 3. We considered also using phylogenetic generalized least-squares (PGLS) analysis, but prior studies have found negligible differences between phylogenetic RMA and PGLS results or conclusions [43,44], so we opted to apply only phylogenetic RMA to our data, thereby retaining a closer methodological equivalence to our standard scaling analysis.

Figure 3.

Cladograms illustrating phylogenetic relationships (from Bininda-Emonds et al. [41]) among the artiodactyl (a) and macropod (b) species used to perform phylogenetic-independent contrast calculations.

Table 3.

Bone length and mid-shaft second moment of area phylogenetically informed RMA scaling exponents against body mass, where l ∝ M or Imax ∝ M. The scaling exponent (slope, b) is indicated, along with the coefficient of determination (R2) indicating the strength of the correlation between either l or Imax and body mass values, while p denotes the p-value for the isometry test (logarithm of estimated slope is compared to hypothesized value using a t-test, as described in [45]). λ denotes Pagel's lambda, a parameter quantifying phylogenetic signal, while L denotes log-likelihood of the fitted model.

	artiodactyl phylogenetic RMA scaling analysis					macropod phylogenetic RMA scaling analysis
bone	b	R2	p	λ	L	b	R2	p	λ	L
l ∝ Mb
humerus	0.354	0.917	0.468	1.000	−23.071	0.480	0.780	0.023	0.657	−19.970
radius	0.401	0.855	0.108	1.000	−27.833	0.489	0.843	0.005	0.562	−20.015
ulna	0.373	0.880	0.270	0.982	−25.363	0.471	0.853	0.007	0.553	−19.015
metacarpal	0.449	0.806	0.040	1.000	−27.234	0.340	0.726	0.919	0.514	−13.808
femur	0.314	0.894	0.455	0.548	−27.826	0.346	0.894	0.695	0.592	−11.671
tibia	0.330	0.858	0.917	1.000	−25.230	0.446	0.949	0.001	0.409	−11.396
metatarsal	0.401	0.824	0.158	1.000	−24.973	0.432	0.935	0.004	0.818	−9.837
Imax ∝ Mb
humerus	1.740	0.951	0.001	0.842	−43.370	1.726	0.640	0.164	0.669	−41.410
radius	1.688	0.909	0.017	0.827	−47.003	1.609	0.645	0.279	0.539	−44.176
ulna	1.777	0.914	0.005	0.976	−46.334	1.668	0.717	0.156	0.587	−42.659
metacarpal	1.725	0.956	0.002	0.675	−37.438	1.409	0.880	0.685	0.429	−23.174
femur	1.399	0.942	0.424	0.000	−54.017	1.496	0.908	0.199	0.591	−32.535
tibia	1.552	0.914	0.078	0.697	−48.396	1.469	0.927	0.243	0.419	−30.572
metatarsal	1.651	0.950	0.008	0.633	−37.981	1.564	0.941	0.038	0.407	−31.343

Results

Comparisons of bone length versus body mass, and 50 per cent length Imax versus l and M, are presented in table 2 and figure 4. Imax, l and M are strongly correlated in all bones with high R2 (0.737–0.995) and p < 0.001, excluding the null hypothesis of no scaling relationship between the variables. Humerus, radius and ulna lengths scale with positive allometry (b > ) in macropods, but with isometry (b not significantly different from ) in artiodactyls. Artiodactyl metacarpal bones are much longer than in macropods of similar mass, indicated by the high elevation (1.47 versus 0.86). In the hindlimb, femur and metatarsal lengths scale similarly in macropods and artiodactyls, with the macropod femur having a higher elevation than artiodactyls and the metatarsals' slopes and elevations not significantly different. Tibia length scales isometrically in artiodactyls and with strong positive allometry in macropods. Comparing stylopod (humerus and femur), zeugopod (radius, ulna and tibia) and autopod (metacarpal and metatarsal) elements between limbs within each of the two clades, there is a high degree of overlap between the confidence limits of scaling exponents in all the limb segments (table 2), meaning that bone length proportionality between fore- and hindlimb segments is maintained within clades. Mid-diaphyseal Imax scales against M with similar exponents between clades but with differing elevations in all bones except the humerus, indicating constant proportionality between clades with increasing animal size. Mid-diaphyseal Imax scales against l with different exponents between clades in the humerus, radius, ulna and tibia, and with different elevations in the metacarpal and metatarsal. Notably, the femora are indistinguishable at 50% length in their Imax ∼ l scaling.

Table 2.

Summary statistics for bone length (l) scaling against body mass (M), second moment of area (lmax) at 50% length versus M and Imax at 50% length versus l regressions for all bones. The scaling exponent (slope, b) is indicated alongside its upper and lower 95% confidence limits (±b), along with the coefficient of determination (R2) and p indicating the strength of the correlation between bone length and body mass values (i.e. the p-value of the F-test of the correlation coefficient against zero). Isometric scaling exponents (b) are listed for each comparison. The likelihood of equality of scaling exponents between artiodactyl and macropod bones is indicated by the result of a Wald test, pb. Where slopes are not significantly different, elevations (log10(a)), their 95% confidence limits (±log10(a)) and equality of elevations between artiodactyls and macropods (pa) are reported. Statistical estimates were generated by R using smatr calls in the form: >summary(sma(log10(length) ∼ log10(mass) * order)); >summary(sma(log10(length) ∼ log10(mass) + order, type = ’elevation’)).

	artiodactyl scaling exponent					macropod scaling exponent					equal slopes	artiodactyl elevation			macropod elevation			equal elevations
bone	b	−b	+b	R2	p	b	−b	+b	R2	p	pb	log10(a)	−log10(a)	+log10(a)	log10(a)	−log10(a)	+log10(a)	pa
l = aMb (bi = 13)
humerus	0.317	0.270	0.372	0.928	<0.001	0.463	0.384	0.558	0.892	<0.001	0.003
radius	0.357	0.275	0.464	0.804	<0.001	0.487	0.416	0.571	0.923	<0.001	0.042
ulna	0.333	0.272	0.408	0.884	<0.001	0.468	0.401	0.547	0.927	<0.001	0.009
metacarpal	0.425	0.337	0.535	0.863	<0.001	0.323	0.219	0.477	0.800	<0.001	0.198	1.470	1.298	1.642	0.858	0.739	0.977	<0.001
femur	0.299	0.258	0.346	0.917	<0.001	0.353	0.311	0.400	0.951	<0.001	0.788	1.735	1.661	1.810	1.858	1.817	1.899	<0.001
tibia	0.282	0.227	0.350	0.857	<0.001	0.447	0.408	0.491	0.974	<0.001	<0.001
metatarsal	0.369	0.294	0.462	0.870	<0.001	0.392	0.349	0.442	0.957	<0.001	0.608	1.545	1.440	1.651	1.519	1.472	1.566	0.500
Imax = aMb (bi = 43)
humerus	1.660	1.494	1.851	0.970	<0.001	1.651	1.295	2.104	0.833	<0.001	0.954	0.872	0.512	1.232	0.970	0.673	1.268	0.535
radius	1.586	1.344	1.871	0.923	<0.001	1.545	1.211	1.972	0.816	<0.001	0.855	0.755	0.253	1.258	0.281	−0.027	0.590	0.012
ulna	1.691	1.377	2.076	0.881	<0.001	1.630	1.320	2.013	0.862	<0.001	0.794	−0.859	−1.447	−0.27	0.143	−0.169	0.454	<0.001
metacarpal	1.655	1.484	1.845	0.979	<0.001	1.352	1.041	1.756	0.912	<0.001	0.134	0.546	0.197	0.894	−0.620	−0.928	−0.306	< 0.001
femur	1.399	1.237	1.583	0.942	<0.001	1.533	1.363	1.723	0.958	<0.001	0.268	1.212	0.927	1.497	1.701	1.540	1.862	<0.001
tibia	1.472	1.283	1.689	0.942	<0.001	1.499	1.338	1.679	0.964	<0.001	0.832	1.135	0.823	1.447	1.494	1.336	1.652	0.002
metatarsal	1.574	1.400	1.770	0.965	<0.001	1.550	1.399	1.717	0.968	<0.001	0.833	0.699	0.423	0.974	0.921	0.775	1.068	0.023
Imax = alb (bi = 4)
humerus	5.251	4.838	5.700	0.981	<0.001	3.539	3.222	3.887	0.975	<0.001	<0.001
radius	4.437	3.828	5.142	0.939	<0.001	3.171	2.703	3.720	0.922	<0.001	0.003
ulna	5.080	3.754	6.875	0.737	<0.001	3.480	3.067	3.949	0.951	<0.001	0.024
metacarpal	3.895	3.185	4.764	0.896	<0.001	4.188	3.291	5.329	0.926	<0.001	0.614	−5.338	−6.724	−3.952	−4.059	−4.860	−3.258	<0.001
femur	4.685	4.237	5.181	0.960	<0.001	4.342	4.176	4.514	0.995	<0.001	0.152	−6.337	−6.732	−5.943	−6.412	−6.763	−6.060	0.190
tibia	5.223	4.443	6.141	0.920	<0.001	3.352	3.149	3.569	0.989	<0.001	<0.001
metatarsal	4.270	3.576	5.100	0.920	<0.001	3.952	3.549	4.400	0.965	<0.001	0.433	−5.541	−6.438	−4.643	−5.210	−5.919	−4.501	0.013

Figure 4.

Bone length (l) versus body mass (M; a–g); second moment of area (Imax) at 50% length versus M (h–n) and max at 50% length versus l (o–u) regressions for all bones. Bone lengths and body masses are presented in table 1. Scaling exponents (slopes, b), elevations (log10(a)), R2 and p-values are presented in table 2.

Normalized Imax versus per cent length plots (figure 5) reveal that artiodactyls’ cross-sections become relatively more robust with increasing body mass, indicated by the larger animals' traces tending towards the top of the range. Meanwhile, macropods show the opposite trend, with normalized Imax decreasing with increasing body mass so that traces from the larger animals appear at the bottom of the range, indicating increased gracility with increasing body mass. In general, and in common with prior studies on cats and birds [25,26], the artiodactyl diaphysis occupies a decreasing proportion (and metaphyses and epiphyses increasing proportions) of bone length with increasing body mass, but this relationship is not maintained consistently in macropods. Notably, the trochlear notch and coronoid processes of the ulna drift distally in larger artiodactyls, but proximally in larger macropods (figure 5e,f).

Figure 5.

Normalized second moment of area (Imax1/4/l) at each fraction of length. Dispersal of traces indicates bone shape differences among species. Higher traces indicate relatively more robust bone geometry, seen in larger artiodactyls and smaller macropods. Note the increasing proportion of length occupied by epiphyseal and metaphyseal components in larger artiodactyl species (a,i), and the distal drift of the ulna's trochlear notch (dashed line) in larger artiodactyls (e) and smaller macropods (f).

Scaling exponents (figure 6) and elevations (figure 7) for Imax versus M reveal near-identical scaling exponents between clades for all regions of all the bones, and overlapping elevations for all bones in all regions except for the proximal tibial and femoral metaphyses, indicating very similar bone cross-sectional scaling against body mass. Positive allometry (exponent above the isometry line) is strongest in the proximal metaphyses, possibly relating to their increasing relative length mentioned above, and this is amplified by increased elevations (i.e. larger value of Imax at a given M) in these regions in macropods (figure 7f,l,n). Imax versus l scaling reveals positive allometry for much of the length of artiodactyl bones. The wide confidence interval of artiodactyl ulna (figure 6e) probably reflects the variability of fusion to the radius, reducing the strength of the body size signal. By contrast, macropod Imax scales with negative allometry against l for much of the length of humerus, radius, ulna and tibia, with positive allometry in the femur and isometry in the metacarpal and metatarsal. The raised elevation of macropods relative to artiodactyls in the Imax versus l plots (figure 7) is difficult to interpret because the scaling exponents are markedly different between clades in the regions where elevations are different. Despite their orders of magnitude difference, the elevations may not relate to functional differences, which may be more strongly indicated by differing scaling exponents.

Figure 6.

(a–n) Scaling exponents for second moment of area (Imax) versus bone length (l) and body mass (M) for all bones and both clades. Light grey regions indicate the 95% confidence interval; dark grey regions occur where the confidence intervals overlap and where substantial may be interpreted as no significant difference in scaling exponent in that region of the bone, between clades. Scaling exponent estimates at 50% length are presented in table 2.

Figure 7.

(a–n) Scaling elevations [log10(a)] for second moment of area (Imax) versus bone length (l) and body mass (M) for all bones and both clades. Elevations are directly comparable only where slopes (scaling exponents, figure 6) are equal. Light grey regions indicate the 95% confidence interval; dark grey regions occur where the confidence intervals overlap and where substantial may be interpreted as no significant difference in scaling elevation between clades. Scaling elevation estimates at 50% length are presented in table 2.

Despite there being a highly variable amount of phylogenetic signal (Pagel's λ ranged between 0.407 and 0.818 for macropod bones, while it covered the entire range of 0–1 for artiodactyl bones), scaling exponents for both artiodactyl and macropod bone length and mid-shaft Imax corrected for phylogenetic effects using phylogenetically informed RMA remained comfortably within the 95% confidence limits of the scaling exponents calculated without phylogenetic adjustment (tables 2 and 3). Correlations calculated using phylogenetic RMA were similarly strong as their non-corrected counterparts (unadjusted R2 = 0.800–0.979, adjusted R2 = 0.640–0.956).

In our phylogenetically informed RMA, macropod femur and third metacarpal lengths scale isometrically, while all other macropod bone lengths scale with positive allometry. Macropod mid-shaft Imax does not scale differently from isometry in all bones except the metatarsal (p: humerus 0.164; radius 0.279; ulna 0.156; metacarpal 0.685; femur 0.199; tibia 0.243; metatarsal 0.038).

Phylogenetically informed RMA suggests that artiodactyl bone lengths scale isometrically except for the metacarpal (slope not different from , p: humerus 0.468; radius 0.108; ulna 0.27; metacarpal 0.04; femur 0.455; tibia 0.917; metatarsal 0.158). Artiodactyl mid-shaft Imax scales with positive allometry in the humerus, ulna, metatarsal and metacarpal (p < 0.005), with a similar tendency in the radius (p = 0.017), but not in the femur or tibia (p: femur 0.424; tibia 0.078).

In summary, scaling exponents calculated using phylogenetic adjustments for bone length, mid-shaft CSA and mid-shaft Imax concur with the results from the unadjusted scaling relationships. Thus, phylogeny has varying explanatory value for bone geometry scaling (as quantified by Pagel's λ) within the two clades of mammals studied here, but does not influence the interpretation of the scaling exponents, because the estimates are essentially no different between phylogenetically informed and non-phylogenetic scaling analyses.

Discussion

Scaling of the forelimb and hindlimb segments is similar within clades, except the stylopod, in which the Imax versus length scaling displays positive allometry in the artiodactyl and macropod femur and artiodactyl humerus, but negative allometry in the macropod humerus, meaning that in macropods the humerus becomes more gracile with increasing length while the femur becomes more robust. Unlike Bennett [24], who found that tibial second moment of area scales more strongly positively in kangaroos than quadrupeds (b = 1.52 versus 1.28), our data show that tibial cross-sections scale similarly against body mass between clades. This may be a consequence of comparing macropods to artiodactyls only, and not to a more diverse sample of quadrupeds, because it is known that artiodactyls' bones scale differently to other mammalian clades [19,20]. Tibial cross-sections scale with strong negative allometry in macropods and positively allometrically in artiodactyls against length. This means that larger kangaroos’ tibiae are relatively less robust—they are relatively longer and more slender, consistent with a relatively reduced ability to resist bending moments. This apparent reduction in relative bending strength is surprising, considering that bending stresses predominate over compressive stresses due to the off-axis component of the muscular forces, with a stress range of −110 to −60 MPa and 90–110 MPa [1].

The Imax versus body mass scaling elevation is higher in the macropod proximal femur and tibia than the same region in the artiodactyl femur and tibia, indicating increased robustness around the greater and lesser trochanters and tibial crest, which are the bony insertions for the massive gluteal, iliopsoas and quadriceps muscle groups that drive bipedal hopping in macropods. Positive allometry of tarsal joint moment arms potentially ameliorates the musculotendinous compressive force on the tibia during tarsus extension [22], allowing the distal half of macropods' tibial cross-sections to remain within similar parameters as artiodactyls'. Like McGowan et al. [22] demonstrated in macropods, we find that the femur is more robust in larger macropods [22], which is consistent with a proposal of a strong, general relationship between stylopod cross-sectional parameters and body mass [21]. We find that humeral and femoral lengths scale significantly differently against body mass between macropods and artiodactyls, in contrast to suggestions of common mammalian femur length to body mass scaling [46], which may have implications for mid-shaft bending stresses.

The largest extant artiodactyls are an order of magnitude more massive than the largest extant macropods, while the smallest of both clades included in this study are approximately 1–2 kg. It would be unwise to extrapolate macropod scaling trends beyond the current series, because bipedal hopping was probably not a feature of the extinct giant kangaroos and may not be physiologically possible beyond approximately 160 kg [7-10]. Janis et al. suggested that large, extant kangaroos are functionally specialized for hopping in contrast to their larger extinct kin that may not have hopped, somewhat similar to the medium-sized, gracile and hyper-athletic cheetah (Acinonyx jubatus, M = 35–70 kg) compared with bigger and more robust felids such as lion (Panthera leo, M = 120–250 kg) [9].

We found that the trochlear notch of the ulna is relatively more distal in larger artiodactyls, but that an opposite trend of a relatively shortened olecranon process and proximally drifting trochlear notch is observed with increasing mass in macropods. We first noted a trend to a more centrally placed trochlear notch in the ulna of large felids [25] and proposed that this may be a mechanism that allows reduced muscle forces by increasing the lever arm of the olecranon process and increasing the elbow extensor muscles' effective mechanical advantage. The relatively shortening olecranon in larger macropods may relate to forelimb use in the low-intensity pentapedal gait and lack of loading in high-intensity bipedal hopping, or to reaching and combat behaviours favouring a longer forelimb. Interclade differences in metacarpal dimensions relate to their functional specializations for grasping, or plantigrade or unguligrade locomotion in macropods and artiodactyls, respectively. Similar isometric Imax scaling exponents against length indicate the maintenance of overall bone shape that may relate to specialized manus function, whereas positive allometry against body mass in artiodactyls but isometry in macropods may reflect an influence of locomotor loading on artiodactyl metacarpal robustness that is absent or reduced in macropods. The femur's cross-sectional scaling is remarkably conserved between artiodactyls and macropods (figures 4s, 6j and 7j), illustrated by low Pagel's λ (table 3), despite their phylogenetic separation and different gaits. The digits and zeugopod (radius, ulna and tibia) appear to be the dominant sites of evolutionary plasticity in limb bones, and this is borne out by their higher Pagel's λ values.

Bones respond anabolically, that is, by increasing bone tissue formation and decreasing bone resorption, when they experience a small number of novel high strain and high strain rate events with a rest period between bouts of loading [47,48]. Repetitive loading has a saturation or habituation effect, in which tissue is no longer responsive to mechanical loads after a few tens of cycles [48]. The lack of a difference in femoral and tibial Imax versus body mass scaling exponents between bipedal hopping macropods and quadrupedal artiodactyls suggests that the occasional very high load of hopping may not be sufficient to overcome the mechanobiological saturation engendered by frequent but lower intensity loading in crouching and pentapedal walking. Alternatively, bipedal hopping may be a no more intense stimulus to the hindlimb than quadrupedal galloping, leading to little discernible difference between clades. Galloping and hopping engender similar muscle stresses at preferred speeds in rats and kangaroo rats, respectively, despite a fourfold difference in ground reaction force [49]. We found similar Imax ∼ M scaling exponents in the forelimb bones, despite macropods' relatively unloaded forelimbs during bipedal hopping.

Limb bone adaptation occurs throughout the life of the individual due to bone's phenotypic plasticity, and at the population level over evolution due to selection pressure relating to skeletal development and functional specialization. Variable safety factors among species and bones [50-52] suggest that peak strains from uncommon or high energy gaits are not necessarily the dominant stimuli for bones' phenotypic adaptation, or a critical selection pressure in evolution. Bone is sensitive to its mechanical environment during growth: altered muscle and gait forces on growing bones relate to bone deformities in children with cerebral palsy and in experimentally manipulated chick and crocodile embryos [53,54], while 12 h per week of throwing practice leads to substantial cortical bone adaptation in children's throwing arms [55]. In vivo bone adaptation experiments have shown that the mouse tibia receives less than 300 µɛ (microstrain) during walking and less than 600 µɛ when landing a 30 cm jump (approx. 3 N physiological load), yet requires at least 1000 µɛ from a 10 N experimental load to stimulate further bone formation [56-58]. Only high-intensity race training, well beyond horses' ad libitum locomotor behaviour, was sufficient to engender an increase in cortical area [59]. Sciatic neurectomy removes daily habitual loading in the mouse hindlimb, sensitizing the tibia to subsequent load-induced (2000 µɛ) bone deposition [60], suggesting that the removal of background stimulus can rescue bone's load responsiveness [61], or in other words, that daily stimulus might saturate bone's ability to respond to further applied loads. Simple body mass support and the ground reaction forces incurred by a slow gait may be sufficient to maintain cross-sectional bone geometry, while infrequent high-intensity quadrupedal gaits might offer little further selection pressure or modelling stimulus to increase diaphyseal size over and above that provided by standing and walking.

There are few data on the daily numbers of stride cycles in each gait for the species in the study, which limits our ability to calculate bone loading histories and infer which gaits relate most strongly to bone structural scaling; however, in those species that have been studied, low-intensity behaviours predominate. In large macropods, the most frequent behaviour is lying down or standing still, followed by slow locomotion and only very occasional hopping [3,5]. Basic data on macropod locomotor activity exist in addition to that of the tammar wallaby already mentioned [6]. Locomotion comprises only 5–10% of the behavioural repertoire of the parma wallaby (Macropus parma) [62]. The larger red and grey kangaroos (Macropus rufus and M. giganteus) spend the day alternating between lying, standing, crouching, grazing and licking [5]. Agile wallabies' (Macropus agilis) most common behaviour is foraging (73%), followed by ‘vigilance’ (23%) and locomotion (0–6%) [63]. In grey kangaroos, over 90% of daily activity is crouching and lying, with only 0.0–3.3% accounted for by ‘moving’ [64]. Artiodactyls are similarly slow most of the time: wildebeest (Connochaetes spp.) travel only 2–3 km daily [65]; red deer (Cervus elaphus) move on average 100–400 m h−1 [66], while giraffe may walk for 5 h daily [67] and can canter for only a few minutes at a time [68].

A scaling trend in gait preference might exist, such that small animals hop, trot or gallop more frequently, or display different relative loading among forelimbs, hindlimbs and tail, than large animals, which could influence the interpretation of our results. Macropod species that live in open country generally have a shorter period of suspension than those that live in dense forests or rocky hills, with a potential phylogenetic contribution to duty factor [2] and thus peak ground reaction forces and bone strains. Our phylogenetically corrected scaling analysis found only limited effects of phylogeny on skeletal scaling parameters, suggesting little relationship between behavioural ecology, locomotor style and bone geometry scaling within macropod and artiodactyl clades. Comprehensive behavioural ecology, activity pattern, kinematic, ground reaction force and tissue strain data would help to place the skeletal scaling that we have identified into the context of functional loading. Examining the skeletal scaling of closely related quadrupeds from the diprotodont order such as wombats (Vombatidae), koala (Phascolarctos) and possums (Phalangeridae) in the context of hopping macropods' skeletal scaling may help to further separate phylogenetic from mechanobiological effects. Kinematic data exist for sheep, goats [69,70], pigs [71], cattle and a small number of other artiodactyls during walking [72] and for a small number of macropods [73-75], but bone strain data are missing in all but a few species [70,76].

The lack of differential cross-sectional scaling in the macropod hindlimb despite their hopping behaviour led us to the speculation that their bones might have enhanced fatigue damage repair by increased remodelling, thus reducing the need for extra bone mass. We failed to find secondary osteonal remodelling in a Macropus giganteus femur sample, which was somewhat unusual for an animal of 33 kg body mass [77]. The absence of secondary osteons may relate to the single sample failing to include any by chance, a load-related suppression of remodelling protecting bone from local weakening due to osteoclastic resorption [59], or infilling of existing osteons as occurs in horses after moderate-intensity training [78]. The current and other studies of bone organ allometry assume no size-related variation in bone microstructure or physiology and that all mammalian bone has similar biomechanical and mechanobiological behaviour. Our recent work demonstrated that secondary osteons are wider in larger animals and narrower in smaller animals [77], and that trabeculae are thicker and more widely spaced in larger animals [79], indicating that biophysical constraints or cellular behaviour may vary among mammals and potentially interact with whole-bone level scaling. Integration of macro- and micro-level perspectives in future scaling studies could be particularly informative.

Forelimb–hindlimb and bipedal–quadrupedal comparisons of scaling relationships have revealed very similar cross-sectional scaling against mass in the primary weight-bearing limb bones in artiodactyls and macropods, despite differences in their high-intensity gaits, suggesting that habitual low loads rather than occasional high loads may be the dominant stimuli for bone modelling (i.e. scaling) in individuals and across evolution. Cross-sectional scaling against length meanwhile appears to relate to clade-related specializations such as macropods’ long, gracile forelimb used in low-speed weight-bearing and grasping food, and artiodactyls' more robust forelimb bones specialized for cursorial locomotion. Our conclusions that high-speed bipedal hopping and quadrupedal galloping do not seem to leave a strong signal in limb bone geometry are some cause for caution in the interpretation of extreme behaviours in other, especially extinct, taxa from bone scaling. However, identification of locomotor patterns such as bipedalism/quadrupedalism is best done via multiple lines of evidence, which could still include bone scaling analysis [80,81].