Development of limb bone laminarity in the homing pigeon (Columba livia)

Sampling and histology

We acquired salvaged carcasses of 17 homing pigeons (Stromberg’s Chicks and Game Birds Unlimited; Pine River, MN, USA). The sample comprises a postnatal growth series of known mass (15–563 g) (Table S1). Although the age range of the sample is 0–9 weeks, the precise age of death for most individuals was not recorded. The right fore- and hindlimb of each individual was dissected to reveal the humerus, radius, ulna, femur, and tibiotarsus (Fig. 2A). The length of each element was measured (Tables S2–S6), and a 1-cm mid-diaphyseal block from each bone was excised using a rotary tool (Dremel 4000; Dremel, Mt. Prospect, IL, USA). We followed an established protocol for preparing plastic-embedded undecalcified bone (Lee & Simons, 2015). Specifically, two transverse 700-μm wafers were cut from each specimen at mid-diaphysis using a precision saw (Isomet 1000; Buehler, Lake Bluff, IL, USA). Wafers were mounted (Gorilla Epoxy; Gorilla Glue Inc., Cincinnati, OH, USA) to glass slides and thinned to 100 ± 10 μm using a grinder/polisher (Metaserv 250; Buehler, Lake Bluff, IL, USA).

Imaging

Sections were acid-etched and stained with toluidine blue (Eurell & Sterchi, 1994) to improve contrast of in-plane primary vascular canals (Fig. 2B). The stain also highlights secondary osteons (specifically cement lines) and their (Haversian) canals, which are traditionally excluded from measurements of laminarity (De Margerie, 2002). Whole section images pre- and post-staining were captured with a motorized microscope (Ni-U; Nikon, Tokyo, Japan) as previously described by Lee & Simons (2015). Once imaging was completed, each section was mounted (Permount; Fisher Scientific, Hampton, NH, USA) with a glass coverslip (#1; Fisher Scientific, Hampton, NH, USA) for preservation (Lee & Simons, 2015). Sections are housed in the Arizona Research Collection for Integrative Vertebrate Education and Study (ARCIVES) at Midwestern University.

Bone profiles

A bone profile was prepared from each montage following procedures described by Lee & Simons (2015). Specifically, montages were sharpened in Photoshop (CS5; Adobe, San Jose, CA, USA) with the “Unsharp Mask” filter (5 px). The area bounded between the periosteal and endosteal surfaces was filled with white to represent bone. The surrounding non-bone area as well as in-plane vascular canals and resorption spaces were filled with black. Bone profiles were exported to ImageJ (1.51d; National Institutes of Health, Bethesda, MD, USA) for further analysis. We measured the periosteal circumference and vascular porosity (ratio of in-plane primary and secondary vascular area to total cortical area) of each bone profile (Tables S2–S6). Montages and bone profiles can be viewed at the Paleohistology Repository (Lee & O’Connor, 2013) and downloaded at Dryad.

Ontogenetic scaling of polar section modulus

Bone profiles were imported into ImageJ, and the BoneJ plugin v1.4.1 (Doube et al., 2010) was used to calculate geometric properties such as second moment of area (I), polar moment of area (J), and polar section modulus (Z_p). These properties describe the distribution of material around the centroid of a given bone section and are inversely proportional to bending stress, overall (bending and torsional) stress, and maximum overall stress, respectively (Schoenau et al., 2001; Ruff, 2002; Habib & Ruff, 2008; Young, Fernández & Fleagle, 2010; Hedrick et al., 2020). In other words, a bone section with large I, J and Z_p experience less stress for a given load, giving it greater strength and rigidity to withstand larger loads before failure. Because those properties are closely related, we only present the Z_p data (Tables S2–S6), which are the most relevant to test the torsional resistance hypothesis.

Z_p is proportional to the maximum torsional stress that occurs at the outermost surface of a bone loaded in pure torsion (Craig, 2000). It is calculated by taking the ratio of J and the maximum distance between the centroid and the periosteal surface (r_max), where failure is most likely to occur. BoneJ calculates Z_p without assuming circular or elliptical geometry by applying the following general relationship (Eq. (1)) to pixel data in the actual bone profiles: (1) $$Z_p={\displaystyle \frac{J}{r_{max}}={\displaystyle \frac{\int r^{2}dA}{r_{max}},}}$$where J is the integral sum of the area of each pixel dA (2.36E−7 mm²) representing bone that is a distance r from the centroid, and r_max is the maximum distance between the centroid and the periosteal surface (Doube et al., 2010). Z_p is an appropriate proxy for torsional strength and rigidity when the cross section of a long bone is nearly circular (Craig, 2000). To test the suitability of this proxy to each cross section, we used BoneJ to calculate the aspect ratio (I_max/I_min), which equals 1 in a circular cross section. Values for the bone sections range from 1.03 to 1.80 (Tables S2–S6). When compared to a reference figure (Daegling, 2002), the values of I_max/I_min in our sample suggest errors in torsional rigidity less than 6.7%. Therefore, we find no major problem in using this proxy.

Using R (R Core Team, 2019), Z_p and body mass were log₁₀-transformed to linearize their relationship before performing separate Type I linear regression for each element. The scaling coefficients (b) and exponents (a) of the regression analysis are presented in Table 1. Because Z_p and body mass (M) are each proportional to length³, allometric scaling of the log-log model, log₁₀(Z_p) = log₁₀(b) + a log₁₀(M), was inferred if the 95% confidence interval of the scaling exponent (a) excluded the value of 1 (isometry). Note that we chose to use Type I instead of Type II (RMA) regression for two reasons. First, Z_p and body mass were measured precisely with low error. Second and more importantly, there is natural individual variation in Z_p for a given body mass. Such variation, whether influenced by intrinsic or extrinsic factors, is known to weaken Type II regression, which has a tendency to detect a steeper relationship than actually exists (Kilmer & Rodríguez, 2017).

Robust principal component analysis and beta regression

In this study, we had to address the issue of multicollinearity among our explanatory variables: mass, bone length, and Z_p. Principal component analysis (PCA) enables the formation of new, uncorrelated predictors (principal components) through linear combinations of the original variables. As such, we were able to resolve the issue of multicollinearity while still being able to assess the effect of each variable (Fekedulegn et al., 2002). PCA, however, is known to be highly sensitive to non-normal data. Therefore, we used robust PCA, which is appropriate for skewed data (Hubert, Rousseeuw & Verdonck, 2009), as implemented by the R package “rospca” (Reynkens, 2018). We standardized mass, bone length and Z_p by median and median absolute deviation with the function “RobScale” (Signorell, 2019) in R. Data were grouped by homologous element, and a separate robust PCA was performed for each element (Table 2).

We used regression analysis to relate the minimum number of principal components (PCs) that account for at least 95% of the variation in the original variables with laminarity. However, LI values do not satisfy assumptions required of traditional linear regression because they are not normally distributed and are bounded between 0 and 1. To overcome these problems, we used the following beta regression model with a logit link function (Eq. (2)) to connect mean LI to the linear predictor: (2) $$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}\left(\mathrm{L}\mathrm{I}\right)=ln\left({\displaystyle \frac{\mathrm{L}\mathrm{I}}{1-\mathrm{L}\mathrm{I}}}\right)={\beta }_0+{\beta }_{\mathrm{i}}\mathrm{P}{\mathrm{C}}_{\mathrm{i}}+\dots +{\beta }_{\mathrm{k}}\mathrm{P}{\mathrm{C}}_{\mathrm{k}},\mathrm{i}=1,\dots ,\mathrm{k}$$where logit(LI) is the logit link function for the mean of LI and the linear predictor is defined by PC_i, …, PC_k as the scores of each principal component, β₀ as the intercept, β_i, …, β_k as coefficients corresponding to each principal component, and k as the number of principal components (Ferrari & Cribari-Neto, 2004). No 0-or 1-values were observed in laminarity, so we did not need to apply either a transformation to slightly shift boundary values or sophisticated zero-one inflated beta regression (Smithson & Verkuilen, 2006; Douma & Weedon, 2019). Instead, traditional beta regression was performed with the R package “gamlss” (Rigby & Stasinopoulos, 2005). To ease interpretation of fitted models on the scale of observed laminarity (0, 1), the linear predictor was transformed using the inverse logit function. Thus, the resulting expression (Eq. (3)) becomes a relationship between absolute changes in PC scores and mean LI: (3) $$\mathrm{L}\mathrm{I}={\displaystyle \frac{e^{\eta }}{1+e^{\eta }},}$$where the linear predictor η is β₀ + β_i PC_i + … + β_k PC_k (Douma & Weedon, 2019). Standardized coefficients of the original predictors (mass, bone length and Z_p) were calculated by multiplying the vector of β-coefficients by the matrix of eigenvectors from the PCA (Fekedulegn et al., 2002: Eq. 23).

Histological description

At mid-diaphysis, the cortical bone in the limbs of the sampled pigeons become increasingly compact with growth (Tables S2–S6). In very young individuals ranging from 0 to 2 weeks (0–40% of adult mass), bone walls are largely cancellous (porosity > 30%) with irregular vascular spaces. That cancellous structure is consistent with rapidly-growing juvenile bone as seen in other avian species (De Margerie et al., 2004). Older individuals show compact bone with vascular canals. For each bone, peak laminarity (i.e., proportion of circumferentially oriented canals) occurs in pre-fledge juveniles (Figs. 2 and 3) that range approximately from 2 to 4 weeks of age and 40–70% of adult mass. These juveniles also show large disparity in laminarity with much greater values occurring in the humerus, ulna and femur than in the radius and tibiotarsus. In post-fledge individuals (4–9 weeks of age and 70–100% of adult mass), the outer bone cortex is poorly vascularized with a subperiosteal region of avascular parallel-fibered bone in the oldest individuals of the sample (Figs. 3A–3J). The remaining deep cortex is highly vascularized, but canals are predominantly longitudinal with only slight differences in laminarity among elements (Figs. 3A–3J). Secondary osteons are generally uncommon in the sample. The notable exception is the tibiotarsus in which secondary osteons are abundant in the deep cortex towards the end of the sampled growth period.

Z_p scaling analysis

The polar section modulus (Z_p) of the humerus, radius, ulna, femur, and tibiotarsus increases with growth (Tables S2–S6). When scaled to log₁₀-transformed body mass, log₁₀-transformed Z_p shows significant positive allometry for all five sampled elements (Fig. 4) with 95% confidence intervals of the allometric exponents (a) each exceeding and excluding the isometric value of 1 (Table 1). The allometric exponent of the tibiotarsus is slightly (but not significantly) shallower than those of the other elements, in part reflecting the relatively large size of the tibiotarsus at hatch (Tables 1; Tables S2–S6).

Robust principal component analysis

Robust principal component analysis is generally consistent across the five limb elements (Table 2). PC1 captures at least 95% of the variance in the original predictors: 98% for the humerus, 95% for the radius, 97% for the ulna, 98% for the femur and 96% for the tibiotarsus. We ignored the residual variance (approximately 2–5%) that is absorbed by PC2 and PC3, thereby reducing data dimensionality from three components to one. Mass, Z_p (torsional rigidity), and bone length each have positive loadings on PC1. Z_p and bone length have strong effects on PC1, but dominance varies depending on the element. In the humerus and femur, Z_p is dominant or codominant with bone length, respectively, whereas in the remaining elements, length is dominant (Table 2). Taken together, the loadings are consistent with PC1 representing an ontogenetic axis. Small PC1 scores are associated with juvenile features (small mass with short bones that are relatively compliant to torsion), whereas large PC1 scores are associated with adult features (large mass with long bones that are relatively rigid to torsion).

Beta regression

The fitted beta regression models with a logit link have significant β-coefficients (Table 3). They predict that an absolute one-unit increase in PC1 leads to a relative change by a factor of exp(β₁) in the ratio of laminarity (LI) to non-laminarity (1-LI). Thus, there is a relative decrease in the ratio of laminarity to non-laminarity by: 22% in the humerus (exp(−0.249) − 1); 22% in the ulna (exp(−0.245) − 1); 18% in the radius (exp(−0.197) − 1); 22% in the femur (exp(−0.244) − 1); and 13% in the tibiotarsus (exp(−0.137) − 1) (Table 3). To better interpret these results as absolute changes in laminarity, the fitted beta regression model for each element was back-transformed with the inverse logit function and plotted. Each element shows a significant negative non-linear relationship between LI and PC1 (Fig. 5). Two groups are apparent. The first group, consisting of humerus, ulna, and femur, is characterized by relatively strong goodness-of-fit (pseudo-R² exceeds 0.70), relatively positive intercept (i.e., larger mean LI across ontogeny as calculated by Eq. (3)), and steep negative slope. In contrast, the second group, consisting of radius and tibiotarsus, is characterized by relatively weak goodness-of-fit (pseudo-R² < 0.55), relatively negative intercept (i.e., smaller mean LI across ontogeny), and shallow negative slope. Although laminarity generally decreases with ontogeny in the homing pigeon, laminarity in the radius and tibiotarsus may be influenced by additional unknown factors.

To express the PC1 coefficient in terms of the original predictors (mass, length and Z_p), the eigenvectors and coefficient of PC1 were multiplied by each other (Fekedulegn et al., 2002: Eq. 23). The result is a set of principal component estimators of the standardized effects of the original predictor on laminarity (Table 3). In all but two elements, laminarity is most influenced by the effect of bone length, which is at least 1.5 times greater than the effects of either Z_p or mass. Different patterns occur in the humerus and femur. In the former, the strongest effect is Z_p, and in the latter, the effects of bone length and Z_p are similar in strength. Contrary to expectations, none of the standardized coefficients are positive indicating that, at least in the homing pigeon, growth in mass, bone length, and torsional rigidity have inverse effects on laminarity.

Positive allometric growth of torsional rigidity at midshaft

This study demonstrates that some structural changes in midshaft cortical bone are predictable responses to delayed locomotor development. Juvenile pigeons receive extended parental care within protected nests and have limited mobility for over half of their postnatal growth period (Levi, 1962; Janiga & Kocian, 1985; Johnston & Janiga, 1995; Vriends & Erskine, 2005; Liang et al., 2018). Under those conditions, selection for relatively robust juvenile limbs is likely relaxed and reflected in “underbuilt” midshaft cortices. As juvenile pigeons approach adult size and become fully mobile, torsional rigidity increases rapidly in forelimb and hindlimb bones as indicated by significant positive ontogenetic allometry of midshaft Z_p (Fig. 4).

Other volant birds also show positive ontogenetic scaling of geometric properties in forelimb bones. Like pigeons, larids and mallards hatch with altricial wings and are unable to fly until juveniles are nearly adult size. Late development of wings involves accelerated growth in pectoral muscle mass, wing surface area, and breaking strength of bones. The latter is attributed in part to disproportionately large increases to forelimb bone midshaft width and second moment of area (Carrier & Leon, 1990; Bennett, 2008; Dial & Carrier, 2012). Late investment in wings may allow resource allocation to develop other metabolically expensive organ systems first (e.g., digestive, nervous, cardiovascular, and integumentary) (Carrier & Leon, 1990; Johnston & Janiga, 1995; Dial & Carrier, 2012; Altimiras et al., 2017). A tentative conclusion from these results is that altriciality and positive allometric growth are coupled at least in avian wing development. However, we are aware that semi-intensively farmed Japanese quail from China show positive allometric growth of midshaft width in forelimb bones (Ren, Wang & Zhang, 2016), which is odd given the early locomotor capabilities of the wings. Positive allometry of midshaft width was also found in the femur and tibiotarsus, the latter being inconsistent with results from the closely-related chicken (Biewener & Bertram, 1994). Because the analysis on those Japanese quail used Type II regression, which has tendency to upward bias allometric slopes (Kilmer & Rodríguez, 2017), we suspect that negative allometry or isometry was mistaken for positive allometry. In any event, additional sampling is needed to test whether avian taxa with precocial wings capable of flight shortly after hatching (e.g., the Megapodidae (Dial & Jackson, 2011)) show negative allometric or isometric growth of forelimb bones, presumably to keep breaking loads relatively constant as they function in flight while still growing.

Unlike the pigeon, the California gull and mallard have hindlimbs that are functional for walking or swimming shortly after hatching. Their juveniles have relatively robust hindlimb bones and generally experience negative allometric or isometric growth in midshaft width and second moment of area, presumably to maintain locomotor performance and bone strength comparable to those of adults (Carrier & Leon, 1990; Dial & Carrier, 2012). The glaucous-winged gull also has functional hindlimbs at hatching. However, results from an ontogenetic study are not presented in terms of allometric scaling to body mass (Hayward et al., 2009). Instead, they show that sigmoidal growth in midshaft width in hindlimb bones peaks relatively early and is completed before forelimb bones. That pattern is consistent with negative allometry and supports the hypothesis that growth in limb bone geometry reflects locomotor needs during development.

Although distinct allometric growth trajectories of hindlimb bones may, to an extent, indicate locomotor ability in juveniles, there are at least two examples that caution against oversimplification. Like previously discussed larids (Carrier & Leon, 1990; Hayward et al., 2009), the Black noddy shows negative allometric growth in second moment of area of the tibiotarsus (Bennett, 2008). However, the Black noddy is tree- or cliff-nesting, which is a behavior that provides refuge to juveniles but limits their locomotion during the growth period. Nevertheless, juveniles still grow tibiotarsii with relatively thick cortices that compensate for low mineralization to maintain bending strength and stiffness comparable to adults (Bennett, 2008). Allocation of resources to develop relatively robust locomotor-capable hindlimb bones in nest-bound juveniles is counterintuitive, and negative allometric growth in hindlimb bones may simply be a relic of ancestral selection for locomotor ability in juveniles shared among larids (Bennett, 2008).

Emus are capable of walking within two days after hatching (Sales, 2007). Yet femora and tibiotarsii of chicks and early juveniles tend to be “underbuilt” relative to adults as indicated by positive allometry in polar moment of area (Main & Biewener, 2007). This ontogenetic scaling relationship differs from the negative allometric growth reported in some precocial mammals (e.g., jackrabbit and goat), which increases bone safety factors in juveniles when performing locomotion similar to adults (Carrier, 1983; Main & Biewener, 2004, 2007). The reason for the difference in hindlimb bone growth allometry among precocial species is not clear, but captivity is not likely a strong factor. Although emus and goats were studied under similar captive conditions with free access to exercise, they show divergent allometric patterns in hindlimb bone growth of polar moment of area (Main & Biewener, 2004, 2007). Differences in behavioral ecology was noted as a potential explanation; when compared to young goats, juvenile emu spend more time foraging at slower speeds instead of evading threats at faster speeds. Therefore, the “underbuilt” cortices of hindlimb bones in juvenile emu may reflect relaxed selection for high-performance locomotion (Main & Biewener, 2006). In a similar way, juvenile pigeons may experience relaxed selection for torsionally rigid limb bones given that they delay development of aerial and terrestrial locomotion until nearly full-grown.

Limb bone laminarity decreases with maturity

Contrary to previous work on domestic turkey (Skedros & Hunt, 2004) and emu (Kuehn et al., 2019), our results in the homing pigeon do not support the hypothesis that laminarity increases as a developmental response to locomotor-induced torsion. Unlike midshaft cortical geometry, which appears to respond to locomotor maturation by growing progressively rigid to torsion, laminarity decreases such that locomotor-capable adults have low values (2.4–15.8%; Figs. 3 and 5). Adults have low laminarity bone in part because the growth of the medullary cavity resorbs all record of earlier juvenile bone with higher laminarity. The remaining cortex with fewer circumferential canals (i.e., lower laminarity) reflects gradual reduction in periosteal deposition with skeletal maturation (Videos S1–S5). According to the laminarity hypothesis, low laminarity is a feature of limb bones adapted to locomotion involving reduced torsional loading. For example, low laminarity in the forelimb bones of birds with long and narrow wings may reflect habitual loading in bending rather than torsion (De Margerie et al., 2005; Simons & O’Connor, 2012). The juvenile and young adult homing pigeons in our sample were raised in enclosed lofts with minimal crowding and free access to exercise, but without in vivo bone strain data, we recognize that biomechanical inferences based on our sample of homing pigeons remain speculation. Nevertheless, there is consilience between cross-sectional geometry from the current study and in vivo off-axis principal strains from other studies of adult birds, including feral pigeons (Lanyon & Rubin, 1984; Biewener, Swartz & Bertram, 1986; Biewener & Dial, 1995; Carrano & Biewener, 1999; Main & Biewener, 2007), that suggests substantial locomotor-induced torsion at least in the humerus, ulna, femur and tibiotarsus. A circumspect conclusion drawn from limited data is that bones adapted to resist torsional loading do not necessarily increase laminarity during development.

Recent work on adult feral pigeons reveals “high” laminarity at least in the humerus (Skedros & Doutré, 2019). The contrasting results suggest that laminarity in homing pigeons may not be as representative of the wild-type as assumed by the present study. However, the assessment of “high” laminarity in feral pigeons is based on casual inspection rather than a canal-by-canal count. Therefore, we cannot exclude the possibility that it is an overestimate. Indeed, when canal-by-canal counts are performed in other avian taxa either with histological (Lee & Simons, 2015) or computed tomographic methods (Pratt & Cooper, 2017), laminarity that is previously reported to be “high” changes to “low” (Pratt et al., 2018). Similarly, we expect a reassessment to find low laminarity in adult feral pigeons.

The focus on circumferential canals in the laminarity hypothesis needs better biomechanical justification. If vascular canals align circumferentially to reflect torsional loading on a given bone, then the canals should be roughly concentric with an approximate angle of 45° from the longitudinal axis. This angle corresponds to the orientation of peak principal strain during pure torsional loading (Craig, 2000). Empirical data from in vivo bone strain experiments on birds during locomotion demonstrate that the orientation of peak principal strain varies with limb bone and species but is still reasonably close to 45° (Biewener, Swartz & Bertram, 1986; Biewener & Dial, 1995; Carrano & Biewener, 1999; Main & Biewener, 2007). In contrast, circumferential canals are defined as “in-plane” features (De Margerie, 2002) that are aligned closer to the transverse plane of section. Because vascular canals are approximately cylindrical (Fig. 1), orientation can be estimated from a histological section by applying the in-plane aspect ratio of the canal (i.e., greater than 3: De Margerie, 2002) to a trigonometric equation (De Boef & Larsson, 2007)—angle in degrees = 180/π cos⁻¹(aspect ratio⁻¹). Using that equation, we calculate that circumferential canals are angled in excess of 45° and range between 70.5° and 90° from the longitudinal axis. This discrepancy suggests that circumferential canals are not as adapted to torsion as originally thought. Instead, canals once classified as “longitudinal” may be better torsion-resisting features. For example, we estimate that “longitudinal canals” with the aspect ratio of 1.41 would be aligned with peak principal strains in bone under torsion. The revised biomechanical interpretation may explain why “longitudinal canals” are so abundant not only in homing pigeons but generally in amniotes. Further testing is needed but requires shifts from: (1) assumptions to experiments; (2) two-dimensional to three-dimensional imaging; and (3) categorical to continuous data.

The developmental approach used by the current study may inform how the nanostructural organization of collagen fibers also contribute to the torsional rigidity of bone. In adult birds, collagen fibers with oblique-to-transverse orientation (i.e., spiraling 45°–90° to the longitudinal axis) are especially abundant throughout the cortex of bones shaped to resist torsion (De Margerie et al., 2005). Similar collagen fiber orientation evolved independently in adult birds and at least one species of fruit bat (Skedros & Doutré, 2019) suggesting that it may be a fundamental adaptation of vertebrate flapping flight. If so, we expect collagen fiber obliquity and torsional rigidity of wing bones to increase with locomotor maturity. Preliminary evidence suggests that the predicted trend occurs at least in the ulna of growing turkey (Skedros et al., 2003; Skedros & Hunt, 2004). Future investigations should apply the developmental approach across a broader phylogenetic sample to test whether loading has similar effects on collagen fibers and vascular canals.