Size And Locomotor Ecology Have Differing Effects on the External and Internal Morphologies of Squirrel (Rodentia: Sciuridae) Limb Bones

Synopsis Mammals exhibit a diverse range of limb morphologies that are associated with different locomotor ecologies and structural mechanics. Much remains to be investigated, however, about the combined effects of locomotor modes and scaling on the external shape and structural properties of limb bones. Here, we used squirrels (Sciuridae) as a model clade to examine the effects of locomotor mode and scaling on the external shape and structure of the two major limb bones, the humerus and femur. We quantified humeral and femoral morphologies using 3D geometric morphometrics and bone structure analyses on a sample of 76 squirrel species across their four major ecotypes. We then used phylogenetic generalized linear models to test how locomotor ecology, size, and their interaction influenced morphological traits. We found that size and locomotor mode exhibit different relationships with the external shape and structure of the limb bones, and that these relationships differ between the humerus and femur. External shapes of the humerus and, to a lesser extent, the femur are best explained by locomotor ecology rather than by size, whereas structures of both bones are best explained by interactions between locomotor ecology and scaling. Interestingly, the statistical relationships between limb morphologies and ecotype were lost when accounting for phylogenetic relationships among species under Brownian motion. That assuming Brownian motion confounded these relationships is not surprising considering squirrel ecotypes are phylogenetically clustered; our results suggest that humeral and femoral variation partitioned early between clades and their ecomorphologies were maintained to the present. Overall, our results show how mechanical constraints, locomotor ecology, and evolutionary history may enact different pressures on the shape and structure of limb bones in mammals.

Synopsis Mammals exhibit a diverse range of limb morphologies that are associated with different locomotor ecologies and structural mechanics. Much remains to be investigated, however, about the combined effects of locomotor modes and scaling on the external shape and structural properties of limb bones. Here, we used squirrels (Sciuridae) as a model clade to examine the effects of locomotor mode and scaling on the external shape and structure of the two major limb bones, the humerus and femur. We quantified humeral and femoral morphologies using 3D geometric morphometrics and bone structure analyses on a sample of 76 squirrel species across their four major ecotypes. We then used phylogenetic generalized linear models to test how locomotor ecology, size, and their interaction influenced morphological traits. We found that size and locomotor mode exhibit different relationships with the external shape and structure of the limb bones, and that these relationships differ between the humerus and femur. External shapes of the humerus and, to a lesser extent, the femur are best explained by locomotor ecology rather than by size, whereas structures of both bones are best explained by interactions between locomotor ecology and scaling. Interestingly, the statistical relationships between limb morphologies and ecotype were lost when accounting for phylogenetic relationships among species under Brownian motion. That assuming Brownian motion confounded these relationships is not surprising considering squirrel ecotypes are phylogenetically clustered; our results suggest that humeral and femoral variation partitioned early between clades and their ecomorphologies were maintained to the present. Overall, our results show how mechanical constraints, locomotor ecology, and evolutionary history may enact different pressures on the shape and structure of limb bones in mammals.

Introduction
Mammals exhibit a variety of locomotor modes to transverse across a wide range of habitats ( Hildebrand 1985 ). Adaptations to these various locomotor modes are repeatedly observed in mammalian limb bones ( Polly 2007 ). For example, scansorial mammals tend to exhibit more elongate, gracile limb bones ( Alexander et al. 1979 ;Burr et al. 1989 ;Kimura 1991 ;Polly 2007 ), whereas fossorial mammals exhibit relatively shorter and more robust limb bones in response to stress enacted on the bone when digging ( Peterka 1936 ;Gasc et al. 1985 ;Straehl et al. 2013 ;Montoya-Sanhueza and Chinsamy 2017 ). Furthermore, different locomotor modes used in varying environments necessitate long bones being able to resist possible deformation to locomotion-specific forces. For example, species that fly, glide, or leap exhibit relatively long limbs with circular-shaped cross sections in the diaphyses to resist high torsional and multidimensional bending forces ( Burr et al. 1989 ;Swartz et al. 1992 ;Patel et al. 2013 ;Hunt et al. 2021 ). In contrast, highly fossorial mammals exhibit large amounts of compact cancellous bone in the forelimb with elliptical-shaped cross sections to withstand uniaxial bending loads associated with digging ( Amson et al. , 2022. Limb bone morphology is also influenced by body size ( Alexander et al. 1979 ;Biewener 1983 ;Christiansen 1999 ). As body size increases, more mechanical support is needed to compensate for an increase in the mechanical forces exerted during locomotion ( McMahon 1973 ;Biewener 1990 ;Christiansen 1999 ), which may result in adaptations in the shape and structural properties of the limb bones. Increasing body size is associated with increasing bone robustness in multiple lineages of tetrapods ( Alexander et al. 1979 ;Demes and Jungers 1993 ;Christiansen 1999 ;Doube et al. 2009 ;Ryan and Shaw 2013 ;Mielke et al. 2018 ). In quadrupedal terrestrial tetrapods, body size rather than locomotor mode has a stronger influence on conserved limb bone structural traits such as the minor diaphyseal circumference of weight-bearing bones and the ratio between the stylopodial circumference and body mass ( Campione and Evans 2012 ). Body size may also impose limits on bone shape; some cursorial carnivorans without the ability to supinate their forelimbs appear to have smaller maximum body sizes, while carnivorans with more generalized bone shapes exist at much larger body sizes ( McMahon 1973 ;Polly 2007 ;Fabre et al. 2015 ).
Unsurprisingly, a plethora of work has examined the influence of locomotor behavior or body size on the overall shape of limb bones (e.g., Fabre et al. 2013 ;Martín-Serra et al. 2014a, 2014bHedrick et al. 2020 ;Etienne et al. 2021 ), and bone structural characteristics-such as bone compactness, diaphysis elongation, or cross-sectional shape-that are most directly impacted by different mecha nical loads (e.g., Currey and Alexander 1985 ;Schaffler et al. 1985 ;Burr et al. 1989 ;Patel et al. 2013 ;Kilbourne and Hutchinson 2019 ;Scheidt et al. 2019 ;Wölfer et al. 2019 ;Amson et al. 2022 ). In this study, we further our understanding of limb bone adaptations by testing how the interaction of locomotor ecology and scaling influences the shape and bone structure of the humerus and femur, using squirrels (Sciuridae) as a model clade. The squirrel family consists of approximately 280 species that inhabit a variety of microhabitats, from far above ground in tree canopies to deep underground in burrow systems. Squirrel species can be categorized into four ecotypes with distinct locomotor ecologies and behaviors: ground squirrels that dig, tree squirrels that climb, gliding squirrels that glide between trees, and more versatile chipmunks that both dig and climb. Previous work has examined relationships between limb lengths and ecotypes, finding, for example, that gliding squirrels exhibit relatively longer forelimbs than all other ecotypes ( Peterka 1936 ;Bryant 1945 ;Thorington and Heaney 1981 ;Grossnickle et al. 2020 ;Linden et al. 2023 ), whereas ground squirrels exhibit relatively shorter forelimbs with increasing body elongation ( Linden et al. 2023 ). Furthermore, across Sciuromorpha (squirrel-like rodents), cross-sectional characteristics (e.g., cross-sectional area and second moment of area) of the femur scale with positive allometry with respect to body mass ( Scheidt et al. 2019 ). Nevertheless, whether the scaling patterns of external shape and bone structure of the humerus and femur differ among ecotypes remains to be tested within this clade. We predicted that humeral and femoral morphology is best explained by an interaction of ecotype and size, where both limb bones would be more robust, compact, and exhibit a more oval-shaped cross section with increasing bone size in grounddwelling (i.e., ground, tree, chipmunks) squirrels to resist the relatively greater mechanical loads that occur at larger body sizes ( McMahon 1973 ;Biewener 1990 ;Christiansen 1999 ). Conversely, gliding squirrels would exhibit opposite trends. Within ecotypes, we predicted that ground squirrels would have more robust, compact limb bones with oval-shaped cross-sections to reinforce the bone in the cranial-caudal plane during scratch digging ( Hildebrand 1985( Hildebrand , 1995Lieberman et al. 2004 ;Lagaria and Youlatos 2006 ). We predicted that chipmunks and tree squirrels would exhibit more generalized limb bone shapes in comparison to other ecotypes, as they employ running, climbing, and digging behaviors. Finally, we predicted that gliding squirrels would exhibit more elongate, gracile limb bones with circularshaped cross sections to increase patagia surface area for gliding.

Mor pholog ical data
We acquired three-dimensional scans of humeri and femora from 76 squirrel species (one specimen per species) through computed tomography (CT) scanning with Skyscan1172 μCT, Skyscan1173 μCT, GE Phoenix Nanotom M, and NSI X5000 scanning systems. Scans were performed with a resolution of 26.05 μm. Scan data were reconstructed using NRecon and exported and segmented in 3D Slicer ( Kikinis et al. 2014 ). All specimens were sourced from the collections of 11 museums (Table S1). We used female, male, and sex-unknown individuals for our measurements to achieve the largest number of species. We determined that each specimen was fully mature by verifying that the cranial exoccipital-basioccipital and basisphenoidbasioccipital sutures were fused and that all vertebrae and limb bones were ossified.

Ecotype classification
Following Linden et al. (2023) , we categorized species into four ecotypes-chipmunk (n = 14), gliding (n = 8), ground (n = 25), or tree (n = 29)-based on locomotion, evolutionary grouping, and reproductive behavior ( Fig. 1 ). Species that display fossorial locomotion and reproduce in underground burrows were categorized as ground squirrels, species that display both arboreal and scansorial locomotion and reproduce in nests in trees as were categorized tree squirrels, and species that have patagia were categorized as gliding squirrels. Our fourth ecotype group was chipmunks (genus Tamias ), which display the broadest range of locomotor and nesting behaviors; species are considered terrestrial, semifossorial, or semi-arboreal depending on the source consulted, but none are considered fully fossorial or arboreal.

External shape analyses
All phylogenetic analyses were conducted using a pruned version of the Upham et al. (2019) phylogeny.
All statistical analyses were performed in R (2022). We quantified humeral and femoral size and shape using three-dimensional geometric morphometrics ( Rohlf and Slice 1990 ;Zelditch et al. 2012 ). We first transformed all specimens to correspond to the left side (an arbitrary choice) of the body using the Transforms module in 3D Slicer ( Fedorov et al. 2012 ). We generated 194 pseudolandmarks on a humerus template and 152 pseudolandmarks on a femur template using the PseudoLMGenerator module of the SlicerMorph extension  in 3D Slicer ( Fedorov et al. 2012 ). The most average squirrel in our dataset (determined visually on the phylomorphospace as the Western gray squirrel, Sciurus griseus ) was used as our template. Following Diamond et al. (2023) , we transferred the pseudolandmarks from the template model to each individual humerus and femur 3D model using the ALPACA module  in SlicerMorph ( Fig. 1 ). We conducted a General Procrustes Analysis (GPA) using the gpagen function in the R package geomorph ( Baken et al. 2021 ;Adams et al. 2022 ). We then verified that there were no significant shape differences between limb bones originating from the left or right side with Procrustes ANOVAs (humerus: F = 1.10, Z = 0.49, p = 0.319; femur: F = 1.38, Z = 1.03, p = 0.156).
We used humeral and femoral centroid size as our metrics of humeral and femoral size, respectively. We tested whether humeral and femoral size differed among ecotypes using phylogenetic analyses of variance (ANOVAs), which jointly estimate Pagel's λ to account for phylogenetic covariance present in the model residuals, using the R package phylolm v2.6.2 ( Tung Ho and Ané 2014 ). We assessed differences in limb sizes among ecotypes by first generating 1000 bootstrap replications of the model coefficients of limb size. We then computed the observed difference between the mean limb sizes of each ecotype pair and created its 95% confidence interval using the bootstrap replications. Confidence intervals that encompassed zero indicated that ecotype pairs were not significantly different from each other.
We visualized the phylomorphospace of humeral and femoral shape by performing principal component analyses (PCA) using the gm.prcomp function in geomorph ( Baken et al. 2021 ;Adams et al. 2022 ). We also visualized differences between PC extremes by creating average shape models using the Morpho package in R ( Schlager 2017 ) and used a custom Python script ( Diamond et al. 2023 ) to extract the mesh distances from the areas with pseudolandmark points.
We then tested how size, ecotype, and their interaction influenced humeral and femoral shape, respectively, using a series of phylogenetic Fig. 1 Over vie w of the species and morphometrics used. The pruned phylogeny is overlaid with an ancestral state reconstruction of squirrel ecotype. A total of 194 pseudolandmarks on the humeri and 152 pseudolandmarks on the femora were generated and applied to our sample using the PseudoLMGenerator  ) and ALPACA modules  in SlicerMorph. Bone structure variables were calculated using the SegmentGeometry module ( Huie et al. 2022 ) in Slicer. Cg = global compactness; DE = diaphysis elongation; CSS = crosssectional shape.
penalized-likelihood multivariate generalized least squares (mvPGLS) models. These mvPGLS models use a penalized-likelihood to fit linear models to highdimensional data sets in which the number of variables is much larger than the number of observations ( Clavel et al. 2019 ). We fit our first model to test how each bone shape scaled with bone size (i.e., mvPGLS size model: bone shape ∼ ln bone size). With the second model, we tested the relationship between bone shape and ecotype (i.e., mvPGLS ecotype model: bone shape ∼ ecotype). With the third model, we tested if size and ecotype influenced bone shape (i.e., mvPGLS size + ecotype model: bone shape ∼ ln bone size + ecotype). With the fourth model, we tested how the interaction of size and ecotype influenced bo ne shape (i.e., mvPGLS size*ecotype model: bone shape ∼ ln bone size*ecotype). The last model we fit was a null model (mvPGLS null model) in which size and ecotype have no influence on bone shape (i.e., bone shape ∼ 1). We assessed the importance of size, ecotype, their combined effect, and their interaction on limb shape variation by calculating the relative support for each of the five models through computation of Extended Information Criterion (EIC) weights (EICw) with 1000 bootstrap replications. We also assessed the scaling pattern of bone shape within each ecotype (i.e., bone shape ecotype ∼ ln bone size ecotype ). All mvPGLS models jointly estimated Pagel's λ with the model parameters using the mvgls function ( Clavel et al. 2019 ;Clavel and Morlon 2020 ) in the R package mvMORPH v.1.1.6 ( Clavel et al. 2015 ). We performed the multivariate test using Pillai's trace with 1000 permutations using the manova.gls function. We found that ecotype was an important predictor of bone shape (see Results section); however, we were unable to perform post hoc pairwise testing among ecotypes because the mvMORPH function pairwise.glh is only suited for binary data thus far. Therefore, we performed post hoc pairwise permutation tests using the function pairwise in RRPP v1.3.1 ( Collyer et al. 2015 ;Adams and Collyer 2018 ). These pairwise tests required us to repeat our analyses among ecotypes using Procrustes ANCOVAs in the R package geomorph v4.0.4 ( Baken et al. 2021 ;Adams et al. 2022 ). Because our initial mvPGLS results indicated there was phylogenetic signal in the residuals of the model (Pagel's λ = 0.18), we performed both (a) Procrustes ANCOVAs that assume no phylogenetic structure and (b) phylogenetic Procrustes ANCOVAs that assume Brownian motion to assess how the strength of phylogenetic signal in the residuals influenced our results. Procrustes AN-COVAs were performed with 1000 random residual permutation procedure (RRPP) using the functions procD.lm and procD.pgls in geomorph v4.0.4. Lastly, we estimated how well ecotype can predict limb shape using canonical variate analyses (CVA) with a jackknife cross-validation procedure in the R package Morpho ( Schlager 2017 ).

Bone structure analyses
We quantified three bone structural traits-global compactness (Cg), diaphysis elongation (DE), and cross-sectional shape (CSS)-that serve as proxies for biomechanical function (e.g., Patel et al. 2013 ;Kilbourne and Hutchinson 2019 ;Hedrick et al. 2020 ;Amson and Bibi 2021 ). Cg can reflect resistance towards axial compression where more compact limb bones (i.e., higher Cg) tend to resist more compressive loads ( Berman et al. 2015 ). We used Cg instead of cross-sectional area because we aimed to quantify how much material is invested in the diaphysis of the limb bones while accounting for cross-sectional size. DE can be informative of the mechanical advantage of the limb movement, where relatively longer limb bones (i.e., higher DE values) reduce the outlever whereas relatively shorter limb bones (i.e., lower DE values) increase the outlever ( Smith and Savage 1956 ). CSS can reflect torsion and uni-or multidirectional bending of the humeral and femoral diaphyses, where bones with more elliptical cross sections (i.e., higher CSS values) are hypothesized to be adapted for bending loads in one direction, whereas bones with more circular cross sections (i.e., lower CSS values) are hypothesized to be adapted for the multidirectional bending loads that are associated with using diverse locomotor behaviors or resisting torsion during aerial locomotion ( Swartz et al. 1992 ;Patel et al. 2013 ;Amson et al. 2022 ). We quantified these three traits by creating slice-by-slice proximodistal profiles of each scanned specimen in the Slicer module SegmentGeometry ( Huie et al. 2022 ). Each profile consisted of length (in millimeters), crosssectional area (in square millimeters), second moment of area around the major principal axis, second moment of area around the minor principal axis, and compactness (ratio between sample cross-sectional area and total cross-sectional area). Outputs from SegmentGeometry were used to calculate three bone structure traits: Cg, defined as the mean compactness of all included slices; DE, defined as the ratio between the functional length of the bone and the square root of the cross-sectional area at the midshaft slice; and CSS, defined as the ratio of the second moment of area around the major principal axis to the second moment of area around the minor principal axis at the midshaft slice ( Fig. 1 ). Each bone structure trait was computed across the middle 40% of the bones' functional length, which captures the greatest distance between the proximal and distal articular surfaces of the bones and ensures only the diaphysis is measured.
We tested how size, ecotype, and their interaction influenced each bone structure trait (i.e., Cg, DE, CSS) using a series of phylogenetic generalized least squares (PGLS) models. Similar as above, we assessed the importance of size, ecotype, their combined effect, and their interaction on bone structure traits by fitting five models (i.e., PGLS size , PGLS ecotype , PGLS size + ecotype , PGLS size*ecotype , PGLS null ) and calculating the relative support for each model through computation of Akaike information criterion weights (AICw). All models with AIC below 2 were considered to be supported by the data ( Burnham and Anderson 2004 ). Regression coefficients for all models were estimated simultaneously with phylogenetic signal as Pagel's λ in the residual error using the R package phylolm v2.6.2 ( Tung Ho and Ané 2014 ). We also generated 1000 bootstrap replications of the scaling slopes and means of each variable in each PGLS model. In instances where the best fitting model incorporated ecotype, we assessed if bone structure variables differed among ecotypes by computing the observed difference between the mean values of each ecotype pair and creating its 95% confidence interval using the bootstrap replications. Confidence intervals that encompassed zero indicated that ecotype pairs were not significantly different from each other. To maintain consistency with the external shape analyses, we also assessed if bone structure variables differed among ecotypes using (a) ANCOVAs that assume no phylogenetic structure and (b) phylogenetic AN-COVAs that assume full Brownian motion. All ANCO-VAs were performed with 1000 RRPP iterations using the function lm.rrpp in the R package RRPP . In instances where the best fitting model incorporated size, we determined if bone structure scaling was significantly different from isometry across all squirrels. To do this, we used the bootstrap replications to create 95% confidence intervals for the mean scaling slopes across all ecotypes as well as for each individual ecotype. Confidence intervals greater than zero were considered to represent positive allometry, while confidence intervals less than zero were considered to represent negative allometry. Confidence intervals that included zero were considered isometric due to the dimensionless quality of the bone structure variables. In instances where the best fitting model incorporated the interaction between ecotype and size, we assessed if scaling patterns of bone structure variables differed a mong ecotypes by computing the observed difference between the mean slopes of each ecotype pair and created its 95% confidence interval using the bootstrap replications. Confidence intervals that encompassed zero indicated that ecotype pairs were not significantly different from each other. We also assessed if bone structure scaling was significantly different from isometry within each ecotype by using the bootstrap replications to create 95% confidence intervals for the mean scaling slopes for each individual ecotype. Confidence intervals greater than zero were considered to represent positive allometry, while confidence intervals less than zero were considered to represent negative allometry. Confidence intervals that included zero were considered isometric due to the dimensionless quality of the bone structure variables.

Comparisons between bone shape and structure, and between humeral and femoral traits
We tested if there were relationships between the external shape and each of the bone structure variables of the humerus and femur across squirrels and within each ecotype using two-block partial least squares (PLS) analyses in the R package geomorph version 4.0.4 ( Baken et al. 2021 ;Adams et al. 2022 ). We also tested if there were relationships between humeral traits and femoral traits across all squirrels and within each ecotype. We used PLS analyses in geomorph to test relationships between humeral shape and femoral shape, and phylogenetic paired t-tests in the phytools function phyl.pairedttest to test if each bone structure variable differed between the humerus and femur.
Humeral external shape PC1 accounted for 36.8% of the humeral shape variation and largely separates the morphospaces among ground squirrels, chipmunks and tree squirrels, and gliding squirrels ( Fig. 3 ). PC1 primarily describes the lengthening of the humerus. Positive PC1 scores are associated with more robust humeri with more pronounced supinator crests, deltoid tuberosities, and medial epicondyles, whereas negative PC1 scores are Fig. 3 Phylomorphospaces of (A) humeral shape and (B) femoral shape with PC1-2 extremes for visualization. The humeral and femoral models reflect the computationally determined PC1 and 2 extremes, and their colors reflect the ecotype associated with extreme ends of the first two PCs. Blue circles represent g liding squir rels, yellow circles represent ground squirrels, navy circles represent chipmunks, green circles represent tree squirrels, and open circles represent tree nodes within the phylomorphospace. associated with more elongate humeri with less developed supinator crests, deltoid tuberosities, and medial epicondyles. PC2 accounted for 12.5% of the humeral shape variation and is associated with shape variation at the ends of the humerus ( Fig. 3 ). Positive PC2 scores are associated with shape variation at the medial epicondyle and negative PC2 scores are associated with changes at the condyle. Much of this variation across PC2 explains within-ecotype differences.
The mvPGLS ecotype model was the best supported model for humeral shape (EICw = 1.00; Table 1 ) and showed a statistically significant relationship between humeral shape and ecotype (Pagel's λ = 0.19, Pillai's trace = 2.83, P < 0.001). Because we were unable to perform post-hoc pairwise tests among ecotypes using the PL-MANOVA model (see Methods section), we performed pairwise permutation tests using both phyloge-netic Procrustes ANCOVA that assumes Brownian motion and Procrustes ANCOVA that does not incorporate phylogenetic structure. Unlike the mvPGLS ecotype model, the phylogenetic Procrustes ANCOVA model indicated that there was no significant relationship between humeral shape and ecotype (R 2 = 0.04, Z = −0.137, P = 0.550). In contrast, the nonphylogenetic Procrustes ANCOVA model indicated that there was a significant relationship between humeral shape and ecotype (R 2 = 0.37, Z = 6.551, P < 0.001), and posthoc pairwise tests revealed that there were significant differences in humeral shape among all ecotypes (all P < 0.002) except between chipmunks and tree squirrels (P = 0.075). This is consistent with the CVA with jackknife cross-validation, which reclassified humeral shapes in their correct ecotype with 93.4% accuracy overall and accurately reclassified 100% of gliding and Table 1 Comparisons of the best-fitting phylogenetic generalized least squares (PGLS) models in humeral and femoral shape and bone structure. Extended Inf or mation Criterion (EIC) and Akaike Inf or mation Criterion (AICc) were used to assess model fits on bone shape and bone structure, respecti vel y ground squirrels, 85.7% of chipmunks, and 89.7% of tree squirrels. Within ecotypes, gliding squirrels (Pagel's λ = 0.00, Pillai's trace = 1.00, P = 0.016) and ground squirrels (Pagel's λ = 0.21, Pillai's trace = 0.98, P = 0.009) exhibited significant humeral shape allometry whereas chipmunks (Pagel's λ = 0.00, Pillai's trace = 0.98, P = 0.382) and tree squirrels (Pagel's λ = 0.26, Pillai's trace = 0.94, P = 0.742) did not ( Fig. 4 ).
Femoral external shape PC1 accounted for 29.3% of the femoral shape variation and separated the morphospace of gliding squirrels from a morphospace largely containing ground squirrels, chipmunks, and tree squirrels ( Fig. 3 ). Positive PC1 scores are associated with more robust femora with more pronounced femoral heads, patellar grooves, greater trochanters, third trochanters, and lesser trochanters whereas negative PC1 scores are associated with more elongate, gracile femora ( Fig. 3 ). PC2 accounted for 12.8% of the femoral variation and loosely separated tree squirrels from ground squirrels and chipmunks. Positive PC2 scores are associated with shape variation around the greater trochanter and one side of the lesser trochanter, whereas negative PC2 scores are associated with shape variation around the femoral head, greater trochanter, and the opposite side of the lesser trochanter (i.e., one extreme reflects shape Fig. 4 Allometric trends between external shape and size of the (A) humerus and (B) femur. Limb bone images show the difference in shape from the smallest to biggest species in each ecotype (chipmunk: Tamias alpinus to Tamias townsendii ; gliding: Glaucomys volans to Petaurista petaurista ; ground: Callospermophilus madrensis to Marmota caligata ; tree: Tamiops maritimus to Ratufa macroura ).
The mvPGLS null model was the best supported model for femoral shape (EICw = 1.00; Table 1 ), suggesting that neither size nor ecotype are important predictors of femoral shape. Nevertheless, the CVA accurately reclassified femoral shapes in their correct ecotype with 90.8% accuracy overall and accurately reclassified 100% of gliding squirrels, 96.6% of tree squirrels, 92.9% of chipmunks, and 80.0% of ground squirrels. Consistent with the CVA, results from the Procrustes ANCOVA model that does not incorporate phylogenetic structure indicated that there was a significant relationship between femoral shape and ecotype (R 2 = 0.31, Z = 6.88, P < 0.001), and post-hoc pairwise tests revealed that there were significant differences in femoral shape among all ecotypes (all P < 0.004). However, incorporating phylogenetic structure resulted in no significant relationship between femoral shape and ecotype (mvPGLS ecotype model with Pagel's λ = 0.32: Pillai's trace = 2.57, P = 0.194; phylogenetic Procrustes ANCOVA under Brownian motion: R 2 = 0.04, Z = −0.446, P = 0.676).
Ecotype is a more important predictor of femoral CSS than size, as PGLS ecotype was the best supported model (R 2 = 0.40, λ = 0.00, AICw = 0.64; Table 1 Table S2). Mean differences in f emoral CSS a mong ecotypes were corroborated by the ANOVA (R 2 = 0.40, P = 0.001; all pairwise comparisons P < 0.006 except for the chipmunk-ground and gliding-tree comparisons) but not the phylogenetic ANOVA (R 2 = 0.01, P = 0.834).

Relationships between external shape and bone structure
We found that humeral shape exhibited significant relationships with DE (r-PLS = 0.88, Z = 4.547, P < 0.001) and CSS (r-PLS = 0.62, Z = 3.257, P < 0.001) but not Cg (r-PLS = 0.39, Z = 0.956, P = 0.186; Table S3). Within ecotypes, only tree squirrel humeral shape exhibited significant relationships with all three bone structure traits, and ground squirrel humeral shape exhibited a significant relationship with humeral DE and humeral CSS (Table S3).

Discussion
Size, locomotor ecology, and phylogenetic history are often identified as the leading factors that underlie limb bone variation in mammals (e.g., Fabre et al. 2013 ;Kilbourne and Hoffman 2013 ;Martín-Serra et al. 2014b, 2014aMielke et al. 2018 ;Scheidt et al. 2019 ;Wölfer et al. 2019 ;Etienne et al. 2021 ). Here, we found that these factors exhibit different relationships with the external shape and structure of squirrel humeri and femora: ecotype influenced variation of the external shape of these bones, whereas interactions between size and ecotype influenced variation of their structure.
Our results suggest that differences in locomotor behavior have a stronger effect on humeral shape evolution and, to a lesser extent, femoral shape evolution rather than size (contrary to our predictions; Table 1 ). Using our model selection approach, we found that humeral shape is best described by the mvPGLS ecotype model ( Table 1 ), which revealed a significant relationship between humeral shape and ecotype. Gliding squirrels tend to exhibit gracile, elongate humeri; ground squirrels tend to exhibit robust humeri; and chipmunks and tree squirrels tend to exhibit intermediate humeri ( Fig. 3 ; Table 1 ). However, the statistical relationship between ecotype and humeral shape was lost when using a model that assumes Brownian motion. In contrast, we found that femoral shape is best described by the mvPGLS null model ( Table 1 ), suggesting that femoral shape is not well predicted by ecotype, size, or their interaction. Nevertheless, the CVA indicated that femoral shape can be accurately reclassified by ecotype, and the Procrustes ANOVA also supported a significant relationship between femoral shape and ecotype. However, like in the humeral shape analyses, the statistical relationship between femoral shape and ecotype is lost when incorporating phylogenetic signal as Pagel's λ or full Brownian motion (i.e., λ = 1). That accounting for phylogenetic relationships confounded these relationships is not surprising considering squirrel ecotypes are phylogenetically clustered ( Fig. 3 ); chipmunks and gliding squirrels are each monophyletic, ground squirrels consist of two clades, and tree squirrels are ancestral ( Fig. 1 ). These results also indicate that there was no evolutionary transition from relatively elongate, gracile humeri and femora to relatively shorter, robust ones with increasing size across squirrels, a pattern that is found in other mammalian clades (e.g., Fabre et al. 2013 ;Martín-Serra et al. 2014a, 2014bEtienne et al. 2021 ). Instead, the lack of shape allometry in the proximal limb bones is consistent with what is found in myomorph and geomyoid rodents ( Hedrick et al. 2020 ). This provides additional evidence that only large mammals require major allometry-driven shape changes ( Biewener 2005 ) or allometric changes occur in other aspects of the limb morphology such as bone structure (see next section). Overall, our results emphasize that locomotor behavior has a stronger influence on humeral shape compared to femoral shape, and that the diversity of squirrel humeral shapes and to a lesser extent, femoral shapes, were ecologically partitioned early between clades and these ecomorphologies were maintained to the present.
In contrast, bone structure variables of the humerus and femur exhibited various scaling trends across all squirrels. Global compactness (Cg) of both the humerus and femur was best described by the PGLS size model ( Table 1 ), which shows Cg increased with increasing limb size ( Fig. 5 ). This pattern is consistent with findings that cortical bone thickness ( Currey and Alexander 1985 ) and cross-sectional properties of the femur ( Mielke et al. 2018 ;Scheidt et al. 2019 ) scale with increasing size. These positive allometric relationships may be due to the increase in mechanical stress at larger body sizes ( McMahon 1973 ;Biewener 1990 ;Christiansen 1999 ). Positive allometry towards more compact bones may also free the constraints of allometry-driven changes in limb bone shape; the mechanical demands against gravity may be met with positive allometry of Cg, enabling other ecological factors to have greater influence on limb bone shapes. CSS of the humerus and femur was best described by models that included both size and ecotype ( Table 1 ). These models revealed that scaling patterns differed between humeral CSS and femoral CSS; humeral cross sections became more ovalshaped as size increased through positive allometry, whereas femoral CSS exhibited isometry. DE of both the humerus and femur was best described by models that included both size and ecotype and their interaction ( Table 1 ). These models revealed that scaling patterns of humeral and femoral DE did not statistically differ from isometry across all squirrels. Like the bone shape analyses, the statistical relationships between bone structure variables and ecotype were lost when using phylogenetic ANOVAs that assume Brownian motion but not with the non-phylogenetic ANOVAs. These results are consistent with the bone shape analyses, providing additional support that humeral and femoral morphologies were ecologically partitioned early between clades and their ecomorphologies were maintained to the present.

Scaling and ecomorphology of limb bones within ecotypes
Scaling patterns in external shape and bone structure variables of the humerus and femur are further nuanced by ecological specialization. In ground squirrels, positive allometry towards more robust humeral shapes is consistent with positive allometry of more robust body shapes ( Linden et al. 2023 ), which together may serve as adaptations for the digging and clearing stages of scratch digging during burrow construction. The digging stage consists of the claws of the forelimbs striking the ground while the hind limbs support the body ( Gasc et al. 1985 ;Hildebrand 1995 ), and a plethora of skeletal and muscular adaptations in the forelimb are associated with this first stage ( Lessa and Stein 1992 ;Hildebrand 1995 ;Lagaria and Youlatos 2006 ;Vassallo 2006 ;Steiner-Souza et al. 2010 ). Our finding of positive allometry in humeral shape and CSS provides additional evidence for adaptation towards digging, particularly in larger ground squirrels with more elliptical cross sections that would better resist uniaxial bending loads during digging. Trends towards negative allometry of humeral DE would also lead to more compact, robust humeri that could decrease bone strain during the digging stage. In contrast, humeral Cg scaled isometrically with size, indicating that compact humeri are important for digging across all ground squirrel sizes. These patterns are also consistent with the analyses of humeral shape and its phylomorphospace. Ground squirrel humeri occupied regions of phylomorphospace ( + PC1) that are associated with more pronounced deltoid tuberosities ( Fig. 3 A) that would provide greater attachment surface area for deltoids and pectoral muscles, which are needed to increase the force output advantageous for digging ( Hildebrand 1985 ;Samuels and Van Valkenburgh 2008 ). These humeri associated with + PC1 also exhibited more pronounced medial epicondyles that would provide greater attachment sites for the extensor group key for scratch digging ( Samuels and Van Valkenburgh 2008 ;Steiner-Souza et al. 2010 ).
Unlike in the humerus, our results for the relationships between femoral morphology and digging behavior are not as consistent. The clearing stage of scratch digging consists of the hind limb removing excess soil through backward extension while the body is supported by the forelimbs ( Gasc et al. 1985 ;Hildebrand 1995 ). Whereas some have hypothesized that more robust femora in fossorial species would improve stability and load transfer during clearing ( Casinos et al. 1993 ;Hildebrand 1995 ), others have hypothesized that fossorial species exhibit less robust femora to reflect less rigidity in bending and torsion compared to species that jump or run with their hind limbs ( Gambaryan 1974 ;Biknevicius 1993 ;Wölfer et al. 2019 ). Furthermore, some have also hypothesized that robust femora are not a necessary adaptation because the hind limbs are only used to shovel away excess substrate during the clearing stage that had previously been loosened with the forelimbs during the digging stage ( Gambaryan 1974 ;Wölfer et al. 2019 ). Our findings support both these hypotheses: positive allometry of femoral Cg and negative allometry of DE indicate that ground squirrel femora evolved to be more compact and robust with increasing femoral size, which may decrease bone deformation during the clearing phase. Conversely, the lack of allometry in femoral shape and CSS suggests that the femur is less specialized for digging compared to the humerus. Ground squirrel femora occupied regions of phylomorphospace ( + PC1) ( Fig. 3 B) that are associated with more pronounced greater trochanters. These adaptations reflect the larger muscle attachment area for the gluteal muscles, which when enlarged in fossorial species may help the body resist being pushed back when digging, respectively ( Samuels and Van Valkenburgh 2008 ).
Tree squirrels exhibited positive allometry in humeral and femoral Cg with increasing bone size. Additionally, we found that tree squirrels evolved more robust diaphyses (i.e., lower DE) in both the humerus and femur with increasing bone size. This contrasts with prior findings that scansorial or arboreal mammals often exhibit more elongate limb bones as size increases to increase limb span ( Polly 2007 ;Kilbourne and Hoffman 2013 ). However, this may be due to the somewhat unique climbing method employed by tree squirrels, which use their claws to interlock into tree bark and their muscle strength to hold their bodies close to the tree as they climb ( Schmidt and Fischer 2011 ). Furthermore, allometric trends towards relatively shorter humeri may increase mechanical advantage during climbing, given that the elbow functions like a class II lever ( Keith 1919 ;Fleagle 1977 ). At larger body sizes, tree squirrels may need to counter the additional body weight by remaining relatively closer to the tree to minimize opposing forces. No significant relationships were observed between either humeral or femoral CSS and bone size in tree squirrels. This may point to more circular cross-sectional shapes across different sizes as important in resisting multidirectional bending while climbing and jumping between trees ( Patel et al. 2013 ).
The external shape and structure of chipmunk humeri and femora were all isometric. The lack of allometry in chipmunks may be due to their small sizes and limited size range. Prior research in other small mammals has shown that the ability of bone tissue to remodel in response to its mechanical environment (i.e., Wolff's Law) may be restricted by body size ( Dawson 1980 ;Meier et al. 2013 ). This makes specialization of bone structure traits such as compactness less likely to occur, as overall humeral morphology may be sufficient for digging or climbing without further cortical adaptations. Another explanation is that chipmunks primarily dig shallow shelters instead of extensive burrow systems and thus may not require extreme shape specializations ( Samuels and Van Valkenburgh 2008 ). That chipmunk humeri and femora share overlapping regions of phylomorphospace with both tree and ground squirrels ( Fig. 3 ) support these hypotheses.
Gliding squirrels evolved more gracile humeral and femoral shapes with increasing size, but, interestingly, humeral and femoral DE were isometric. The lack of DE allometry in gliding squirrels may be due to wing loading (body mass/patagium area); larger gliding squirrels decrease the detrimental effects of higher wing loading by exhibiting a lower relative body weight instead of evolving relatively longer limbs compared to smaller gliding squirrels ( Thorington and Heaney 1981 ). Nevertheless, gliding squirrels on average exhibited more elongate humeri and femora compared to other ecotypes ( Figs 3 and 6 ), consistent with previous findings ( Peterka 1936 ;Bryant 1945 ;Thorington and Heaney 1981 ;Grossnickle et al. 2020 ;Linden et al. 2023 ). Furthermore, gliding squirrels tended to exhibit the most circular cross sections of the humerus and femur (i.e., low CSS values; Fig. 6 ), which would maximize resistance of torsional stresses during aerial locomotion ( Swartz et al. 1992 ) or like tree squirrels multidirectional bending while climbing and jumping between trees ( Patel et al. 2013 ).
Comparing humeral and femoral shape and bone structure Humeral shape and femoral shape are tightly correlated across squirrels (r-PLS = 0.91, Z = 5.884, P < 0.001), a pattern that is consistent with other mammals (e.g., Martín-Serra et al. 2015 ;Hanot et al. 2017 ;Hedrick et al. 2020 ). Together, this indicates that the relative robustness of the humerus and femur contributes the greatest covariation between the two bones. However, we found that these patterns differed between scansorial and fossorial ecotypes: while gliding squirrels and tree squirrels exhibit correlated humeral and femoral shapes, chipmunks and ground squirrels do not. These differences may not be surprising because the humerus and femur exhibit different functions during the digging and clearing stages, respectively, of burrow construction (see previous section).
Ecotype may also influence whether bone structure traits differ between the humerus and femur. Previous studies suggest that scansorial mammals exhibit similar global compactness values between the forelimb and the hindlimb due to the equal importance of both limbs during climbing, whereas scratch-digging mammals exhibit higher compactness in the forelimb due to the many specializations associated with digging ( Hildebrand 1985 ;Straehl et al. 2013 ). We found higher mean global compactness in humeri than in femora in chipmunks and ground squirrels, an expected result since chipmunks and ground squirrels exhibit forelimb specializations associated with a fossorial lifestyle ( Lessa and Stein 1992 ;Hildebrand 1995 ;Lagaria and Youlatos 2006 ;Vassallo 2006 ;Steiner-Souza et al. 2010 ). Surprisingly, tree squirrels also exhibited higher mean global compactness in humeri than in femora. A possible explanation is that although tree squirrels are scansorial, they are not limited to only climbing. Tree squirrels are often found running on the forest floor, and often bury food by scratch-digging ( Ferron 1981 ). This may help to explain some of the adaptions seen in tree squirrels that have been historically correlated with fossorial locomotion. We also found higher DE values in the femur compared to the humerus in each of the four ecotypes. More elongate femora support previous findings that rodent humeri tend to be the shortest long bone whereas rodent femora are typically one of the longest ( Prodinger et al. 2018 ). Lastly, CSS values do not differ between the humerus and femur, which is consistent with findings that humeral and femoral shape covary along the degree of bone robustness.

Conclusion
We found that size, locomotor ecology, and underlying phylogenetic structure exhibit different relationships with the shape and structure of squirrel humeri and femora. Humeral and femoral shape is best predicted by ecotype, whereas humeral and femoral structure is best predicted by a combination of ecotype, size, and their interaction. Interestingly, the statistical relationships between these morphologies and ecotype were lost when accounting for phylogenetic relationships under Brownian motion. That all squirrel ecotypes are largely phylogenetically clustered invites questions about whether and how morphological variation in the humerus and femur evolved in relation to different ecologies early in squirrel evolutionary history and why they were subsequently maintained within clades. Clade-based evolutionary shifts in morphologies have been found in other mammals such as carnivorans ( Law 2021 ;Law et al. 2022 ), and our results provide preliminary evidence that squirrels also exhibit clade-based evolutionary shifts in morphologies. Overall, these results indicate that mechanical and phylogenetic constraints and ecology may enact different pressures on the external and structural aspects of limb bone morphology. Increasing the number of species across the squirrel phylogeny and expanding to other rodents will enable us to tease apart the effects of selection and phylogeny on potentially converging limb morphologies. Together, this work provides a strong morphological foundation for future research investigating the evolutionary biomechanics and ecology of squirrel locomotion.