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

The orientation of vascular canals in primary bone may reflect differences in growth rate and/or adaptation to biomechanical loads. Previous studies link specific canal orientations to bone growth rates, but results between different taxa are contradictory. Circumferential vascular canals (forming laminar bone) have been hypothesized to reflect either (or both) rapid growth rate or locomotion-induced torsional loading. Previous work on the hindlimb biomechanics in the emu shows that the femur and tibiotarsus experience large shear strains, likely resulting from torsional loads that increase through ontogeny. Here, we test how growth rate and biomechanical loading affect bone laminarity in wing and hindlimb elements from growing emu (2–60 wks). If laminar bone is an adaptation to torsion-induced shear strains, it should increase from juveniles to adults. Alternatively, if bone laminarity reflects rapid growth, as has been shown previously in emu, it should be abundant in fast-growing juveniles and decrease with age. Transverse mid-shaft histological sections from the limb bones (femur, tibiotarsus, humerus, ulna, and radius) were prepared and imaged. Growth rates were measured using fluorescent bone labels. Vascular canal orientation was quantified using laminarity index (proportion of circumferential canals). Principal components analysis was performed to convert highly correlated variables (i.e., mass, age, growth rate, and shear strain) into principal components. Random-intercept beta regression modeling determined which principal components best explained laminarity. The fastest growth rates were found in young individuals for all five skeletal elements. Maximum growth rate did not coincide with peak laminarity. Instead, in the femur and tibiotarsus, elevated laminarity is strongly correlated with adult features such as large size, old age, and modest growth rate. This result is contrary to predictions made based on a previous study of emu but is consistent with results observed in some other avian species (penguin, chicken). Shear strain in the caudal octant of the femur and tibiotarsus is positively correlated with laminarity but has a weaker effect on laminarity relative to mass, age, and growth rate. Laminarity in the wing elements is variable and does not correlate with ontogenetic factors (including mass, age, and growth rate). Its presence may relate to relaxed developmental canalization or a retained ancestral feature. In conclusion, ontogeny (including growth rate) is the dominant influence on vascular canal orientation at least in the hindlimb of the emu.


INTRODUCTION
Avian bone tissue is highly vascularized with a fibrolamellar structure, which allows for rapid growth by depositing randomly arranged spicules of woven bone initially, followed by in-filling of the cancellous spaces with centripetal lamellar bone, forming primary osteons (Francillon-Vieillot et al., 1990;de Ricqlès et al., 1991;Curry, 2002). Each primary osteon contains a central canal that houses blood vessels and nerves. These vascular canals vary in orientation and bones can be classified based on the predominant canal orientation. Laminar bone has a higher proportion of canals with circumferential orientation (parallel to the periosteal surface of the bone) relative to other orientations. Additional canal orientations include: radial, those orthogonal to the periosteal surface; longitudinal, those running parallel to the long axis of the bone; and oblique, all other orientations (de Ricqlès et al., 1991).
It has been hypothesized that differences in primary vascular canal orientation might be a reflection of growth rate, biomechanical loads, or phylogenetic relationships (Padian, 2013). Amprino (1947) first suggested that the organization of bone microstructure may be influenced by bone growth rate, such that woven bone is deposited during rapid growth and lamellar bone during slow growth (de Ricqlès et al., 1991). Further studies have investigated whether specific primary vascular canal orientations in fibrolamellar bone are also associated with slow or fast growth by directly comparing microstructure with bone growth rates measured through the use of injectable fluorochromes (Castanet et al., 2000;de Margerie, Cubo & Castanet, 2002;de Margerie et al., 2004). Rapidly growing hindlimb bones of ratites have been found to exhibit structure that is laminar and reticular (bone with numerous obliquely-oriented canals), whereas the more modest-growing wing elements of ratites exhibit reticular and longitudinal canal structure (Castanet et al., 2000). This suggests that laminar bone, in part, may reflect faster growth rates. This result was supported in a recent study of pigeon wing elements, which showed that peak laminarity (proportion of laminar bone) coincides roughly with the growth spurt in each element (Ourfalian, Ezell & Lee, 2016). However, work on mallard long bones showed no relationship between growth rate and predominant vascular canal orientation (de Margerie, Cubo & Castanet, 2002). Additionally, in the king penguin, radially-oriented canals dominated in the fastest growing sections, not circumferential canals (laminar bone) (de Margerie et al., 2004). Likewise, chickens selected for fast growth showed limb bones with predominantly radial canals (Williams et al., 2004;. Laminar bone has been hypothesized to better resist torsional loading. In laminar bone, the bone tissue is arranged in 'sheets' or 'plates' between layers of circumferential canals. Shear strain is thought to flow continuously within these 'sheets', and thus the concentrated stresses on the bone tissue surrounding the canals is reduced (de Margerie et al., 2004). Indeed, bone elements that are hypothesized to predominantly experience torsional loads have been found to exhibit laminar bone. Laminar bone is found to be most abundant in the humerus, ulna, and femur in a large sample of flighted bird species (de Margerie, 2002;de Margerie et al., 2005). In vivo strain gauge studies have shown that these elements experience predominantly torsional loads in at least some species: the humerus in the pigeon during flapping flight, the ulna in the turkey during wing flapping while on the ground, and the femur in the chicken and emu during terrestrial running, (Rubin & Lanyon, 1985;Biewener & Dial, 1995;Carrano & Biewener, 1999;Main & Biewener, 2007).
A limitation of previous studies of laminar bone is the indirect comparison of bone histology in one species with bone growth rates and/or in vivo strain gauge measures taken from different species. In this study, we present an analysis of laminar bone in a species in which bone growth rate and in vivo bone strain data were directly measured. The emu (Dromaius novaehollandiae, Order Struthioniformes, Family Dromaiidae) is a flightless bird endemic to Australia, but widely farmed in the US. The individuals included in this study comprise a growth series that were previously injected with fluorescent bone labels and surgically implanted with gauges to measure in vivo locomotor strains in the femur and tibiotarsus (Main & Biewener, 2007). Shear strains, produced by torsional loads about the long axes of the bones, were the predominant type of strain in the two bones and increased from juveniles to adults (Main & Biewener, 2007). Bone strains were not measured in the wing elements of these individuals. Presumably, shear strains are negligible in the wing elements because emus have extremely reduced wings, which have no known function other than occasionally being raised to aid thermoregulation (del Hoyo, Elliot & Sargatal, 1992;Maxwell & Larsson, 2007).
Therefore, if laminarity is an adaptation to torsion-induced shear strains, we predict that hindlimb bone laminarity will increase from juveniles to adults. Consistent with this hypothesis, we expect no trend in laminarity in the vestigial wing elements from juveniles to adults. Alternatively, if bone laminarity reflects rapid growth (as has been shown previously in emu), it should be abundant in juveniles and decrease with age as growth slows in adults for all elements. Incorporating growth rate measurements, direct biomechanical data, and direct histological classification of laminarity makes this study a first of its kind that will be able to clarify the importance of growth and mechanics on vascular canal orientation in emu limbs.

MATERIALS & METHODS
This study samples forelimb and hindlimb elements from eight emus ranging in age from 2 to 60 weeks (Table 1). Birds used in this study were euthanized as part of a previous study (Main & Biewener, 2007)  female birds were included based on availability. Emus exhibit a minor degree of sexual dimorphism, with females being slightly larger on average (del Hoyo, Elliot & Sargatal, 1992). The difference in size is not large enough to be considered a confounding factor for this study. As a part of a previous study (Main & Biewener, 2007), each bird was given a single intramuscular injection of xylenol orange (80 mg/kg) followed by calcein (30 mg/kg). Injections were given one week apart in birds less than 16 weeks of age and two weeks apart in birds between 16 and 65 weeks of age. Xylenol and calcein are fluorescent labels that are incorporated rapidly into newly mineralizing surfaces of bone at the time of injection (An & Martin, 2003). Thus, the time elapsed and the space between xylenol and calcein labels allows the calculation of periosteal (radial) growth rate. One week after the last injection, surgery was performed to attach strain gauges to the cranial, caudal, and lateral aspects of the left femur and the cranial, caudal, and medial aspects of the left tibiotarsus. Single element strain gauges were used on the lateral femur and cranial and medial tibiotarsus. Rectangular rosette gauges were used on the cranial and caudal femur and caudal tibiotarsus. Rosette strain gauges allow both tensile and compressive principal strains and their orientations to be measured, and were placed so the central element of the gauge was parallel to the long axis of the bone. One day after surgery, the birds were run on a treadmill over a wide range of speeds and gaits. The raw data produced from the strain gauges were converted from voltage to microstrain using a custom MATLAB program. Shear strains were calculated from the rosette strain gauges using standard equations (Biewener & Dial, 1995). High quality shear strain data were most consistently collected from the caudal cortices of the femur and tibiotarsus and that is what is reported here. Trials in which the birds ran with a duty factor near 0.50 are included in the shear strain analysis (mean ± SD: 0.50 ± 0.02). Duty factor is the proportion of the time that the animal's foot spends on the ground during a stride relative to the entire stride time. A duty factor of 0.50 represents the point at which the birds transition to a running gait that incorporates an aerial phase. This is a relatively slow run for emu, but represents the duty factor for which we could maximize the number of animals included in this study based upon successful strain gauge function. Each trial was represented by five footfalls and, generally, two trials were collected for each bird. Following bone strain data collection, animals were euthanized. After death, whole wings were removed from the individuals and stored frozen. Histological sections of the femora and tibiotarsi were prepared (see (Main & Biewener, 2007 for details) and shipped with the frozen wings to Midwestern University.

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 On fluorescent images, growth rate was measured by taking the mean distance (white arrow) between the periosteal extent of the xylenol (red) and calcein (green/yellow) tags divided by number of days between injections. (D) Laminarity was measured by approximating each in-focus primary vascular canal with the best-fitting ellipse, using the best-fitting ellipse to categorize canal orientations, and calculating the proportion of circumferential canals. (Tibiotarsus 17, 4.6 wks).
Full-size DOI: 10.7717/peerj.7616/ fig-2 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, Table 2 Growth rate and bone laminarity for emu specimens. Laminarity Index (LI) was measured in a sample area outlined by the periosteal extent of two bone fluorochromes. Forelimb LI was rescaled following (Smithson & Verkuilen, 2006) to prepare the data with 0 values for beta regression.

Specimen
Age ( Table 3 Caudal shear strain and caudal octant laminarity for emu specimens. Caudal shear strain data were previously collected by Main & Biewener (2007). Caudal laminarity index (LI) was measured in a sample area outlined by the periosteal extent of two bone fluorochromes in the caudal octant only. LI values were rescaled following (Smithson & Verkuilen, 2006) to prepare data with 0 values for beta regression. Only specimens for which caudal shear strain data were available are included.

Specimen
Age ( , 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: 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  Robust principal components analysis of the cardinal octants from the femur and tibiotarsus. An ''ontogenetic axis'' (PC 1) accounts for 95% of the variance with juvenile features to the left (small size, young age, and rapid growth) and adult features to the right (large size, old age, and modest growth). Residual variation in growth rate is absorbed into PC 2, which accounts for another 4.7% of the cumulative variance. Full-size DOI: 10.7717/peerj.7616/ fig-4 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 −5 ; Fig. 5).

Figure 5 Effect of the ''ontogenetic axis'' (PC 1) on laminarity (LI) in the cardinal octants of the hindlimb.
Random-intercept beta regression reveals that the ''ontogenetic axis'' (PC 1) accounts for 89% of the variation in laminarity from the cardinal octants of femur (green triangles) and tibiotarsus (blue squares). Elevated laminarity values are strongly correlated with adult features such as large size, old age, and modest growth rate (p < 5.51e −5 ).
Full-size DOI: 10.7717/peerj.7616/ 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. Random-intercept beta regression reveals that the ''ontogenetic axis'' (PC 1) accounts for 78% of the variation in laminarity from the caudal octants of femur (green triangles) and tibiotarsus (blue squares). Laminarity generally increases along that axis (p = 0.013) to which shear strain has a minor contribution. radius. An ''ontogenetic axis'' (PC 1) accounts for 97% of the variance with juvenile features to the left (small size, young age, and rapid growth rate) and adult features to the right (large size, old age, and modest growth rate). Residual variation from growth rate largely contributes to an ''ontogeny-independent growth rate axis'' along PC 2. This axis only explains 2.3% of the cumulative variance. Full-size DOI: 10.7717/peerj.7616/ fig-8 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 threedimensional vascular canal orientation in the humerus and femur of growth-controlled broiler chickens also found elevated laminarity in a slow-growing (feed-restricted) group . 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;. Instead, they have bones that are poorly-vascularized with a parallelfibered matrix. Because that histology tends to grow very slowly in other species (e.g., <five µm d −1 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 threedimensional 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;. 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 , 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.

CONCLUSIONS
In the emu limb skeleton, ontogenetic factors such as size, age, and growth rate have major effects on vascular canal orientation. The effect of shear strain is relatively weak and suggests that laminar bone is not a good predictor of torsional loading. Even though the forelimb elements likely experience minimal loading, the humerus and ulna show wide variation in laminarity, perhaps due to relaxed selection. Future studies should investigate laminarity in other palaeognathous birds to better understand the effects phylogeny, ontogeny, and torsional loading have on bone laminarity. Other future work should focus on the experimental manipulation of biomechanical loads to observe the effects on vascular canal orientation in limb bones and to better understand to what extent torsional load influences the development of limb bone laminarity. It is also important that variation found between different populations be addressed and studied further. Emu body mass growth rates vary among populations (e.g., Goonewardene et al., 2003), but it is unknown to what extent laminar bone also varies with environmental, dietary, or genetic factors. This study has shown that in emu limb bones, laminarity reflects a complex interplay of ontogeny and biomechanical loads.