The effects of growth rate and biomechanical loading on bone laminarity within the emu skeleton

Histological preparation

Emu wings were thawed and feathers, skin, muscles, and tendons were reflected to expose the skeletal elements. Both right and left wings were used based on availability. Using digital calipers, total length of each bone was measured and recorded. A 37-mm segment was removed using a Dremel tool from the mid-shaft region of the humerus, ulna, and radius. For two and four week old individuals whole elements were harvested due to their small size. Segments were labeled with permanent marker to maintain orientation. Dissected bone segments were placed in 10% neutral buffered formalin for fixation and then dehydrated in a graded ethanol series (70%, 85%, 100%) under vacuum. Segments were cleared with a xylene-substitute (Histo-clear; National Diagnostics, Atlanta, Georgia, USA). The bone segments were then vacuum- infiltrated and embedded in glass vials using Osteo-Bed Plus Resin, a two-part methyl methacrylate (Polysciences Inc.). Vials were placed in a 32 °C bead bath to fully harden.

Once the resin hardened, vials were broken and two roughly 800-µm transverse sections were cut using a diamond blade saw (Isomet 1000; Buehler, Lake Bluff, Illinois, USA). These sections were attached to frosted glass slides using two-ton epoxy (Devcon, Milpitas, California, USA), keeping consistent spatial orientation. Slides were then ground to a thickness of 100 ± 10 µm using a graded scale of grit paper on a stand grinder (Metaserv 250; Buehler, Lake Bluff, Illinois, USA) and coverslipped with Permount (Fisher Scientific). The histological preparation was modified from An & Martin (2003) and closely followed Lee & Simons (2015).

Image collection

The undecalcified sections contain xylenol (orange) and calcein (green) fluorochromes that were incorporated into newly mineralizing bone at the time of injection (see above for injection schedule). These fluorochromes create stable long lasting tags (Van Gaalen et al., 2010) and were examined under bright-field and fluorescent illumination with a motorized epifluorescent microscope (IX73; Olympus). The xylenol (orange) and calcein (green) tags were revealed using TRITC and FITC filter cubes, respectively, and a multichannel (red, green, bright-field) image of each section was generated with imaging software (cellSens, Olympus). Sufficient optical resolution (10X UPlanAPO ≈0.84 µm; 20X UPlan S-APO ≈0.45 µm) allowed a dual color-monochrome camera (DP80, Olympus) to capture high quality images (10X = 1.02 µm/pixel; 20X = 0.51 µm/pixel).

Calculating bone laminarity and radial growth rates

Bright-field and fluorescent images were obtained from the wing and hindlimb elements (Figs. 1 and 2) and divided into equal octants from the estimated bone centroid. Four octants representing the cardinal anatomical positions (wing elements: cranial, caudal, dorsal, ventral; hindlimb elements: cranial, caudal, lateral, medial) were extracted (Fig. 2). Using ImageJ, each extracted octant was then uncurved using the “Straighten” function. The purpose of straightening was to standardize the periosteal tangent line so that appropriate measurements could be made in classifying the orientation of the vascular canals (Lee & Simons, 2015). To ensure there was minimal deformation of the image during the straightening process, known test angles were placed upon the image and measured in relation to the periosteal surface after the straightening function had been applied. Only those images with an average deformation less than or equal to 10° were accepted.

Within each of the four octants, the calcein green and xylenol orange tags were outlined with two reference lines. The distance between reference lines was measured at 10 equally spaced points in each octant. Growth rate was measured by taking the mean distance between consecutive fluorescent tags divided by number of days between injections (Fig. 2C).

Degree of laminarity (Laminarity Index, LI) was measured in the interval of bone between the fluorochrome reference lines across all four octants. Using ImageJ, an ellipse was drawn within each in-focus primary vascular canal in the measurement interval (Fig. 2D). Branching canals were separated at branch points and counted individually. Sharply curving canals were treated as branching. The aspect ratio and angle at which the ellipse sat in relation to the straightened periosteal surface was measured. We used the criteria set forth by de Margerie (2002) to classify the orientation of the vascular canals: (1) “circular” (circumferential) canals are oriented parallel (0° ± 22.5°) to the periosteal surface of the bone; (2) radial canals are orthogonal (90 ± 22.5°) to the periosteal surface; (3) longitudinal canals run parallel to the long axis of the bone and have ellipses with an aspect ratio of less than 3; (4) oblique canals are all other orientations. Only primary vascular canals were measured. Secondary osteons in the sample area were excluded. We used a simple proportion (number of circumferential canals to the total number of canals) to quantify laminarity. To test the growth hypothesis, we used the laminarity index calculated from all sampled octants (Table 2, Table S1). Because consistent high quality shear strain data are only available from the caudal octant, we used the laminarity index calculated from only the caudal octant to test the mechanical hypothesis (Table 3, Table S2).

Robust Principal Components Analysis (RPCA) and beta regression

The explanatory variables thought to affect laminarity in the emu show multicollinearity. For example, mass and age covary with each other (Goonewardene et al., 2003) as well as with growth rate (Montes et al., 2005), and shear strain (Main & Biewener, 2007). If left unaddressed, multicollinearity can decrease precision and reliability when estimating the effect of one variable while holding the others constant (Fekedulegn et al., 2002). Principal components analysis (PCA) accounts for this multicollinearity by forming new uncorrelated variables (i.e., principal components) that are linear combinations of the original explanatory variables while preserving as much of the original variation as possible (Hammer, Harper & Ryan, 2001). However, PCA is highly sensitive to variables with large variances and skewed distributions (Hubert, Rousseeuw & Verdonck, 2009), so we standardized (i.e., centered by the median and scaled by the median absolute deviation) mass, age, growth rate, and shear strain with the function “RobScale” (Signorell, 2019) in R (R Development Team, 2019). This process stabilizes variance and minimizes the effect of absolute scale in the calculation of principal components. Skewed data are often transformed prior to PCA (e.g., logarithmic or Box-Cox), but such transformations may worsen skewness or complicate PCA interpretation (Hubert, Rousseeuw & Verdonck, 2009). Instead, we performed robust PCA (Hubert, Rousseeuw & Verdonck, 2009) as implemented by the R package “rospca” (Reynkens, 2018), which is suitable for skewed data. Three datasets were analyzed separately: (1) cardinal octants from hindlimb elements; (2) caudal octants from hindlimb elements; and (3) cardinal octants from forelimb elements. The results of each robust PCA are presented in Table 4.

For each dataset, the minimum number of principal components (PCs) was selected to cover approximately 95% of the observed variance of the original explanatory variables. We assessed the relationship between PC(s) and mean laminarity index (LI) using beta regression as implemented by the R package “gamlss” (Rigby & Stasinopoulos, 2005). This method is appropriate when the response variable (LI) is a proportion (Warton & Hui, 2011). To accommodate values of 0 in the caudal hindlimb and cardinal forelimb datasets, LI values were rescaled to the effective interval of [0.005, 0.995] (Smithson & Verkuilen, 2006). The logit link function was used to connect mean LI to a linear combination of the PCs. Pseudoreplication is a concern because different bone elements were sampled from the same individual (Hurlbert, 1984; Gillies et al., 2006; Lee & O’Connor, 2013; Jordan, 2018), so we combined the logit link function with a random-intercept model as follows: (1)${\displaystyle {\mu }_{logit\left(LI\right)}={\beta }_0+{\beta }_{1}P{C}_1+{\beta }_{2}P{C}_2+{\beta }_{3}Element+\gamma }$ where PC is the principal component, β is regression coefficient, Element is a dummy variable coding for element type, γ is the random-intercept effect of “specimen ID”, and μ_logit(LI) is the logit link function for the mean of LI (Ferrari & Cribari-Neto, 2004).

For each cardinal dataset, we evaluated two models. The first model includes PC 1 as the sole predictor given that it accounted for approximately 95% of the variance in the original explanatory variables. The second model adds element type as a dummy variable. For the caudal dataset, PC 1 and PC 2 covered at least 95% of the variance in the original explanatory variables. Therefore, we evaluated six models. The first three models involve PC 1 and PC 2 individually as sole predictors as well as together in additive combination. The remaining three models add element type as a dummy variable (Table 5). The small-sample correction of Akaike’s Information Criterion (AIC_c) (Hurvich & Tsai, 1989) was used to compare models within each dataset. In general, the best supported model has the lowest AIC_c value (Burnham & Anderson, 2002). Relative support between the best and alternative models was assessed with difference (ΔAIC_c) values. Alternative models with ΔAIC_c values greater than 3, which is equivalent to a p-value of 0.051 (Taper, 2004), were rejected as having weak support. Raw data and R script for analyses can be found in the Supplementary Files (Tables S1, S2, Code S1).

Regression analysis of cardinal octants from the femur and tibiotarsus

Principal component (PC) 1 (eigenvalue = 8.181), consisting of mass, age, and growth rate, accounts for 95% of the cumulative variance (Table 4). Mass and age loadings have the same sign, whereas growth rate loading has an opposite sign. The loadings suggest that PC 1 represents an “ontogenetic axis” with juvenile features (small size, young age, and rapid growth rate) at one end and adult features (large size, old age, and modest growth rate) at the other (Fig. 4). PC 2 (eigenvalue = 0.403) consists of the residual variation in (i.e., “ontogeny-independent”) growth rate and accounts for 4.7% of the cumulative variance (Table 4).

The random-intercept beta regression model without element type as a predictor has overwhelming support (Table 5). It predicts that ∼89% of the variation in laminarity is explained by the “ontogenetic axis” (Table 5). Juvenile features (e.g., small size, young age, and rapid growth rate) are correlated with lower laminarity values, whereas adult features are correlated with higher laminarity values (p < 5.51e⁻⁵; Fig. 5).

Regression analysis of caudal octants from the femur and tibiotarsus

Table 4 shows that the first two PCs account for at least 95% of the cumulative variance—88% by PC 1 (eigenvalue = 9.435) and 9% by PC 2 (eigenvalue = 0.989). Similar to cardinal octant data, mass, age, and growth rate contribute strongly to PC 1. Their loadings are also consistent with variation along an “ontogenetic axis” with juvenile features at one end (small size, young age, and rapid growth rate) and adult features at the other (large size, old age, and modest growth rate). Although strain has a minor contribution to PC 1, it dominates PC 2, which we interpret as a “loading effect axis” (Fig. 6).

The model with PC 1 as the sole predictor of caudal octant laminarity has strongest support based on ΔAIC_c and explains 78% of the variation in caudal octant laminarity (Table 5). Caudal and cardinal datasets show slightly different estimates for the coefficient of PC 1, which is consistent with slight inter-octant variation in laminarity. Variation aside, the overall ontogenetic trend is similar: higher laminarity values are correlated with adult features, whereas lower laminarity values are correlated with juvenile features (Fig. 7; p = 0.013). Although shear strain contributes to PC 1 and generally increases along the “ontogenetic axis”, the effect is relatively weak. To highlight this, we multiplied the regression coefficient and eigenvectors of PC 1. The resulting standardized coefficients of the original predictors (growth rate = − 0.041; mass = 0.097; age = 0.150; strain = 0.035) suggest that the relative effects of growth rate, mass, and age on caudal octant laminarity in the hindlimb are 1.2 to 4.3 times greater than that of strain.

The caudal hindlimb dataset does not support the “loading effect axis” (PC 2) as a strong predictor for laminarity. Compared to model 1 (PC 1 as the sole predictor), the remaining alternative models with PC 2 only explain an additional 5–6% of the variation in caudal hindlimb laminarity. Each of these alternative models has ΔAIC_C larger than 3 indicating weak to no support (p > 0.05). Furthermore, the model coefficient for PC 2 in model 2 (sole predictor) and model 3 (additive combination with PC 1) is not significant (p-value equals 0.169 and 0.056, respectively; Table 5).

Regression analysis of cardinal octants from the humerus, ulna, and radius

In the absence of strain data for the forelimb elements, mass, age, and growth rate dominate PC 1 (eigenvalue = 7.977), and this “ontogenetic axis” accounts for 97% of the cumulative variance (Table 4, Fig. 8). PC 2 (eigenvalue = 0.187) absorbs residual variation in growth rate and accounts for an additional 2.3% of the cumulative variance.

Although the model with “ontogenetic axis” (PC 1) as the sole predictor has the best support, it only accounts for 27.7% of the variation in forelimb laminarity. Moreover, the model coefficient for PC 1 is not significant (p = 0.089; Table 5).

Does bone laminarity reflect fast growth?

The highest periosteal growth rate in all elements was found in the femur of the 4.6-week old individual (Table 2). As expected, hindlimb elements had higher growth rates than forelimb elements, reaching a maximum of 163 µm/day in the femur and 99 µm/day in the tibiotarsus. The humerus grew the fastest of the wing elements, reaching a maximum rate of 25 µm/day measured in the 2.3 and 4.6 week old individuals. Birds older than 8 weeks experienced a drastic decrease in bone growth rate in both hindlimb and forelimb elements. Previous analysis of emu somatic growth rate (increase in body mass) showed the maximum rate of growth (inflection point) to be about 15–17 weeks of age (Goonewardene et al., 2003). Our results reveal that age at maximum bone growth (approximately 5 weeks) precedes the somatic growth inflection, similar to other vertebrates (Lee & O’Connor, 2013). Therefore, caution is warranted when inferring somatic life-history milestones, such as growth spurts, solely from skeletal data.

Principal components analysis reveals a large proportion of variance lies along an “ontogenetic axis” (Figs. 4, 6 and 8). One end of the axis is represented by juvenile traits such as small size, young age, and rapid growth, whereas adult traits such as large size, old age, and modest growth characterize the other end. This “ontogenetic axis” has a significant influence on laminarity in the femur and tibiotarsus (Figs. 5 and 7), whether analyzed in cardinal octants (p < 0.001) a single octant (p = 0.013). Elevated laminarity in the hindlimb appears correlated with adult features, including modest growth rate. This relationship is consistent with findings in the king penguin that also reported laminar bone to be associated with modest growth rates in four limb bones: femur, tibiotarsus, humerus, and radius (de Margerie et al., 2004). More recently, a study using microCT to assess three-dimensional vascular canal orientation in the humerus and femur of growth-controlled broiler chickens also found elevated laminarity in a slow-growing (feed-restricted) group (Pratt & Cooper, 2018). Interestingly, the effect of the “ontogenetic axis” on laminarity in the wing elements is weak (p = 0.089), only explaining 27.7% of variation (Fig. 9). The weakened “ontogenetic axis” in wing bone laminarity is consistent with relaxed selection on the vestigial wing leading to increased anatomical variability in the species (Maxwell & Larsson, 2007).

Notably, our results differ from those previously reported for young emu bones in which laminar and reticular bone was found in the fastest growing hindlimb bones (Castanet et al., 2000). In particular, Castanet et al. (2000) found laminar bone to be most abundant in the femur and tibiotarsus of emu less than 2 months of age, which corresponds to the youngest individuals in our study. Based on reported body masses, the emus included in our study were about 2–3 times heavier than the emus in the Castanet et al. (2000) study for a given age (Table 1). The reason for the differences in size and ontogenetic patterns for bone vascularity types between these two emu samples remain unknown, but could be related to genetic, dietary, or rearing conditions between the two groups. If laminarity is associated with lower growth rates, the youngest emus we studied may have been growing too fast for laminar bone to form. The highest growth rate measured was in the femur of the 4.6-week-old bird (162.62 µm/day), which was about twice the highest growth rate found in the femoral reticular bone tissue reported in the prior study (89.4 µm/day). Our study did not specifically address reticular bone, but by taking the proportion of oblique vascular canals (a “reticular index”), we found the amount of reticular bone in the fastest growing individual to be low in the hindlimb elements (femur and tibiotarsus: 0.17), and moderate to high in the wing elements (humerus: 0.62, ulna: 0.58, radius: 0.45). This result is, at least, consistent with the previous study because Castanet et al. (2000) found reticular bone to be more abundant in the humerus. The data analysis in the earlier emu study (Castanet et al., 2000) was conducted before a more rigorous method for quantifying canal orientation was developed (de Margerie, 2002), which may contribute to the observed differences between the two separate emu populations. Regardless, the results of the present study do not support the hypothesis that elevated laminarity reflects rapid growth in the emu and instead link elevated laminarity with more modest growth rates.

Does bone laminarity reflect biomechanical load?

Robust PCA is consistent with previous scaling analysis showing that shear strains in the femur and tibiotarsus of the emu increase with growth (Main & Biewener, 2007). Shear strain covaries with mass, age, and growth rate along an “ontogenetic axis” such that juveniles experience smaller shear strains in the caudal octants, whereas adults experience larger strains (Table 4; Fig. 6). However, shear strain contributes comparatively less to the “ontogenetic axis” than the other covariates. Consequently, it has a correspondingly minor effect on caudal octant laminarity in the hindlimb. Indeed, transformation of the effect of the “ontogenetic axis” on laminarity back into the scale of the actual covariates reveals that the standardized effects of age (0.150), mass (0.097), and growth rate (−0.041) are relatively larger than that of shear strain (0.035). Elevated laminar bone tissue, in combination with increased bone mineralization and decreased bone curvature during growth, may have collectively helped mitigate shear strains despite the large increase in mass (Main & Biewener, 2007).

Residual variation in shear strain that is not accounted for by the “ontogenetic axis” forms a “loading effect axis” (Table 4; Fig. 6). However, this effect does not correlate with laminarity either as the sole predicator (p = 0.169) or in additive combination with ontogeny (p = 0.056). Put together, the results clearly show that ontogenetic factors largely influence the formation of laminar bone in the caudal hindlimb of the emu, although torsion-induced shear strain is a minor additional factor.

The weak association between laminarity and shear strain limits the predictive potential of this relationship. Our results for emu hindlimb bones are consistent with previous studies of other limb elements presumably loaded in torsion. For example, when comparing wings of similar shape, laminarity in wing bones can be similar despite differences in presumed biomechanical load associated with unique primary flight modes (Simons & O’Connor, 2012; Marelli & Simons, 2014). In addition, preferred flight mode may only have subtle effects on overall loading of the bones, with the dominant loads being the high strains present during take-off (Biewener & Dial, 1995). Furthermore, despite sharing with birds convergent features related to powered-flight such as torsionally loaded bone with relatively thin cortical walls (Swartz, Bennett & Carrier, 1992), bone mineral density (Dumont, 2010), and metabolic rate (e.g., Maina, 2000), bats lack laminar bone entirely (Lee & Simons, 2015; Pratt et al., 2018). Instead, they have bones that are poorly-vascularized with a parallel-fibered matrix. Because that histology tends to grow very slowly in other species (e.g., <five µm d⁻¹ de Margerie, Cubo & Castanet, 2002; Castanet et al., 2004), Lee & Simons (2015) speculated that bats do not grow fast enough to form laminar bone. Indeed, somatic growth is approximately four times slower in bats than in birds of comparable size (Lee & Simons, 2015). Caution is warranted, however, as we show in the emu that interchanging somatic and skeletal growth may be misleading. Whatever the actual cause is, the evidence is clear that elevated laminarity is not a prescriptive feature of torsionally loaded bone.

In this study, age groups are represented by one individual, with the exceptions of the youngest and oldest age groups that contain two. A larger sample of individuals in each age group would allow for investigation of how individual variation may or may not affect the relationship between LI, shear strain, and growth rate. Laminarity indices can be quite variable among individuals in some species. The pigeon humerus, which has been shown to experience large torsional loads, has been documented to exhibit both high and low laminarity in different individuals (Lee & Simons, 2015; Ourfalian, Ezell & Lee, 2016; Skedros & Doutré, 2019). Similarly, a pooled sample of humeri from eight Red-tailed hawks show LI values that range from 0.30–0.70 (Simons & O’Connor, 2012; Marelli & Simons, 2014). Whether these variability patterns are biological or methodological is unclear. Laminarity measured on a histological section is a 2-dimensional representation of a 3-dimensional meshwork of vascular canals in cortical bone. This research is limited by the assumption that one or two closely placed mid-shaft histological sections are an accurate representation of vascular canal structure. MicroCT-based assessment of the three-dimensional network of vascular canals suggests that traditional 2D histological methods may overestimate LI, but also recognizes that these differences may be methodological (Pratt & Cooper, 2017; Pratt et al., 2018). Certainly, future studies should continue to use microCT to assess how well laminarity measured on histological sections represents actual biological structure. In addition, the torsional resistance in bones may more likely be linked to the specific orientation of another histological feature: collagen fibers. Collagen fiber orientation (CFO) has been shown to reflect principal strain distributions (Riggs, Lanyon & Boyde, 1993; Skedros & Hunt, 2004; Skedros, Hunt & Bloebaum, 2004; Skedros & Doutré, 2019). Analysis of CFO is beyond the scope of the current study. However, given the known positive correlation between transversely oriented collagen fibers and bone laminarity (de Margerie et al., 2005), we would expect a similar pattern for the femora and tibiotarsi examined here.

Although there is no direct biomechanical data for the forelimb elements of these birds, the wing elements presumably experience minimal loading. The emu wing is extremely reduced in size, even when compared to other ratites, and has almost no observed function (del Hoyo, Elliot & Sargatal, 1992). Wing muscles of emu contain primarily slow acting tonic muscle fibers that may not allow much wing movement (Maxwell & Larsson, 2007), which suggests the underlying wing elements would experience minimal biomechanical loading. Despite the assumption that the emu wing is under minimal load, a moderate to high degree of laminarity was found in at least the humerus and ulna (Table 2). This laminarity can be attributed to the modest bone growth rate observed in the wing elements and/or to the third factor affecting bone microstructure: phylogenetic relationships. Within the paleognaths, it has been hypothesized that at least three independent flight losses have occurred, with only one order (the tinamous) still retaining the ability to fly (Harshman, Braun & Braun, 2008; Mitchell et al., 2014). The moderate/high wing bone laminarity may be a feature of the flighted common ancestor of paleognaths that is retained in the flightless descendants. Indeed, significant phylogenetic signal has been found in some osteohistological features in a sample of paleognaths (Legendre et al., 2014). Future studies should investigate the histological and in vivo loading of the flighted relatives of emus to better understand the potential influence of phylogeny on bone laminarity.