Weighing homoplasy against alternative scenarios with the help of macroevolutionary modeling: A case study on limb bones of fossorial sciuromorph rodents

Abstract Homoplasy is a strong indicator of a phenotypic trait's adaptive significance when it can be linked to a similar function. We assessed homoplasy in functionally relevant scapular and femoral traits of Marmotini and Xerini, two sciuromorph rodent clades that independently acquired a fossorial lifestyle from an arboreal ancestor. We studied 125 species in the scapular dataset and 123 species in the femoral dataset. Pairwise evolutionary model comparison was used to evaluate whether homoplasy of trait optima is more likely than other plausible scenarios. The most likely trend of trait evolution among all traits was assessed via likelihood scoring of all considered models. The homoplasy hypothesis could never be confirmed as the single most likely model. Regarding likelihood scoring, scapular traits most frequently did not differ among Marmotini, Xerini, and arboreal species. For the majority of femoral traits, results indicate that Marmotini, but not Xerini, evolved away from the ancestral arboreal condition. We conclude on the basis of the scapular results that the forelimbs of fossorial and arboreal sciuromorphs share mostly similar functional demands, whereas the results on the femur indicate that the hind limb morphology is less constrained, perhaps depending on the specific fossorial habitat.


Some of the most instructive examples of homoplasy stem from
the limbs of tetrapods. For example, various mammalian lineages that independently acquired a fully or partly subterranean lifestyle also independently evolved a larger in-lever to out-lever ratio for muscles which retract the arm (teres major, and latissimus dorsi) and/ or extend the forearm (triceps) as compared to cursorial relatives (Hildebrand & Goslow, 1995). This increases the muscles' force output, which was interpreted to be advantageous for digging activities (Hildebrand & Goslow, 1995). Samuels and Van Valkenburgh (2008) could demonstrate that scratch-digging rodents in distantly related taxa also display similarities in the hind limb anatomy. According to them, one shared similarity as opposed to their terrestrial relatives constitutes a larger ratio between the height of the greater trochanter on the femur and the length of the femur and thus a larger mechanical advantage of the gluteus muscles that attach to this trochanter and retract the hind limb. This was interpreted to assist the stabilization of the body against forces produced by the forelimb during digging (Samuels & Van Valkenburgh, 2008).
These examples demonstrate that homoplasy in limb traits of fossorial mammals can be a particularly revealing topic for investigation to foster the understanding of how form reflects function.
Many similar studies have been conducted on fossorial mammals or at least included taxa with independent acquisitions of a fossorial lifestyle (e.g., Carrizo, Tulli, Dos Santos, & Abdala, 2014;Hildebrand & Goslow, 1995;Hopkins & Davis, 2009;Lehmann, 1963;Lessa, Vassallo, Verzi, & Mora, 2008;Morgan, 2009;Morgan & Álvarez, 2013;Piras et al., 2012;Samuels & Van Valkenburgh, 2008;Stein, 2000;Warburton, Grégoire, Jacques, & Flandrin, 2013). However, only a few have been put into the framework of phylogenetic comparative methods (PCMs) with a focus on homoplasy in fossorial taxa (e.g., Meier, Bickelmann, Scheyer, Koyabu, & Sánchez-Villagra, 2013;Piras et al., 2012). One important aspect of a phylogenetic framework is to account for the directionality of evolutionary transformations (Cracraft, 1980). This is important, because examples of homoplasy can only be inferred if the supposedly homoplastic character states are indeed identified as independently evolved novel character states in all focal groups or as at least one independent acquisition of the ancestral character state in case of a reversal. In an adaptive context, this has to be equally true for form and function (Losos, 2011). PCMs represent a class of inferential statistics that utilize this phylogenetic framework to evaluate whether the species' traits within a clade of interest have evolved by chance or resulted from adaptive evolution (Felsenstein, 1985;Hansen, 1997;Harvey & Pagel, 1991;. However, only a few PCMs have been developed to investigate the adaptive significance of homoplasy (reviewed in Stayton, 2015; and see Khabbazian, Kriebel, Rohe, & Ane, 2016 for a more recent method).
The most suitable PCM to assess homoplasy available to date, in our opinion, is the likelihood comparison among a set of macroevolutionary models, because one can compare homoplasy hypotheses to multiple alternative plausible scenarios. A frequently used model is the Ornstein-Uhlenbeck (OU) model which was introduced into a macroevolutionary context by Hansen (1997). The OU model assumes that a trait evolves toward a hypothetical optimal state θ that is best suited to accomplish the functions imposed by the selective regime with the highest selective pressure on that trait. Trait evolution is thereby driven by an adaptive rate α and random perturbations σ which prevent the species to reach the optimum as a result of less influential selective factors and historical constraints (Hansen, 1997). The Brownian motion (BM) model, which was introduced into a macroevolutionary context by Felsenstein (1985), can be regarded as a special case of the OU model with α being zero and hence reflecting a random walk-like trait evolution. The likelihood of competing models can then be compared as suggested by Butler and King (2004).
A very suitable design using macroevolutionary model comparison to evaluate homoplasy in the context of adaptation is presented by Moen, Morlon, and Wiens (2015). They compared the likelihood of frog postcranial morphology reflecting either the independent acquisition of lifestyles, phylogenetic history, or a combination of both factors. They split the selective regimes (lifestyle groups) in order to evaluate the frequency of trait homoplasy on different levels of the phylogeny (homoplasy across the whole phylogeny vs. homoplasy only within biogeographical separated subclades vs. all independent lifestyle acquisitions with different optima). It needs to be emphasized that trait homoplasy in the context of OU model comparison is evaluated in respect to the trait's optima (i.e., trait optimum homoplasy) and not to the species' trait values themselves (i.e., trait homoplasy). Besides models with different levels of homoplasy, Moen et al. (2015) accounted for further plausible scenarios by including models into the comparison that represent nonadaptive evolution or evolution of all species toward a common optimum.
Here, we used this approach by Moen et al. (2015) to evaluate the likelihood of postcranial homoplasy in the limb morphology of extant sciuromorph rodent lineages that independently acquired a fossorial lifestyle. We complemented their approach by the method by Boettiger, Coop, and Ralph (2012) for measuring statistical power.
We only included two of the three fossorial sciuromorph lineages, Marmotini and Xerini, because the third lineage (Aplodontia rufa) is monospecific and we do not consider it meaningful to estimate an optimum trait value using a single species. Marmotini is spread across the Holarctic region, and Xerini is present in Africa and Central Asia (Nowak, 1999;Thorington Jr., Koprowski, Steele, & Whatton, 2012). Species of both taxa dig underground burrows, but search for food above ground (Nowak, 1999;Thorington Jr. et al., 2012).
As homoplasy in the morphology of the skeletal limb elements of fossorial mammalian taxa has been previously demonstrated in both, fore limbs and hind limbs (see above), we compared homoplastic tendencies between the scapula and the femur, (i.e., the most proximal functional skeletal elements of therian limbs, respectively; Fischer, 1994;Fischer & Blickhan, 2006;Fischer, Schilling, Schmidt, Haarhaus, & Witte, 2002;Kuznetsov, 1985). The scapula and femur were also chosen, because they exhibit structures that allow us to derive muscle properties from the morphology of the bone. We analyzed various univariate traits reflecting length, robustness, and the ability to transmit forces by the attaching muscles. This gave us an idea about how these traits conform or differ in the most likely process of trait evolution.
The here-investigated traits of the sciuromorph scapula and femur were recently investigated in terms of how they reflect differences in lifestyle and body mass by  and Wölfer, Amson, et al. (2019), respectively. Both studies revealed that some traits differ significantly between the arboreal and fossorial species. In case of such a significant difference, the trait value of the fossorial group (including A. rufa, Marmotini, and Xerini) was always smaller than that of the arboreal group. In both cases, the authors assumed that this reflects lower demands of a fossorial lifestyle regarding running velocity (a relatively shorter femur leads to shorter limb length), bone robustness, and muscle force output. Both of these studies found the phylogenetic inertia to be low in all traits, suggesting differences to bear an adaptive signal (Wölfer, Amson, et al., 2019;. This leads to the question as to whether this hypothesized adaptedness of proximal limb bone traits in Marmotini and Xerini is reflected in homoplasy. Especially, the differences in clade size of these two tribes render it problematic to infer homoplasy simply from a significant difference between the fossorial and arboreal groups. Marmotini with approx. Two previous studies (Mielke et al., 2018;Scheidt, Wölfer, & Nyakatura, 2019) analyzed further robustness parameters of the sciuromorph femur (trabecular parameter of the femoral head and cross-sectional properties along the proximodistal axis, respectively). Mielke et al. (2018) found no significant differences between the arboreal and fossorial groups. Scheidt et al. (2019) demonstrated that fossorial sciuromorphs are more robust in the distal epiphyseal region than their arboreal relatives. We did not consider the femoral traits investigated by Mielke et al. (2018) and Scheidt et al. (2019), as only four of six xerine species were sampled in both studies what we considered too few to reliably estimate trait optima. Moreover, the dataset of Scheidt et al. (2019) is too large (i.e., 57 variables) for the computationally heavy methods applied herein.
We first explain how the functionally relevant univariate traits were acquired. This is followed by a subtraction of the statistical effect of body mass on those traits. Then, we explain the considered models and how we compare their likelihood to evaluate overall macroevolutionary patterns. We also integrated a simulation study to assess the effect of our sampling design on the likelihood comparison. Finally, in case likelihood comparison suggests a potential case of homoplasy, we utilized pairwise evolutionary model comparison sensu Boettiger et al. (2012) to assess whether the data and phylogeny are powerful enough to unambiguously support this hypothesis in favor of all other considered models.

| Software
Functions from the software R version 3.5.2 (R Development Core

| Phylogenetic tree
We used the compound phylogenetic tree that was assembled by Wölfer, Amson, et al. (2019); displayed in their Supporting Information) as the raw, unpruned phylogeny. It consists of extant species from the maximum credibility consensus tree presented by Zelditch et al. (2015) and from the TimeTree database (Hedges, Marin, Suleski, Paymer, & Kumar, 2015).

| Morphological traits
Our investigation is based on the scapular dataset acquired by  and the femoral dataset acquired by Wölfer, Amson, et al. (2019). The former study used regressions to analyze 15 univariate scapular traits as well as scapular shape (i.e., a multivariate trait) in regard to an allometric effect and whether this effect differs depending on locomotor ecology (termed interaction effect). Wölfer, Amson, et al. (2019) used the same approach to study 11 univariate femoral traits as well as femoral shape. Here, we only investigated univariate traits of both skeletal elements, because likelihood estimation is unreliable when fitting complex multivariate models (Adams & Collyer, 2017).  (Tables S1 and S2, respectively, in Supplementary Material). Specimens of different sex were included depending on availability (Tables S1 and S2). The level of maturity was assessed by comparing the sizes of specimens available at the collection sites and always selecting one of the larger bones. The fusion of the epiphyses was not used as criterion, as this was not observed to be the rule even in the largest specimens. Thus, a sensitivity analysis was conducted to assess the influence of sampling bias on evolutionary model comparison (see below).
According to  and Wölfer, Amson, et al. (2019), four univariate traits of the scapula and the femur, respectively, displayed a significant interaction effect between lifestyle and body mass. This interaction effect might influence the differences among groups if these differ in their body mass. Thus, it would not be possible to discriminate between trait homoplasy between Marmotini and Xerini as a result of independent lifestyle acquisitions and trait similarity as a consequence of the interactive influence of body mass and lifestyle. We tested whether the body mass ranges of the three groups differed significantly (S1). As this was the case between each fossorial group and the arboreal group (Tables S3 and   S4), we decided to only include those traits here that did not display a significantly interaction effect that resulted in a different slope between the fossorial and arboreal groups in the previous studies by  and Wölfer, Amson, et al. (2019).
F I G U R E 1 Investigated traits of the scapula and femur. A: Lateral (a), medial (b), caudal (c), and ventral (d) views on the scapula. Red circles = fixed landmarks. Dark gray circles = curve landmarks that were evenly spaced between the fixed landmarks to capture the outlines of the scapula. Light gray circles = either centroids of traits 6, 8, and 10 (b) or points used to compute trait 4 (c, d). B: Caudal (a), medial (b), dorsocaudal (c), and ventrocranial (d) views on the femur. Red circles = fixed landmarks. Dark gray circles = sliding semi-landmarks. Light gray circle = centroid to compute trait 1. apD, anteroposterior diameter; ca, caudal view; CS, centroid size; EL, effective length; I-L, inlever length; ml, mediolateral; mlD, mediolateral diameter; W, width The univariate traits analyzed here cover various aspects of the scapula and the femur, such as the effective length, the robustness of the articulation sites, and the muscle properties that can be derived from the morphology of the bone (the sizes of the attachment sites as well as the lengths of the in-levers). They will be briefly described in the following, but see  and Wölfer, Amson, et al. (2019) for a detailed description of the acquisition of those traits.
The scapular traits were extracted from landmarks that were placed onto images obtained from four different perspectives orthogonal to each other (lateral, medial, caudal, and ventral; The femoral traits were obtained by Wölfer, Amson, et al. (2019) from 3D surface models. Of those that were used here, four were derived from landmark data and three were measured directly on the surface scans ( Figure 1). Five robustness traits were included. The centroid size of the head was computed using the landmarks surrounding it. The anteroposterior and mediolateral midshaft diameters and the width of the patellar groove were measured directly on the surface model. The widths of the medial and lateral condyles were computed as the distance of the respective most medial and most lateral fixed landmark on the proximal site of each condyle. The only in-lever that did not display an interaction effect concerned the muscles attaching to the lesser trochanter (most importantly the iliopsoas muscle that protracts and adducts the hind limb by flexing the hip joint).
We used the unpruned datasets (186 and 177 species for the scapula and femur, respectively) to remove trait differences that were statistically linked with differences in body mass. Body mass information was taken from the literature (Tables S1 and S2). The lm.rrpp function of the package "RRPP"  with the option SS.type = "I" and 10,000 rounds of permutation was applied to obtain the residuals for the natural log-transformed traits on natural log-transformed body mass. We did not use phylogenetic correction for the estimation of the regression parameters, as it was shown before that phylogenetic inertia is negligible in the traits studied here (Wölfer, Amson, et al., 2019;. The residual datasets were pruned to match the species composition of the pruned phylogeny. Regarding the scapula, 76 arboreal species (representing the monophyletic groups Callosciurinae, Gliridae, Protoxerini, Ratufa, Sciurini, Sciurotamias, and Tamiini), 43 marmotine species, and all 6 extant species of Xerini were investigated (i.e., 125 species overall; Table S1; Figure A1). For the femur, we were able to include 74 arboreal, 43 marmotine, and the six xerine species (i.e., 123 species overall; Table S2; Figure A1). Violin plots including the arithmetic mean were generated to provide a graphical comparison of the distributions of the residuals among the three groups.

| Model definition
The OU2 foss model represents the hypothesis of trait optimum homoplasy between the two fossorial clades, Marmotini and Xerini ( Figure 2

| Ancestral state reconstruction
The selective regimes had to be reconstructed on the two pruned trees for four of the six investigated models in order to evaluate their likelihood (OU2 foss , OU2 arb&Ma , OU2 arb&Xe , OU3; Figure 2). For all four models, we assumed specific changes in selective regimes to occur at the root of the stem lineage leading to Xerini and/or Marmotini ( Figure 2). This assumption is common for the method of "tree painting" proposed by Butler and King (2004). The paintSubTree and paintBranches functions of the package "phytools" (Revell, 2012) were used to define the regimes.

| Evolutionary model comparison via likelihood scoring
We compared the likelihood of all six models for each trait using the Schwarz information criterion (SIC) that penalizes the likelihood of a model with the number of parameters estimated and the number of species included (Burnham & Anderson, 2002). Lower SIC values were interpreted as higher model likelihood (Burnham & Anderson, 2002 (Tables S5 and S6).
This is defined as half the evolutionary time necessary for a lineage to switch to another optimum when entering a novel selective regime (Hansen, 1997) and was used to discuss historical constraints on trait evolution.
We combined the SIC scoring with a sampling simulation to evaluate whether the sampling of a single individual per species affects the outcome of the ranking (S2). In short, we assessed the intraspecific standard deviation of selected species and used them to generate hypothetical populations for all species. From those, we draw samples including a single specimen per species 1,000 times and redid the analysis each time. Finally, it was compared whether the model with the highest frequency of displaying the lowest SIC score among the simulated datasets was also the model that was demonstrated to be the most likely for the original dataset. If this was not the case, the trait was dismissed from further analysis. Accordingly, two traits were removed (the centroid sizes of the teres major and supraspinatus fossae; see Tables 1 and 2).

0.3
Note: Value in bold: lowest Schwartz information criterion score per trait (empirical dataset) or highest frequency (simulated datasets). See Figure 1 for trait abbreviations and text for explanation of model abbreviations. For each pairwise model comparison (OU2 foss vs. one of the five alternative models), the following procedure was applied. Firstly, the saved likelihood values of the model fits outlined above were used to compute the likelihood ratio between both models, an expression of the difference between the two likelihood values of the two models.

| Pairwise evolutionary model comparison
Here, we refer to it as the empirical likelihood ratio δ emp . Secondly, the saved parameters of the two models were used to simulate two distributions of data (5,000 datasets) along the phylogenetic tree that could have evolved according to the two models. This ensures that the structure of the phylogenetic tree is accounted for when assessing the likelihood between the two models. For this purpose, the function mvSim of the package "mvMORPH" was used. Finally, both models were fit to all datasets of both distributions using the functions and options outlined for the evolutionary model comparison via likelihood scoring. This resulted in four distributions of parameter and likelihood values. A likelihood ratio distribution was computed for the two likelihood distributions that were fit to the same simulated data distribution. This lead to two likelihood ratio distributions (one belonging to each of the two models) and allowed us to evaluate which of the two likelihood ratio distributions better reflects δ emp and hence which of the two models better fits the data. According to Boettiger et al. (2012), one way to do so is to check whether δ emp lies within the 95% confidence interval (CI) of the respective likelihood ratio distribution. We encountered three scenarios: either δ emp lay only inside the 95% CI of the distribution generated from the OU2 foss model, inside the 95% CIs of the distributions generated from both models, or outside of the 95% CIs of the distributions generated from both models (exemplified in Figure 3). Only the first-mentioned case ( Figure 3B) would have indicated that the distribution of the residuals across species and the structure of the phylogenetic tree are powerful enough to discriminate between the likelihood of the two models. The latter two cases would have not allowed us to discriminate between the two models. Thus, we decided that the homoplasy hypothesis can only be considered favorable in the first case. The  (Tables S5 and S6).

| Evolutionary model comparison via likelihood scoring
Only three of six models (OU1, OU2 foss , or OU2 arb&Xe ) scored lowest according to the SIC (Tables 1 and 2)

48.0
Note: Value in bold: lowest Schwartz information criterion score per trait (empirical dataset) or highest frequency (simulated datasets). See Figure 1 for trait abbreviations and text for explanation of model abbreviations.
parameters were reliable (i.e., they fell within their 95% CI), except for the nonadaptive evolutionary rate σ concerning the BM1 model (Tables   S5 and S6). This rate was not reliably estimated for any trait that we studied. The largest estimate for t 1/2 among all traits and all OU models fit to them was 9.5 million years. This means that the switch from one optimum to another took at most ~19 million years. According to Regarding the scapular traits, the OU1 model most frequently yielded the lowest SIC values (Table 1). This was the case for the effective length, the size of the glenoid cavity from the caudal view, and the lengths of the in-levers of the rotator cuff muscles.
The OU2 foss model displayed the lowest SIC score for the size of the coracoid process and the mediolateral in-lever of the muscles attaching to this process' tip. The OU2 arb&Xe model yielded the lowest SIC score for the two centroid sizes of the infraspinatus and subscapularis muscle attachment sites.
Regarding the femoral traits, the OU2 arb&Xe model most frequently yielded the lowest SIC values (Table 2). This was the case for the centroid size of the head, the two midshaft diameters, the width of the patellar groove, and the in-lever of the muscles attaching to the lesser trochanter. The OU1 model was the most likely one for the width of the medial condyle, and the OU2 foss model was the most likely choice for the width of the lateral condyle. The violin plots appear to reflect the modeling results ( Figures 4 and 5).

| Pairwise evolutionary model comparison
According to the previous steps of the analysis, only two scapular traits (the centroid size of the coracoid process, the mediolateral inlever of the muscles attaching to the process' tip) and one femoral trait (the width of the lateral condyle) were qualified for the investigation using the power analysis suggested by Boettiger et al. (2012). For none of these three traits, the suggested hypothesis of trait optimum homoplasy (OU2 foss ) was favoured over all the other models of trait evolution, that is, at least one alternative model could not be discriminated from the homoplasy hypothesis in terms of model likelihood (Table 3).

| Historical constraints on homoplasy
It has been noted by previous authors that the skeletal morphology of sciuromorph clades is much conserved (Black, 1963;Bryant, 1945;Cardini, Hoffmann, & Thorington, 2005;Cardini & O'Higgins, 2004;Emry & Thorington Jr., 1982;Moore, 1959) or reflecting phylogenic relationships rather than ecological differences (Cardini, 2003). In our study, the role of historical constraints limiting the realization of homoplastic trait optima can be generally questioned on the basis of the low likelihood of the BM1 model for all traits (Tables 1 and   2; also reflected in the bad estimation of σ according to the CIs) and the estimated phylogenetic half-life t 1/2 for all OU models (assuming that at least one OU model considered in our study is somewhat representative of the true underlying evolutionary process). The longest switch from one optimum to another took at most ~19 million years among all traits and all OU models, which is shorter than the time since Marmotini and Xerini acquired a fossorial lifestyle. Thus, all three groups (including the arboreal group, as it is the oldest one) should have already reached their optimum, independent of which OU model was most likely. Consequently, we have to assume that adaptation has played a major role in the evolution of the proximal limb elements, but did not result in homoplasy. Our sampling method (one individual per species) did only bias the modeling outcome of two traits that we consequently dismissed from further consideration. Still, for the three traits that qualified for the power analysis, the number of species per group might have favored the rejection of the homoplasy hypothesis. Although we already sampled all six extant xerine species, this number might still be too low to statisitcally assess the similarity between their trait optima on those of Marmotini.
Including extinct species, which were not available for the current analysis, could help to achieve a better assessment of the dynamics of shape optimum shifts in future studies (Cooper, Thomas, Venditti, Meade, & Freckleton, 2016;Mahler, Weber, Wagner, & Ingram, 2017;Slater, Harmon, & Alfaro, 2012).

| Trends in trait evolution according to likelihood scoring
As the OU1 model displayed the lowest SIC value regarding five out of nine scapular traits, we assume that the requirements on this skeletal element to resist and transmit forces are very similar between the fossorial and arboreal species. This is in agreement with Stalheim-Smith (1984), who demonstrated on the basis of comparative physiological experiments on selected forelimb muscles of an arboreal and a fossorial sciuromorph rodent species that both have similar force-generating potentials. In the few instances, in which differences were found, the arboreal species exceeded the fossorial one in its potential (Stalheim-Smith, 1984). This also reflects the fact that the two fossorial lineages only evolved toward smaller trait optima in case they departed from the arboreal optimum. The only scapular traits in which both, Marmotini and Xerini, likely departed from the arboreal optimum concerned the two traits associated with the coracoid process. A longer mediolateral in-lever of the muscles attaching to it allows for the generation of larger adduction forces, which may not be as important during digging as compared to climbing. Similarly, the strength of the articulation created by the coracoid process and the clavicula as determined by the size of the process appears to be of less importance. We think that this can be related to the relevance of mediolateral movements of the forelimb which are more likely to be required during climbing than digging.
Since the OU2 arb&Ma model never scored lowest and the OU2 arb&Xe model was most frequently suggested to be favorable among all traits according to the SIC (seven out of 16 cases), Marmotini rather than Xerini evolved away from the ancestral arboreal trait optimum. In five of these cases, it concerned femoral traits, indicating that differences in trait evolution between Marmotini and Xerini are more frequently found in the hind limb. As Marmotini displayed a smaller trait optimum value compared to Xerini and the arboreal group, this might hint at a relaxed functional constraint on the hind limb concerning its ability to resist and transmit forces.
Perhaps this is related to differences in the utilization of the hind limb during digging between these two fossorial groups, leading to differing requirements on this skeletal element. Differences in how the hind limb is used (e.g., stabilization of the body during forelimb digging, hind limb digging, or pushing the body forwards when using the forelimbs to shovel away the soil) have been observed in other fossorial rodent taxa (Lehmann, 1963). Gambaryan (1974) argued that kicking the loose dirt is not a very strenuous behavior, but the forces which need to be generated during hind limb digging or resisted when the forelimbs are utilized might depend on the characteristics of the habitat. For example, it was shown that the soil hardness is reflected in the humeral morphology of different populations of the caviomorph rodent species Ctenomys minutus (Kubiak et al., 2018). Xerine species are found in arid habitats in Africa and Central Asia, whereas Marmotini is spread across the Holarctic region (Nowak, 1999;Thorington Jr. et al., 2012).

| CON CLUS IONS
The independent acquisition of a fossorial lifestyle in mammals is a textbook example of homoplastic evolution of the postcranium (Hildebrand & Goslow, 1995). This is not supported by our study that utilized a rigorous phylogenetic comparative framework to study the most proximal limb bones of sciuromorph rodents. Our findings rather indicate a mosaic pattern of evolutionary processes depending on the skeletal element and trait under consideration. We provide evidence that the scapula tends to be constrained in more traits than the femur. This appears to be a result of shared functional demands between fossorial and arboreal species and not because of phylogenetic constraints. Contrary to this, the femur might be less constrained in its ability to resist and transmit forces. This suggests that its morphology depends on the requirements of the specific fossorial habitat, perhaps manifested in the biogeographical differences between Marmotini and Xerini. Comparative studies integrating habitat structure, behavior, performance, and biomechanics (as suggested by, e.g., Bock, 1980;Losos, 2011;Wainwright & Reilly, 1994) are necessary to test these conclusions and evaluate a more indepth form-function relationship in Sciuromorpha.

ACK N OWLED G M ENTS
We are grateful to M. Zelditch for providing the sciurid phylogeny, S.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
JW and JAN conceptualized the study. JW collected the data, conducted the data analysis, and drafted the manuscript. Both authors interpreted the results and revised previous versions of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data and R files associated with this manuscript are archived on