Locomotory Adaptations in 3D Humerus Geometry of Xenarthra: Testing for Convergence

Three-dimensional (3D) models of fossil bones are increasingly available, thus opening a novel frontier in the study of organismal size and shape evolution. We provide an example of how photogrammetry can be combined with Geometric Morphometrics (GMM) techniques to study patterns of morphological convergence in the mammalian group of Xenarthra. Xenarthrans are currently represented by armadillos, sloths, and anteaters. However, this clade shows an incredibly diverse array of species and ecomorphotypes in the fossil record, including gigantic ground sloths and glyptodonts. Since the humerus is a weight-bearing bone in quadrupedal mammals and its morphology correlates with locomotor behavior, it provides an ideal bone to gain insight into adaptations of fossil species. A 3D sample of humerii belonging to extant and fossil Xenarthra allowed us to identify a signiﬁcant phylogenetic signal and a strong allometric component in the humerus shape. Although no rate shift in the evolution of the humerus shape was recorded for any clade, fossorial and arboreal species humerii did evolve at signiﬁcantly slower and faster paces, respectively, than the rest of the Xenarthran species. Signiﬁcant evidence for morphological convergence found among the fossorial species and between the two tree sloth genera explains these patterns. These results suggest that the highly specialized morphologies of digging taxa and tree sloths represent major deviations from the plesiomorphic Xenarthran body plan, evolved several times during the history of the group.


INTRODUCTION
Species morphology varies in size and shape. These two components can be strongly correlated to each other (Shingleton et al., 2007;Figueirido et al., 2011;Voje et al., 2014;Klingenberg, 2016) and somewhat limited by the existence of evolutionary constraints (Gould, 1989;Brakefield, 2006;Arnold, 2015;Meloro et al., 2015a). In this regard, the vertebrate skeleton has been intensively investigated, because the shape of its components is greatly influenced by body size and by the constraints impinging on specific adaptations linked to body support and other essential organismal functions (e.g., locomotion, feeding). The skeleton also allows sampling ancient diversity that in many clades can greatly overcome variation of extant taxa. In fact, the appreciation of fossil diversity provides strong support for the existence of size-induced shape changes (allometry) across different taxonomic scales and several components of the skeleton (Speed and Arbuckle, 2016).
The combined effects of such size-related shape changes and the evolutionary pressure originated by adaptation might generate patterns of morphological convergence in distantly related clades (Harmon et al., 2005;Mahler et al., 2010;Losos, 2011;Meloro et al., 2015b). Convergence is more likely to take place when adaptation is at the most extreme, and it can be identified in both extant and fossil taxa (i.e., the skulls of marsupial and placental carnivores including sabertooth morphologies, Wroe and Milne, 2007;Goswami et al., 2011).
Extant Xenarthra are currently limited to 31 species falling within the two clades Cingulata (armadillos) and Pilosa (sloths and anteaters; Simpson, 1980;Engelmann, 1985;Springer et al., 2003;Delsuc et al., 2004). Yet, in the past they showed much greater phenotypic and taxonomic diversity, encompassing some 700 species overall, including the gigantic late Pleistocene ground sloths and armadillos (Prothero, 2016). 1 Although recent advances in proteomics and genomics provide new insights into Xenarthra phylogenetic history, their position within Placental mammals is still a matter of controversy (Gibb et al., 2016).
Thanks to 3D modeling, it is currently possible to build precise replicas of fossil bones and investigate their size and shape variation with better accuracy than ever before. Applying such 3D modeling on Xenarthran limb elements is particularly welcome, given the great diversity of size and lifestyle the group experienced in its recent past .
We take advantage of the newly developed photogrammetry technique (Falkingham, 2012) to build a dataset of 51 Xenarthran humerus 3D models belonging to 29 species (16 extant plus 13 extinct). The advantage of photogrammetry is that it minimizes specimen handling (which is convenient for their fragile status) and allows a relatively quick data collection based on museum specimens (taking pictures for photogrammetry models might take between 5 and 10 min; Giacomini et al., 2019). On the other side, software post-processing time can still be quite long, although the development of professional software (e.g., 1 https://paleobiodb.org/#/ Agisoft) and novel open access sources are making the process increasingly quicker. 2 There has been a lot of research focusing on the adequacy and the accuracy of photogrammetry method for GMM analyses (see Giacomini et al., 2019, for a recent overview). Particularly for long bones, Fau et al. (2016) demonstrated that photogrammetry provides a good level of accuracy compared to other laser scanners (e.g., structured light Artec laser or Breukmann) on relatively medium-sized vertebrate long bones.
We explored Xenarthran humerus 3D morphology using GMM within a comparative framework. This technique is now well established with a long record of research applications also on fossil mammals (Adams et al., 2004(Adams et al., , 2013. The use of homologous points (landmarks) facilitates a comparison between species belonging to the same clade and additionally provides a powerful tool for separating and visualizing size and shape variations. GMM application to the study of Xenarthran functional morphology includes works from Monteiro and Abe (1999) focusing on the scapula or Milne et al. (2009Milne et al. ( , 2012 on humerus and femur shape. Interestingly, the implementation of the comparative methods into GMM datasets is a relatively more recent phenomenon. Comparative methods were first introduced by Felsenstein in his seminal paper on the independent contrasts (Felsenstein, 1985). This method explores the assumption that interspecific data are not independent, since species might share a different degree of common ancestry in any macroevolutionary dataset. The effect of shared inheritance can be assessed by estimating the "phylogenetic signal" in the data. When fossil species are concerned, two major limitations occur: fossil phylogenies are morphology based, so that testing morphological hypotheses using these trees might generate a circular argument (in the majority of cases fossil dataset might exhibit higher phylogenetic signal); fossil phylogenies are generally incomplete with taxonomic confirmation generally scattered between different publications. In spite of this, early attempts demonstrated that fossil phylogenies can be incorporated to test macroevolutionary hypotheses, and their inclusion provides stronger statistical power (for early examples see: Finarelli and Flynn, 2006;Meloro et al., 2008). We provided on several occasions examples on how comparative methods can be implemented in macroevolutionary studies incorporating fossils and GMM Raia et al., 2010;Meloro and Slater, 2012;Piras et al., 2012). More recently, the development of new R packages (including geomorph Adams et al., 2019;and RRphylo Raia et al., 2020) allows to detect evolutionary rates with a high degree of accuracy (in spite of phylogenetic fossil uncertainty, e.g., Smaers et al., 2016;Castiglione et al., 2018).
We take advantage of the most recently published phylogenies for fossil sloths and armadillos in order to test hypotheses about the influence of size and locomotor behavior on humerus shape and rate of its morphological change in Xenarthra. Since behavioral and morphological convergence has been proposed for extant sloths, as well as fossil Megatherium and extant armadillos (Nyakatura, 2012;Billet et al., 2015), we explicitly tested convergence in humerus morphology using a novel approach that can well be applied to multivariate shape data and phylogenies inclusive of extant and fossil species .

Studied Specimens
We built 3D models of Xenarthra humerii belonging to 29 species housed at the following museum institutions: the Natural History Museum (NHLM, London) and the Muséum National d'Histoire Naturelle (MNHN, Paris, see Supplementary Table  S1 for details).
For each specimen, we took about 200 photos based on dorsal, ventral, and lateral views using a dark background and a standard digital SLR (Nikon D5300, lens 18-140 mm). Most of the pictures were taken with the 55 mm lens setting on a fixed focus (see Giacomini et al., 2019). The software Agisoft Metashape was then employed to build 3D models, while MeshLab software (v1.3.3, 2014) was used for scaling them. The scaling was based on one single measurement (generally the maximum bone length), as recommended by Falkingham (2012). Sensitivity analyses testing relative bone proportion were performed on selected specimens of similar-sized mammalian humerii, generally providing a difference smaller than 5% when comparing measurements taken with digital caliper and those with the MeshLab software (v1.3.3, 2014).

Landmarking
The software Landmark (v. 3.0) was employed to identify on each specimen 28 landmarks (Figure 1). The landmarks were designed to cover main anatomical regions of the Xenarthran humerus including proximal and distal epiphyses as possible. The Landmark descriptions are in Table 1 following anatomical nomenclature (Figure 1). Most of the landmarks were type 2, since humerus epiphyses do not allow to identify type 1 points. However, these were previously evaluated in Milne et al. (2009), which we followed as a baseline for our configuration (see Table 1 and Figure 1). Landmarking was performed twice on a subsample of 20 specimens to detect the level of error in size and shape using Procrustes ANOVA, which in all cases turned out to be nonsignificant explaining less than 5% of inter-individual variation.

GMM and Comparative Methods
The 3D landmarked coordinates of the scaled models were subject to GPA (Generalized Procrustes Analysis, Rohlf and Slice, 1990). This technique removes the non-shape information related to size, position, and orientation of the specimens. GPA returns a new set of coordinates subsequently subjected to Principal Component Analysis (PCA, Rohlf and Slice, 1990). PCA decomposes the shape variation into orthogonal axes of maximum variation named PCs (Principal Components). PC vectors are used as variables in subsequent analyses. Species mean shapes were calculated before performing analyses so that our shape data represented inter-specific variation only. Humerus size was quantified using the natural logarithm (Ln) of the centroid size (=the square root of the sum of the squared distances from each landmark to the centroid of each configuration).
We assembled two different phylogenetic trees. The first was based on molecular data following Presslee et al. . We decided to use two different phylogenies because molecular vs. morphological trees may be conflicting (Cohen, 2018). For the phylogenetic position of species for both trees and references see Supplementary Table S2. The phylogenetic trees were calibrated by using the scaleTree function in the RRphylo package , species last appearance and internal node ages used for calibration are in Supplementary Table S2).
By using plotGMPhyloMorphoSpace function in the geomorph package (Adams et al., 2019), the trees were mapped into PCA in order to generate a phylomorphospace. Phylogenetic signals were quantified by using the K (for size) and the K multiv (for shape) statistic (Adams, 2014).
Phylogenetic Generalized Least Squares (PGLS) regression was employed assuming Brownian motion as the mode of evolution to test for macroevolutionary allometry and differences among locomotory categories. We partitioned species into discrete stance categories as proposed by  based on the main lifestyle of extant taxa. Tree sloths were considered as fully "arboreal." The anteaters, capable of both climbing and having a unique digging style, were classified as "intermediate" (Hildebrand, 1985;Kley and Kearney, 2007), whereas armadillos were divided into fully "terrestrial" and "fossorial." For fossil species, we carried out an extensive revision of the literature (see Supplementary Table S3). Extant members of Cingulata are specialized diggers as inferred by their limb morphology (Vizcaíno and Milne, 2002;Milne et al., 2009;Marshall, 2018;Mielke et al., 2018). Other specialized diggers can be found also among Pilosa. Glossotherium robustum was demonstrated to be a specialized digger (Bargo et al., 2000;Vizcaíno et al., 2001;de Oliveira and Santos, 2018). On the contrary, other ground sloths were best adapted to a terrestrial lifestyle (Bargo et al., 2000;Vizcaíno et al., 2001;de Oliveira and Santos, 2018). Extant sloths (Bradypus and Choloepus) are known as tree sloths for their strictly arboreal lifestyle (Montgomery, 1985;Chiarello, 2008;Toledo et al., 2012), as well as the anteaters (White, 2010), but the latter are capable of above branch locomotion (Nyakatura, 2011). Locomotor categories were equally employed to test for convergence.
We assessed the rate of humerus size and shape evolution by using the RRphylo function (Castiglione et al., 2018). RRphylo returns a vector of evolutionary rates for all branches in the tree and a vector (or a matrix if the phenotype is multivariate, i.e., with the shape data) of ancestral states estimated for each node. We applied RRphylo on both size (log-transformed) and shape (as PC scores and by using size as covariate).
To search for possible shifts in the evolutionary rates between clades or locomotory categories we used the function search.shift (Castiglione et al., 2018). First, we applied search.shift on shape and size under the "clade" condition. In this case, the function compares the average absolute evolutionary rate of a specific clade with the rest of the tree. The significance level is assessed by randomizations. Then, we applied search.shift on shape and size under the "sparse" condition to test for differences in evolutionary rates among locomotor categories. Under this condition, the function tests if species having the same state evolve differently from the others.
Eventually, we applied the new function search.conv  to test for morphological convergence. This function tests morphological convergence by assessing the angle between phenotypic vectors (vectors of PC scores) between species and comparing this angle to a random expectation. Given two phenotypic vectors (here PC scores), the cosine of angle θ between them represents the correlation coefficient (Zelditch et al., 2012). Under the Brownian Motion, θ is expected to be proportionally related to phylogenetic distance. Yet, convergence violates this assumption. The new method search.conv, calculates the θ angles between entire clades ("automatic") or species under the same state ("state"). It tests whether the mean θ between the species evolving under a specific state or belonging to a single clade is smaller than expected by chance, and whether θ divided by the mean phylogenetic distance among convergent tips is smaller than expected. We applied search.conv to test for morphological convergence between fully arboreal sloths and the intermediate anteaters, and species evolving under fossorial lifestyles. Then, we tested convergence between the two tree sloths: Bradypus and Choloepus.
To test the robustness of our results to phylogenetic uncertainty and sampling effects, we applied the new implemented function overfitRR Melchionna et al., 2020) in the RRphylo package . Under this function, the original tree is trimmed of a predetermined number of species (here 5%) and the position of tips is changed randomly by up to two nodes away from its original placement. For instance, a simple phylogenetic tree [(A, B), C], would be changed in [(C, B), A] or [(A, C), B] tree topology. In addition, the function changes the node age, randomly, within a range between the age of the focal node immediate ancestor node and the age of its older daughter node. We ran overfitRR with 100 iterations. At each iteration, the function swapped 5% of the tree size of tips and changed in age 5% of the tree size of tree nodes, then used a new phylogenetic tree to perform search.conv to either confirm or reject any instance of significant morphological convergence. All the analyses were performed both with the molecular and the morphological phylogenies. Table S4), of which the first 13 explained up to 95.80% of the shape variation. More in detail, PC1 and PC2 explained 53.58% and 15.72%, respectively, and they show a degree of separation in humerus shape between Cingulata and Folivora, while Vermilingua clade overlap with the latter (Figure 2).

GMM returned 28 PCs (Supplementary
Deformations described by PC1 are related to the relative elongation of the humerus. Arboreal species occupy positive scores of PC1 and are characterized by a relatively longer and slender humerus morphology. Both proximal and distal epiphyses are reduced. This configuration is typical of the suspensory tree sloths (genera Bradypus and Choloepus). On the contrary, PC1 negative scores describe short and stocky humerus shape. The head and both the greater and lesser tuberosities are more expanded, as well as the trochlea, capitulum, and epicondyles (Figure 2).
Changes along PC2 axis described differences among Xenarthra clades. Folivora and Vermilingua occupy positive PC2 scores. In these clades, the humerus head is slender, the curvature between the deltoid crest and the distal part is wider and developed backward. On the negative scores of PC2, Cingulata shows distinct humerus head components with the distal part longer and slender (Figure 2).

Comparative Methods
When applying the molecular tree, the phylogenetic signal for both shape and size was statistically significant. The observed K (size) and K multiv (shape) are 0.997 and 0.453 (p value = 0.001), respectively. Very similar results were obtained when the morphological tree was used (K size = 0.939; p value = 0.001, K multiv shape = 0.497, p value = 0.001).
Differences in shape were explained by differences in size, as resulted by using Procrustes ANOVA without (Table 2A) and with the addition of the phylogenetic effect of the morphological tree ( Table 2C). The allometric effect remained significant in Folivora and Cingulata when the phylogenetic relationship was accounted for (Table 3 and Figure 3). Procrustes ANOVA showed a significant impact of locomotor categories on shape, either without considering phylogenetic effects or by using the molecular tree (Table 2). However, locomotion had no impact on humerus size variation ( Table 2).
By applying search.shift to the humerus size, we found a significant increase in evolutionary rates in the clade, including Mylodontidae and Choloepus (average rate difference = 0.080, p = 0.001; Figure 4A), and a significant decrease in the clade, including Panoctus, Doedicurus, Neosclerocalyptus, and Hoplophorus (average rate difference = −0.025, p = 0.006; Figure 4A). This same negative shift applied when the morphological tree was considered (average rate difference = −0.020, p < 0.001; Figure 4B), whereas a positive rate shift was found for the clade including Euphractinae and   Chlamyphorus (average rate difference = 0.047, p < 0.001; Figure 4B) by using this tree. We did not find significant shifts in the rate of humerus shape evolution with either tree. By testing for rate shifts in humerus size per locomotory state, we found a positive and significant difference in rates pertaining to strictly arboreal species (tree sloths) by using molecular trees (rate difference = 0.123, p value = 0.001, Table 4). Similarly, tree sloths showed significantly higher rates of shape evolution as compared to the rest of the taxa (rate difference = 11.259, p value <0.001, Table 4). With this same molecular tree, fossorial species had slower rates of humerus shape evolution as compared to the species falling in different locomotory states (rate difference = −2.040, p = 0.002, Table 4). When applying the morphological tree, these differences were less apparent although fossorial species still showed slower, and arboreal species higher rates as compared to the rest of the taxa, either by analyzing humerus size or shape (Table 4).
By using search.conv, we tested the convergence in humerus shape between species living on tree branches although ascribed to different categories (intermediate and arboreal). search.conv returned a non-significant angle of 98.253 • (p value = 0.867) between the species under the two states, the same applied when the angle was tested per unit time (p value = 0.099; Table 5A). A strong evidence for convergence appears within the strictly arboreal species, that is the Bradypus and Choloepus clades (θ real = 0.235, pθ real = 0.010; (θ real+θace )/time = 0.383, p( θreal+θace )/time = 0.001; Table 5B and Figures 5B,D). Morphological convergence was additionally found for species evolving in the "fossorial" category (mean angle among fossorial species = 68.148 • ; p angle state = 0.002, Table 5A and Figure 5C). This notion is confirmed when the time distance between species is accounted for (θ time = 0.786, pθ time = 0.001, Table 5A).
overfitRR returned 0% of significant simulations when convergence between the arboreal and intermediate states is tested for both phylogenies. This is in agreement with search.conv results.
For the issue of convergence between Bradypus and Choloepus overfitRR returned 100% (for both θ real and θ real+θace )/time) instances of significance by using molecular tree. When the same analysis is performed by using morphological tree, overfitRR returned 71% for θ real and 74% for θ real+θace )/time instance of significant convergence.
For the issue of convergence between fossorial species, overfitRR returned 100% (p.θ state ) and 99% (p.θ time ) of significant p values when molecular tree is used. The same figures for the morphological tree were 100% (p.θ state ) and 98% (p.θ time ).

DISCUSSION
The taxonomic and phenotypic diversity of extant Xenarthrans represents only a small fraction of their past variations. Perhaps unsurprisingly, giant extinct glyptodonts (Doedicurus, Glyptodon, and their allies) experienced a shift toward a slow rate of size evolution, in keeping with their uniformly large body size. This stands true irrespective of whether molecular or morphological trees are used. In contrast, a significant positive shift in the rate of humerus size evolution appears in Euphractinae plus Chlamyrophus by using the morphological tree only, with Eutatus placed outside the clade of Glyptodontinae plus Euphractinae ( Figure 4B). This probably related with a decrease in size of Euphractinae plus Chlamyrophus in contrast with the plesiomorphic condition of Eutatus, whose estimated body size is supposed to range between 36.8 and 71.7 kg (Vizcaíno and Bargo, 2003). On the contrary, extant Euphractine plus Chlamyrophus are only 2 kg in body size on average.
We found no significant instance of shape evolutionary shifts pertaining to the Xenarthran trees, meaning that major shape differences were channelled through phylogeny as further supported by GMM results. Phylomorphospace showed a neat separation between the two clades Cingulata and Pilosa along PC2 (Figure 2).
However, the humerus changed significantly relative to locomotory styles. The rate of shape evolution in fossorial   Xenarthrans is significantly smaller than in other species. Such fossorial habitus characterizes most Cingulata and is possibly plesiomorphic to the group (Vizcaíno and Milne, 2002;Milne et al., 2009;Marshall, 2018), although it is surprisingly present among giant ground sloths like Glossotherium (Bargo et al., 2000;Vizcaíno et al., 2001;de Oliveira and Santos, 2018). Conversely, arboreal species evolved at faster rates than in any other Xenarthra. Tree sloths were noted for their highly derived and convergent morphologies, canalizing change from a plesiomorphic, fossorial condition (Nyakatura, 2011;Nyakatura and Fischer, 2011), although it must be emphasized that scansorial species (which were not included in the analyses here) were present since the Miocene (Pujos et al., 2012;Toledo et al., 2017). Overall, our results suggest that the acquisition of arboreal and fossorial lifestyles are loaded with great functional demands, leading to constrained, little variable morphologies (fossorial) or to significant changes in shape toward a particular configuration (from the ancestral conditions) to match a demanding lifestyle (arboreal sloths). These challenging adaptations present ideal cases to test convergence. We provided this test for both fossorial species and tree sloths, separately, and confirmed they do represent significant instances of morphological convergence (Figure 5). The same applies if phylogenetic uncertainty is accounted for. The results further support the observation that phylogenetic signal in humerus shape is as much high as for humerus size (Vizcaíno et al., 1999;Vizcaíno and Bargo, 2003), but also suggest that tree-living is the room for morphological change among Xenarthran taxa.
Limb proportions in fossil Eutatus were similar to extant Euphractine (Vizcaíno and Bargo, 2003). These species present limbs specialized in building borrows rather than in gathering food (Vizcaíno and Bargo, 2003). Thus, the high size evolutionary rates can be related with more robust humerii specialized for digging. It was demonstrated that Miocene glyptodonts weighed about 100 kg, but Pleistocene species may have weighed up to 1 ton (Vizcaíno et al., 2010;Milne et al., 2012). Our results indicate that the lowest evolutionary rate in humerus size occurs in Panoctus, Doedicurus, Neosclerocalyptus, and Hoplophorus suggesting humerus size in these taxa was linked to their large but uniform body size.
Xenarthra originated in South America about 62.5 million years ago (Presslee et al., 2019). During the Great Biotic Interchange in late Pliocene (Marshall et al., 1979;Webb, 2006) FIGURE 5 | (A,C) PC1 vs. PC2 scatter plot. Convex hulls are colored according to different states. Circular plots represent the mean angle between states, the gray area is the random angle range. (B,D) Traitgram plot with colored branched representing convergent clade. PC1 vs. PC2 scatter plot are produced, convergent clades are color coded. In both, traitgram and PC1 vs. PC2 scatter plot asterisks represent the ancestral phenotypes of the individual clades.
several species migrated to the North (Bocherens et al., 2017). Studies of fossil species have demonstrated that none of the known fossil sloths had arboreal lifestyle (White, 2010;Nyakatura, 2011). The last common ancestor of sloths probably was terrestrial or semi-arboreal (White, 2010;Nyakatura, 2012). Indeed, fossil sloths appear morphologically closer to extant Vermilingua (i.e., Tamandua and Myrmecophaga, herein classified as "intermediate") rather than to extant tree sloths . We did not find instance of convergence among intermediate and arboreal species (Figure 5). This supports the idea that modern tree sloths acquired the suspensory habitus secondarily, which explains their higher shape evolutionary rates as compared to the humerii of species ascribed to different locomotor categories. Similarly, the long branch separating the extant tree sloth genera (Presslee et al., 2019) are suggestive of secondary adaptation. Extant sloths present a forelimb-dominated locomotion. Bradypus moves up to 10 m only using his forelimbs, and Choloepus hind limbs lost their primarily propulsive elements (Mendel, 1985;Nyakatura et al., 2010).
Similarly, since digging kinematics is one of the most demanding behaviors the mammalian skeleton could be designed for Sansalone et al. (2019), the pervasive call for convergence to a functionally optimal design linked to digging was expected (Sansalone et al., 2020).
One obvious caveat we urge to consider is that, although nearly one-half of the species we considered in our tree are extinct, the history of Xenarthra cautions against giving too much faith to phylogenetic analyses using a tree devoid or otherwise scarce in terms of fossil species representation. The inclusion of fossil phenotypes is, and must carefully be, considered in trait evolution inference, especially when major patterns such as morphological convergence are sought after.

DATA AVAILABILITY STATEMENT
The datasets for this study are included in the Supplementary Material.

AUTHOR CONTRIBUTIONS
CS, PR, and CM conceived the project. CS and CM prepared samples and conducted the data analyses. All the authors discussed the results and wrote the manuscript.

FUNDING
CS was supported for this study by LJMU Ph.D. studentship fund.