Analyzing Xenarthran Bone Microstructure

For all regressions and flexible discriminant analyses (FDAs), we chose to use both phylogenetic and nonphylogenetic methods because ecology and phylogeny are closely aligned in extant xenarthrans, likely making it difficult to cleanly isolate phylogenetic and ecological signals (Fig. 1) (Zack et al. 2021 ). Using both phylogenetic and nonphylogenetic methods will allow us to directly compare an analysis that assumes that phylogeny has no effect to an analysis that assumes the similarity of organisms is due to their relatedness rather than to convergence (Revell 2010 ). The phylogenetic regressions we used do not allow for more than one specimen from each species, so we had to use the average of each. Because of this, we also calculated regressions using the species means for each vertebral position. If there is no impact of using the species means, then we would expect the species mean regressions to be exactly the same as the nonphylogenetic specimen-based regressions. Comparing these three methods will allow us to better understand how the inclusion of phylogenetic information is affecting the slope estimation in these regressions.
To statistically analyze the impact of body size on each bone microstructure metric, we used GLS regressions on log-transformed metrics for the entire dataset (Table 2, Fig. 5). We additionally calculated the linear regressions for species averages (Table S3, Fig. S3). BV.TV and GC are ratios and are therefore unitless with an isometric slope of 0. Tb.N is measured in numbers/mm, so the isometric slope is −1. Tb.Th is a length value with an isometric slope of 1; CSA is a measure of area with an isometric slope of 2. Isometric slope for Conn.D is −3 because it is measured in numbers/mm³ (see Mielke et al. 2018 for further explanation of isometry for these metrics; Plasse et al. 2019 (link)). We also calculated confidence intervals (CIs) for all regressions using the confint function in R (Team RC 2022 ). A slope was considered allometric if the isometric slope fell outside the CI. We calculated the regressions for each of the three clades examined (Cingulata [armadillos], Vermillingua [anteaters], and Folivora [sloths]) and for each of the ecologies (arboreal, hook-and-pull digging, and scratch digging). We then compared the slopes using the standardized major axis estimation and testing function in the smatr package in R (Tables 3 and 4) (Warton et al. 2012 ; Team RC 2022 ). We also used PGLS regressions to determine the impact of phylogenetic covariance on the scaling of TBA metrics (Table 2, S2). We pruned the time-scaled tree of Gibb et al. (2015) (link) to include only the species in our dataset. We prepared the data by calculating the species means of each metric by vertebral position. We used the GLS function and corBrownian in the R package ape to calculate the PGLS of each metric (Table 2, S2) (Paradis and Schliep 2019 (link)). We calculated Blomberg's K using the R package phytools to further analyze the phylogenetic signal of each variable (Revell 2012 ). Because we analyzed each vertebral position separately, we used the average of the PGLS regressions of each metric along the vertebral column to compare to the individual linear regressions.
To further quantify the impact of size, ecology, and phylogeny on TBA, we used both pFDA and FDA. We determined ecology groups using the primary locomotor ecology of each genus (Rood 1970 ; Ramsey 1978 ; Navarrette and Ortega 2010 ; Hayssen 2011 ; Hayssen et al. 2012 ; Gaudin et al. 2018 ; Attias et al. 2020 ), and we determined size class based on the groups from GLS regressions (Table 2). Using the same data as the PGLS, we performed the pFDA using code from Motani and Schmitz 2011 (link) and Smith et al. 2018 ). This package can only use up to three metrics, so we could not undertake a fully multivariate analysis of our dataset. Therefore, we used three subsets of metrics to complete each analysis: the most size-correlated metrics (Tb.Th, CSA, and Conn.D), the least size-correlated metrics (BV.TV, GC, and DA), and the most phylogenetically-correlated metrics (DA, Tb.Th, and CSA). We chose to use the least size-correlated metrics in our analyses because this model most accurately resolved ecology (Table S4, Fig. S4). For the FDA, we used the mda and nnet R packages to complete the analysis and visualize the results (Venables and Ripley 2002 ; Leisch and Hornick 2022 ).