Integrative Approach Uncovers New Patterns of Ecomorphological Convergence in Slow Arboreal Xenarthrans

Identifying ecomorphological convergence examples is a central focus in evolutionary biology. In xenarthrans, slow arboreality independently arose at least three times, in the two genera of ‘tree sloths’, Bradypus and Choloepus, and the silky anteater, Cyclopes. This specialized locomotor ecology is expectedly reflected by distinctive morpho-functional convergences. Cyclopes, although sharing several ecological features with ‘tree sloths’, do not fully mirror the latter in their outstandingly similar suspensory slow arboreal locomotion. We hypothesized that the morphology of Cyclopes is closer to ‘tree sloths’ than to anteaters, but yet distinct, entailing that slow arboreal xenarthrans evolved through ‘incomplete’ convergence. In a multivariate trait space, slow arboreal xenarthrans are hence expected to depart from their sister taxa evolving toward the same area, but not showing extensive phenotypical overlap, due to the distinct position of Cyclopes. Conversely, a pattern of ‘complete’ convergence (i.e., widely overlapping morphologies) is hypothesized for ‘tree sloths’. Through phylogenetic comparative methods, we quantified humeral and femoral convergence in slow arboreal xenarthrans, including a sample of extant and extinct non-slow arboreal xenarthrans. Through 3D geometric morphometrics, cross-sectional properties (CSP) and trabecular architecture, we integratively quantified external shape, diaphyseal anatomy and internal epiphyseal structure. Several traits converged in slow arboreal xenarthrans, especially those pertaining to CSP. Phylomorphospaces and quantitative convergence analyses substantiated the expected patterns of ‘incomplete’ and ‘complete’ convergence for slow arboreal xenarthrans and ‘tree sloths’, respectively. This work, highlighting previously unidentified convergence patterns, emphasizes the value of an integrative multi-pronged quantitative approach to cope with complex mechanisms underlying ecomorphological convergence.


Introduction
Convergent evolution is defined as the independent acquisition of similar features in phylogenetically distant lineages (Stayton 2015a). Some astonishing examples of adaptation are possibly explained by the functional convergence of morphologies in distant clades occupying similar ecological niches, a pattern known as ecomorphological convergence (Wainwright and Reilly 1994;Schluter 2000;Muschick et al. 2012). In this regard, mammals, characterized by wide ecological and morphological variety, have represented one of the most extensively investigated groups, yielding the identification of outstanding instances of ecomorphological convergence (McGhee 2011). Textbook examples of functional convergences in mammals are the short and heavily-built limbs adapted for burrowing in fossorial species (Hildebrand 1985), the evolution of the patagium that enables lift in gliders (McGhee 2011;Futuyma 2013;Reece et al. 2014), and or a tongue morphology that allows anteating taxa to extract prey from nests (Griffiths 1968;Redford 1987;Reiss 2001). A deep understanding of mechanisms behind convergence is crucial for explaining mammal morphological diversification. Yet, several aspects of convergence remain to be clarified.
Despite implying phenotypic predictability from ecology (Losos et al. 1998;Collar et al. 2014), convergence may not manifest with the expected patterns and magnitude. Analyses of ecomorphological convergence yielded only partial confirmation for functional traits expected to be shared by mammals occupying the same ecological niche (e.g., Meloro et al. 2015;Grossnickle et al. 2020;Alfieri et al. 2021). Several factors may be behind these outcomes, such as historical contingency/constraints (Harvey and Pagel 1991;Losos and Miles 1994;Zelditch et al. 2017) or stochastic evolution (Stayton 2008(Stayton , 2015b. Another such factor could relate to the concept known as 'many-to-one mapping of form onto function', which is the idea that distinct morphological configurations might achieve the same function with equal efficiency (Hulsey and Wainwright 2002;Alfaro et al. 2005;Wainwright et al. 2005;Collar et al. 2014;Zelditch et al. 2017).
A central principle of convergent evolution states that greater morphological similarity is expected in clades more strongly convergent in lifestyle (i.e., subjected to stronger selective pressure) (Conway Morris 2010). A corollary is that clades that only partially converge in ecology are hypothesized to morphologically converge to a lesser extent. These intermediate degrees of convergence, characterized by different patterns of morphological convergence, may be recognized in multivariate trait spaces ('morphospaces', hereafter). These different patterns do not effectively correspond to different natural processes; they only represent abstractions introduced by evolutionary biologists studying morphospaces to consider nuanced degrees of convergence (e.g., Stayton 2006;Grossnickle et al. 2020). Cases of particularly strong morphological similarity (i.e., resulting in extensive overlap in the morphospace; e.g., Melville et al. 2006;Cooper and Westneat 2009;Meachen-Samuels 2012;Arbour and Zanno 2020) are referred to as 'complete' convergence (Losos 2011;Stayton 2006: fig. 3c). This type of convergence is evident when all convergent taxa occupy a distinctly smaller region of the morphospace, departing from that of close relatives. Another scenario may occur if putatively convergent taxa occupy a smaller sub-region of the morphospace compared to their ancestral relatives but do not all overlap. In other words, convergent taxa overall acquire some similarities, but the non-overlapping ones retain a unique morphotype. Such a pattern is referred to as 'incomplete' convergence (Herrel et al. 2004;Stayton 2006Stayton : fig. 3b, 2015aGrossnickle et al. 2020: fig. 3c).
Assessing the degree of morphospace overlap to discriminate complete and incomplete convergence involves an element of subjectivity. Indeed, complete convergence can be considered as a merely theoretical pattern. Identical phenotypes are extremely rare in nature, and a more detailed inspection of clades completely overlapping on morphospaces reveals more differences on a finer scale. Statistical tools allow coping with this issue. To decrease subjectivity in the identification of stronger similarity and degree of overlap, one can use measures of convergence strength (the measure of Stayton 2015a is used in this study; see also Castiglione et al. 2018) and/or morphological disparity (e.g., Stayton 2006;Arbuckle et al. 2014;McLean et al. 2018;Grossnickle et al. 2020). Moreover, evolutionary trajectories on phylomorphospaces may suggest different convergence patterns (Stayton 2006(Stayton , 2015a. Incomplete convergence can be recognized by evaluating how convergent taxa deviate from sister clades in the phylomorphospace. Trajectories of convergent taxa that are directionally similar but do not overlap may indicate incomplete convergence. The lower degree of overlapping may also be caused by divergent trajectories for some of the putatively convergent taxa (e.g., Grossnickle et al. 2020). Due to theoretical problems when trying to objectively distinguish between complete and incomplete convergence, the combination of the aforementioned lines of evidence makes the identification of these two models of convergence less subjective.
Analysing several levels of morphological complexity, it is possible to identify additional patterns of morphological convergence. The suite of investigated functional features may crucially bias the detection of convergence, since this evolutionary process potentially follows mosaic patterns. Indeed, only some traits might converge even within the same anatomical element (Spear and Williams 2020). Moreover, different scales of anatomical detail of the same element, relying on different structures and mechanisms, could evolve through different convergence patterns (Watanabe et al. 2020). As a result, analyses based on one or few traits and/or focusing on a specific scale of investigation may fail to or only partially reveal convergence.
To study ecological and morphological covariation in mammals, the superorder Xenarthra provides a suitable context. Xenarthrans are a clade of placentals nowadays represented by 'tree sloths', anteaters, and armadillos (members of Folivora, Vermilingua, and Cingulata, respectively, with the first two grouped in Pilosa; Fig. 1). They provide compelling examples of functional adaptations, with their morphology showing high evolutionary plasticity (Billet et al. 2012;Amson et al. 2017;Amson and Nyakatura 2018a, b). Major locomotor adaptations of extant xenarthrans reflect the group's taxonomic diversity. Armadillos are generally characterized by a fully terrestrial locomotor ecology dominated by digging habits (Vizcaíno and Milne 2002;Attias et al. 2016;Amson et al. 2017). Most anteaters are at least semi-arboreal and the three vermilinguan genera represent a cline from fully terrestrial (Myrmecophaga) to fully arboreal (Cyclopes) locomotion, with Tamandua spanning the two ends of the behavioral spectrum (Young et al. 2003;Orr 2005;Amson et al. 2017). 'Tree sloths' are known for their distinctive suspensory slow arboreal locomotion (Nyakatura 2012).
The silky anteater Cyclopes didactlyus (the only species of the genus, referred to as Cyclopes hereafter) has generally been described as slow arboreal (van Tyne 1929;Hayssen et al. 2012;Nagy and Montgomery 2012;Granatosky et al. 2014) but not recognized as mirroring 'tree sloths' in behavior. It may be due to its elusive and nocturnal habits (Hayssen et al. 2012;Nagy and Montgomery 2012), which have prevented detailed observations. Moreover, previous analyses on Cyclopes have mainly been focused on locomotor habits shared with other anteaters (e.g., Amson et al. 2017). Nevertheless, several ecological/biomechanical adaptations are present in both 'tree sloths' and Cyclopes, which can justify ascribing them to a slow arboreal ecology (following Alfieri et al. 2021). Bradypus, Choloepus, and Cyclopes are overall characterized by cautious movements on trees, dominance of rest and quiescence in their daily activity balance, and extremely low basal metabolic rates (van Tyne 1929;Hayssen et al. 2012;Nagy and Montgomery 2012;Nyakatura 2012;Pauli et al. 2016;Vendl et al. 2016). From a strictly locomotor perspective, Cyclopes is not as adapted to inverted quadrupedalism (i.e., suspensory locomotion) as 'tree sloths' (Nyakatura 2012), and it can be more accurately referred to as a vertical climber (van Tyne 1929;Hayssen et al. 2012;Nagy and Montgomery 2012). However, the locomotion of the silky anteater may also be categorized as antipronograde (sensu Granatosky et al. 2014), in that forelimbs and/or hindlimbs are mainly loaded in tension, a feature shared by suspensory and vertical climbers (Preuschoft 2002;Nyakatura and Andrada 2013;Schmitt 2017, 2019;Hanna et al. 2017). Several traits of this species have previously been identified as possible adaptations to fully arboreal vertical climbing (White 1993;Fonseca et al. 1996;Nowak 1999), some of which are shared with Choloepus and Bradypus (White 1993;Lewton and Dingwall 2016;Alfieri et al. 2021;Toledo et al. 2021). Hereafter, we will consider 'tree sloths' and the silky anteater as slow arboreal species, meaning that this type of locomotion evolved three times independently in xenarthrans (Fig. 1). Our aim is to identify morphological convergences driven by slow arboreal locomotor ecology in xenarthrans. Hence, in addition to slow arboreal xenarthrans, it is crucial to examine taxa adapted to others types of locomotion, collectively described as non-slow arboreal. The latter includes, for instance, some medium-sized extinct sloths from the early Miocene of Patagonia (Santa Cruz Formation;White 1993;Perkins et al. 2012; Fig. 1). We base this locomotor ecology assignment on previous ecological reconstructions (as detailed in Materials and Methods).
Here, we study convergence in the humerus and the femur, which have a well-known capacity to respond to biomechanical loadings and functionally adapt to locomotion style (Pearson and Lieberman 2004;Ruff et al. 2006;Kivell 2016), resulting in a close relationship with locomotor ecology (e.g., Patel et al. 2013;Botton-Divet et al. 2016;Amson et al. 2017;Mielke et al. 2018b;Fabre et al. 2019;Parsi-Pour and Kilbourne 2020). Both bones can be investigated at different scales of organisation (Francillon-Vieillot et al. 1990). We analyzed external shape, diaphyseal anatomy, and epiphyseal internal structure, which were quantified with 3D geometric morphometrics (3D GM), cross-sectional properties (CSP), and trabecular bone architecture, respectively (Fig. 2), due to their strong correlation with locomotion and ecology (e.g., Ryan and Ketcham 2002;Harmon 2007;Patel et al. 2013;Botton-Divet et al. 2016;Amson et al. 2017;Scheidt et al. 2019). These three aspects are also assumed to differently respond to the biomechanical milieu. External shape is argued to be more phylogenetically and anatomically constrained (Kivell 2016). Weak covariation has been reported between diaphyseal and epiphyseal variables Saers et al. 2016), since the former preferentially responds to bending/torsion (Carter and Beaupré 2007), while in the latter, axial loadings predominate (Biewener et al. 1996;Pontzer et al. 2006;Barak et al. 2011). Despite the fact that data provided by these three aspects could be complementary, rarely have they been combined in an integrated manner. To our knowledge, it has been done by Shaw and Ryan (2012), Sylvester and Terhune (2017), Amson and Nyakatura (2018a) and Saers et al. (2021), including only two of the aforementioned levels. Here, we combine the three types of variables to gain insights into functional morphological convergence driven by slow arboreal locomotor ecology. Several traits are shared by the humeri and femora of slow arboreal xenarthrans (e.g., White 1993;Straehl et al. 2013;Toledo et al. 2013Toledo et al. , 2015Amson and Nyakatura 2018b;de Oliveira and Santos 2018;Marshall et al. 2021) including features related to the anatomical scales and/or techniques employed in this work. Indeed, humeral and/or femoral features shared by slow arboreal xenarthrans have been highlighted analysing the external 3D shape (Milne and O'Higgins 2012;Mielke et al. 2018a), the diaphyseal structure (Patel et al. 2013;Marchi et al. 2016;Amson and Nyakatura 2018a;Montañez-Rivera et al. 2018;Alfieri et al. 2021) Fig. 2 Humerus of Choloepus didactylus NMW B5971 (left) and femur of Bradypus tridactylus ZMB Mam 7614 (right). External shape was quantified with 3D GM, using anatomical landmarks (red) and curve (blue) and surface (green) sliding semi-landmarks. The internal structure of the diaphysis was quantified both averaging CSP on the diaphysis and taking them from the 50% level (at mid-length). Orange rectangles show the cross-sections corresponding to the 30% and 70% levels, while green rectangles highlight the 50% crosssections. The epiphyseal internal structure was analyzed by extracting a spherical volume of interest (VOI) from the humeral head, the humeral capitulum, and the medial and lateral condyles of the distal femur. A hemispherical VOI was used to sample femoral head cancellous bone. Figures are not to scale. Cross-sections and VOIs are not in the same orientation as the whole bones and the epiphyseal internal architecture (Amson et al. 2017;Amson and Nyakatura 2018a). However, among these the studies that deal directly with convergence, they only do so qualitatively it (to our knowledge, only Serio et al. 2020 and Spear and Williams 2020 employed a quantitative approach to humerus external shape in xenarthrans).
In this work, we first partitioned humeral and femoral anatomy into single traits and identified those that correlate with the slow arboreal locomotion. Using this subset of features that are expected to reflect locomotor ecology, we tested for convergence on multivariate sub-datasets, first pooling traits by anatomical levels then collectively for each bone. Since slow arboreality is a highly specialized locomotion, we expect to identify a representative set of functional features discriminating slow arboreal xenarthrans. Furthermore, within slow arboreal xenarthrans, we expect Bradypus and Choloepus to exhibit the highest degree of convergence due to their extraordinarily similar adaptations to suspensory locomotion (Nyakatura 2012). We hereafter refer to the two genera of 'tree sloths' as suspensory slow arboreal xenarthrans, representing a subset within slow arboreal xenarthrans, showing a higher degree of similarity to one another than to other xenarthrans. Cyclopes is expected to share more similarities with 'tree sloths' than with other anteaters because of its slow arboreal locomotion dominated by vertical climbing (van Tyne 1929;Hayssen et al. 2012;Nagy and Montgomery 2012;Granatosky et al. 2014). Consequently, we expect to identify a pattern of incomplete convergence for the three taxa of slow arboreal xenarthrans. Conversely, suspensory slow arboreal xenarthrans are hypothesized to extensively converge, following a pattern of complete convergence. To substantiate these hypotheses, we quantified convergence at several levels, employing metrics from the R package 'convevol' (Stayton 2015a) and traced evolutionary trajectories on phylomorphospaces. We believe that, through the integrative examination of a set of traits stemming from different levels of organisation, it is possible to better appreciate intermediate degrees of morphological convergence and related evolutionary patterns.

Raw Data Collection
We visited mammal collections in Germany (Museum für Naturkunde, Berlin; Zoologische Staatssammlung, Munich; Zoologisches Forschungsmuseum Alexander Koenig, Bonn), Austria (Naturhistorisches Museum, Wien), France (Muséum national d'Histoire naturelle, Paris) and USA (Field Museum of Natural History, Chicago, IL; Yale Peabody Museum of Natural History, New Haven, CT) and sampled 47 humeri (43 complete bones and four isolated epiphyses) and 45 femora (43 complete bones and two isolated epiphyses) from 22 xenarthran taxa (17 extant species + 5 extinct genera). Only wild-caught, non-pathological and skeletally mature individuals (i.e., showing completely fused epiphyses) were selected. The femur of Dasypus novemcinctus FMNH 39307 was included because, although showing a partially unfused lateral condyle (accordingly not analyzed in its internal structure), it showed complete fusion at the other epiphyses. Data from both right and left sides were collected. For each humerus and femur, the side was accounted for in the data extraction procedures, to make data from right and left bones comparable (e.g., mirroring 3D meshes, see below). Then, data extracted from right and left bones were pooled and analyzed together.
The specimens were digitized through micro-focus computed tomography (μCT) (Phoenix | X-ray Nanotom, GE Sensing and Inspection Technologies GmbH; XYLON FF35-CT-System, YXLON GmbH; Microtomograph RX EasyTom 150; Nikon XTH 225 ST; GE v|tome|x). Humeri and femora were scanned with resolution ranging from 0.008 mm to 0.083 mm (Online Resources 4 and 5) and image stacks (16-bit tifs) were obtained. Since trabecular analysis is crucially constrained by resolution (Kivell et al. 2011), a visual assessment and the computation of relative resolution (using trabecular parameters) were performed on each specimen (see below).

Bone Orientation and Processing
Specimens were oriented in standard positions using VG Studio Max 3.3 (Volume Graphics, Heidelberg, Germany), placing the x-axis, y-axis, and z-axis along the mediolateral, anteroposterior, and proximodistal directions, respectively. The humeral head was oriented posteriorly, with the X-axis tangent to the head maximal curvature (following Amson et al. 2017), and the femoral head was oriented medially with the Y-axis tangent to the head maximal curvature. For humeri and femora, the centers of both metaphyses (proximal on the top of the stack) were placed on the Z-axis. Isolated epiphyses were oriented through a visual comparison with complete oriented bones of closely related taxa. Fig. S1 (Online Resource 6) shows orienting steps. Oriented specimens were exported as image stacks for diaphyseal and internal epiphyseal analysis (see below).
Beside stacks orientation, 3D surfaces for external shape analysis were generated. Humeri and femora were converted to mesh in VG Studio, then post-processed and simplified in MeshLab (Cignoni et al. 2008) ('Ambient Occlusion'; 'Remove Vertices wrt Quality' with threshold = 5%) and Geomagic Wrap 2017 (3D Systems, Rock Hill, South Carolina, USA).

3D GM
Due to the paucity of anatomical landmarks in long bones, we also used semi-landmarks sliding on curves and surfaces (Gunz et al. 2005) to capture humeral and femoral shapes Fig. 2). Epiphyseal morphology was represented by 21 anatomical + 195 curve semilandmarks for the humerus and by 22 anatomical + 254 curve semi-landmarks for the femur. Anatomical landmarks and curve semi-landmarks (with the latter equally spaced on curves delimited by two of the former, following Gunz et al. 2005) were placed on meshes using MorphoDig 1.5.4 (Lebrun 2018). Once the appropriate number of curve semilandmarks to represent each anatomical curve in the largest humerus and femur was found, it was applied to the rest of the sample. Landmarking was performed in randomized order of specimens and after mirroring left bones, to have specimens from the same side. Because of biasing damage (i.e., missing parts, also considering specimens only represented by epiphyses) or deformation, eight humeri and five femora were discarded from external shape analysis. Features captured by anatomical landmarks and curve semilandmarks are detailed in additional methods (Figs. S2, S3 and Tables S1, S2; Online Resource 6). To perform sliding semi-landmarks on surfaces  we defined a template for each studied bone. Choloepus didactylus ZMB Mam-102636 and Bradypus sp. ZMB Mam-33806 were chosen for the humerus and the femur, respectively. To optimize the semi-automated placement of surface semilandmarks prior to the sliding process, the templates were inflated in the diaphyseal region, using Blender (Community 2018). The two meshes were decimated in Geomagic, decreasing the number of triangles. Then, surface semilandmarks were positioned on each triangle vertex on the two templates using MorphoDig. The decimation factor was chosen to have approximately the same number of vertices (and, accordingly, surface semi-landmarks) for the humerus and the femur (n = 504 and 533, respectively) evenly covering the bones. For the humeral and femoral templates, semi-landmarks were not positioned in the olecranon fossa and the intercondylar fossa, respectively, since pronounced concavities may cause uneven distribution of surface semilandmarks during sliding. The two template models are provided in online supplementary material (Online Resources 7 and 8).
Surface semi-landmarks were projected from the templates to the whole sample in R (R Core Team 2020), and the result was visually checked for each specimen ('placePatch' and 'checkLM' functions in Morpho package; Schlager et al. 2020). We ran sliding of curves and surfaces semi-landmarks minimizing bending energy of a Thin Plate Spline between the specimens and the Procrustes consensus. This procedure was run iteratively using the Procrustes consensus from the previous iteration ('slider3d' function, iterations = 20, stepsize = 0.5, recursive = TRUE, tol = 1e-8, Morpho package; Schlager et al. 2020). Through generalized Procrustes analysis and principal components analysis (PCA) ('gpagen' and 'gm.prcomp' functions, geomorph package;Adams et al. 2021), the shape information was extracted and represented as PC scores. For the humerus and the femur, the first ten PC scores (cumulatively explaining the 92% and 87% of variability, respectively; Tables S3-S4 in Online Resource 9) were used in the subsequent steps.

Cross-sectional Properties
Oriented stacks were imported in Fiji (Schneider et al. 2012). We focused on the diaphyseal region ranging from 30% to the 70% of the total bone length, as measured from the proximal end (hereafter referred to simply as the diaphysis; Fig. 2). This interval represents the maximum diaphyseal length, allowing epiphyseal elements to be excluded from all specimens. CSP were computed in the diaphysis with a slice-by-slice approach (Amson 2019). The Fiji macro from Amson (2019) performs the automated slice-by-slice computation of global compactness and total cross-sectional area on oriented stacks. We modified this macro by adding the computation of non-directional CSP ('Slice Geometry', BoneJ plugin; Doube et al. 2010; https:// bonej. org/ slice) and we ran it after thresholding ('Optimise Threshold > Threshold Only' routine) and purifying ('Purify' routine). The modified macro is provided as Online Resource 2.
We accounted for two types of bias potentially affecting CSP. First, small porosities occasionally connecting the medullary cavity to the external region of the bone may cause drastic and local misestimation of some parameters (e.g., global compactness; Amson 2019). This type of bias is recognizable identifying spikes in the compactness profile (generated by the macro of Amson 2019) corresponding to single slices or intervals that were not accurately analyzed. Moreover, in fossil specimens, the presence of non-bone material within the medullary cavity affects some CSP. This second type of bias is recognisable using both the compactness profile and a visual assessment of stacks. The position of biased slices was marked in a second run of the macro (using the procedure developed by Amson 2019). To cope with these biased slices, we replaced the biased value in R (or series of values if several neighbouring slices were biased) based on the values of correctly analyzed neighbouring slices (Amson 2019). The process of correction (on isolated or small series of slices) was necessary for 20 humeri and 16 femora. Among them, three humeri and one femur exhibited the biased region at one of the extremities of the diaphysis, which prevented building a sequence based on neighbours. In these cases, the extreme slices were manually restored in Fiji ('Paintbrush' tool), using the original non-binarized slice as a reference to properly restore the binarized one. For each CSP, we extracted an average diaphyseal value and the value from the 50% bone length slice (the latter being considered the most informative level in mammalian limb bones; Laurin 2004). Accordingly, we use the subscripts Parameter Aver and Parameter 50 hereafter. Some specimens were damaged or filled with non-bone material for a preponderant length of the diaphysis. If the 50% diaphyseal level was preserved, we took data only from this level. It occurred for five humeri and three femora. Among them, in two humeri and two femora, some minor cracks were present at the 50% level; they were manually repaired before computing CSP on the 2D section. If even the 50% level was dramatically damaged/filled, the specimens were discarded (necessary for one humerus and two femora). All the aforementioned correction procedures are shown in Fig. S4 (Online Resource 6), and the specimens involved are reported in Online Resources 4 and 5.
We analyzed the global compactness (ResC, %), the second moments of area around the minor and the major axis (Imax and Imin, respectively, mm 4 ), the bone cross-sectional area (CSA, mm 2 ), and the cross-sectional shape (CSS, Imax/ Imin, no unit). These parameters provide information about long bone preponderant mechanical properties (ResC), resistance to axial loadings (CSA), bending rigidity (Imax and Imin), and load history (CSS) (Crowder and Stout 2011;Musy et al. 2017; Parsi-Pour and Kilbourne 2020). Parameters related to torsion rigidity (e.g., polar second moment of area, section moduli) were not included because their examination is recommended only for sections that do not significantly deviate from a circle (Daegling 2002). The presence of processes along the diaphysis in a substantial part of the sample (e.g., armadillos) violates this assumption.

Trabecular Architecture
Using oriented stacks imported in Fiji, two and three spherical volumes of interest (VOIs) of epiphyseal trabecular bone were extracted from each humerus and femur, respectively. VOIs were centered in the humeral head, capitulum, and femoral head, as well as in the lateral and medial femoral condyles (Fig. 2). Hereafter, parameters from the proximal epiphyses are referred to with the subscript Parameter prox . For the distal epiphyses, we use Parameter dist for the humerus and Parameter med.con and Parameter lat.con are to specify which condyle of the distal femur.
The largest sphere encompassing only trabecular bone was extracted from each articular surface using a Fiji macro (Online Resource 3). Once the 2D section corresponding to the 50% level of the proximodistal length of the articulation of interest was found, a rectangle, bounding the articular surface of interest, was manually drawn ('Rectangle' Fiji tool). The macro extracts a spherical volume of cancellous bone of given diameter with the same center as the rectangle. We preferred the sampling of the largest VOI instead of constant or scaled VOI sizes (e.g., Ryan and Shaw 2012) to maximize the representation of cancellous bone in smallsized xenarthrans (showing limited number of trabeculae). VOI diameters range from 1 to 26 mm (Online Resources 4 and 5). The femoral head of specimens from the genera Myrmecophaga and Tamandua show a particularly deep fovea capitis. To avoid introducing a bias in the trabeculae to be sampled, we extracted a hemispherical VOI for the femoral head, representing the lateral half of the sphere initially extracted (Fig. 2). The details of VOI extraction are provided in Fig. S5 (Online Resource 6).
Some fossil specimens showed damaged regions, large areas of broken/absent trabeculae and/or predominant presence of intertrabecular non-bone material which prevented the recognition of trabeculae. These VOIs (from three humeral heads, six capitula, three femoral heads, seven lateral condyles, and seven medial condyles) were excluded from the analysis. The possibility to downscale these VOIs in order to include only preserved and 'clean' trabeculae was considered. However, to sample an analogous region, this process requires subsequent downscaling of all VOIs to the same scaling factor (e.g., Amson and Nyakatura 2018a). We refrained from doing so, because it would have resulted in an even lower sample size, due to the exclusion of VOIs not showing a representative number of trabeculae (< 50; see below). For other fossil specimens (four humeral heads, two capitula and two proximal femora), only a moderate degree of non-bone material was present within intertrabecular spaces, and it was clearly distinguishable. It allowed a manual 'cleaning' of trabecular spaces to be performed in a reliable way (as detailed in Fig. S6, Online Resource 6). Excluded and 'cleaned' specimens are reported in Online Resources 4 and 5. All the other VOIs were automatically thresholded and purified. Using the corresponding BoneJ routines, seven trabecular parameters were computed: degree of anisotropy (DA, no unit), trabecular thickness (Tb. Th., mm), connectivity (Conn., no unit, only used to assess the number of trabeculae), connectivity density (Conn.D., i.e., Conn/Total Volume (TV), mm −3 ), bone volume to total volume (BV/TV, no unit), bone surface to total volume (BS/ TV, mm −1 ) and average branch length (Av. Br. Len., mm). The latter was acquired once the stack was skeletonized ('Skeletonise 3D' routine). Results for Conn.D., BV/TV, and BS/TV were corrected considering a spherical volume (or hemispherical volume for the VOI of a femoral head), since a cubic VOI is the default in BoneJ.
Using trabecular parameters, it was possible to assess the quality of CT-scans, computing the relative resolution (Rel.Res = Tb.Th/scan resolution; Sode et al. 2008). The entire dataset has an average Rel.Res of 8.30. All VOIs show Rel.Res ranging from 4.41 to 19 (Online Resources 4 and 5). Thus, they are appropriate according to the minimum limit recommended by Sode et al. (2008) and Kivell et al. (2011). Only the femoral head VOI of Prepotherium sp. YPM PU-15345 lies out of this range (Rel.Res = 3.72), but we visually validated its quality (assessing resolution and contrast). Following Mielke et al. (2018b), VOIs representing less than 50 trabeculae (with number of trabeculae approximated by Conn., Online Resources 4 and 5) were not further analyzed.

Repeatability Assessment
Several steps involved in data extraction unavoidably imply subjectivity (e.g., landmarking, orienting). Thus, on a subsample of seven humeri and seven femora, data extraction was repeated three times in order to assess repeatability. As for 3D GM analysis, we compared x, y, and z non-superimposed coordinates of anatomical landmarks and curve semi-landmarks among the three repeats, while for CSP and trabecular data, we analyzed results for single parameters obtained in Fiji three times per specimen. With pairwise Pearson's correlation analyses, we tested how strongly data from each repeat correlated with the other rounds. All correlation coefficients (r) exceeded 0.85 (all p-values < 0.01), and for the subset tested for repeatability, data from the second repeat were taken for subsequent analysis. Only Av. Br. Len in proximal femora VOIs showed high fluctuations in results (r min = 0.19, p max = 0.67) and was accordingly excluded from the analysis. Data for repeats are provided in Online Resources 10,11,12,13,14,15,16,17,18,19 and 20 and results in Tables S5, S6, S7 and S8 (Online Resource 9).

Body Mass Proxy
All the variables investigated here potentially correlate with body mass in mammals (Egi 2001;Jungers et al. 2005;Doube et al. 2011). It may dramatically affect results when the body mass range of the sample covers several orders of magnitude. Hence, the effects of body mass were isolated from ecological ones (see below). For each humerus and femur, the log-transformed centroid size, reportedly representing a reliable measure of the specimen size (O'Higgins and Jones 1998), was taken as a proxy of body mass. The centroid size is computed as the squared root sum of the squared distances from each landmark to the centroid, and it was extracted from superimposed configurations (resulting from generalized Procrustes analysis) of anatomical landmarks + curve and surface semi-landmarks. Some damaged/ incomplete specimens were discarded from 3D GM analysis, preventing centroid size to be quantied, though they provided CSP and trabecular architecture data. To inform the analysis of CSP and trabecular parameters with a body mass proxy, the centroid size was estimated ('lm' and 'predict' R functions) from a set of metric measurements taken on each epiphyseal articulation from which a VOI of cancellous bone was extracted. The procedure employed to estimate centroid sizes starting from epiphyseal metrics is detailed in Section S1 (Online Resource 6).

Time-calibrated Phylogeny
For each of the parameters introduced above, we tested the correlation with locomotor ecology. It allowed us to preliminarily identify variables that set apart slow arboreal xenarthrans. To cope with statistical non-independence of observations caused by phylogenetic relationships, the univariate tests of significance (see below) were phylogenetically informed. We used the Maximum Clade Credibility DNA-only node-dated mammal phylogenetic tree of 4098 species taken from the posterior distribution of Upham et al. (2019), which is devoid of polytomies. This tree was adapted to our dataset in Mesquite (Maddison and Maddison 2019) by pruning all non-xenarthran clades and including tips not present in the Upham et al. (2019) phylogeny, which included extinct sloths, Dasypus septemcinctus, and Tamandua mexicana. The latter two were included following the most recent inferences on their phylogenetic relationships and divergences (Feijó et al. 2019;Casali et al. 2020). For extinct sloths, we combined data from Bargo et

Locomotor Ecology Assignment
Locomotor ecology was considered a categorical variable characterized by two possible states, non-slow arboreal or slow arboreal. Slow arboreal was assigned to B. tridactylus, B. variegatus, Ch. didactylus, Ch. hoffmannii, and Cyclopes; non-slow arboreal was assigned to all the other taxa. The assignment of extinct sloths here studied to the non-slow arboreal category was based on previous ecological inferences for these taxa. Santacrucian sloths were initially inferred to be terrestrial (Scott 1903) and subsequently reconstructed as semiarboreal/arboreal with climbing habits (White 1993(White , 1997. More recently, through morphometry and functional indices Toledo et al. 2012), muscle reconstruction (Toledo et al. 2013(Toledo et al. , 2015, as well as descriptive morphology and paleosynecological inference (Toledo 2016), the ecomorphotypes of Santacrucian sloths were shown to be clearly distinct from those of 'tree sloths'. These works have indicated a terrestrial-scansorial locomotor ecology, with digging and climbing adaptations, resembling pangolins, Tamandua, and Myrmecophaga Toledo et al. 2012Toledo et al. , 2013Toledo et al. , 2015. Among Santacrucian sloths, an ecological range has been inferred, from semi-arboreality (e.g., Hapalops and Eucholoeops), with digging adaptations in some cases (e.g., Nematherium and Analcitherium), to almost completely terrestrial with only facultative arboreality (e.g., Prepotherium; Toledo 2016). Moreover, the bone microanatomy of some Santacrucian sloth genera has overall shown patterns contrasting with those of 'tree sloths'. Analysing bone cortical compactness, Montañez-Rivera et al. (2018) have found some similarities between Santacrucian sloths and 'tree sloths'. However, investigating the same trait in a wider sample of Santacrucian sloths Alfieri et al. (2021) have recognized a clearly distinct pattern compared to 'tree sloths'. In addition, the analysis of forelimb trabecular and diaphyseal structure of Amson and Nyakatura (2018a) have suggested that Santacrucian sloths and 'tree sloths' occupied different ecological niches. Pending more detailed ecological characterisations, the aforementioned studies clearly highlight that Santacrucian sloths lacked typical adaptations of 'tree sloths' and occupied distinct ecological niches. For the focus of the present analysis, the non-slow arboreal locomotor ecology is ascribed to these extinct sloths.

Univariate Analysis
The relationship between locomotor ecology and each of the morphological variables was investigated with a series of phylogenetic generalized least square (PGLS) regressions and phylogenetic ANCOVAs. Due to the potential influence of body mass on the investigated parameters, we separated this signal from the ecological one by including the log-transformed centroid size (i.e., body mass proxy, computed as mentioned above) as a covariate in each PGLS (parameter ~ locomotor ecology + body mass proxy). To do so, we preliminarily tested for a possible correlation between the body mass proxy and locomotor ecology ('body mass proxy ~ locomotor ecology' PGLSs) for both humeral and femoral data. Since no significant body mass-locomotion correlation was found (p hum = 0.27 and p fem = 0.20), we could reasonably separate ecological effects from body mass effects in each single-parameter PGLS. For each variable, we preliminarily ran two PGLSs, one using raw values and another using log-transformed values. Log-transformation of PC scores that sometimes include negative values was made possible by adding a constant value (i.e., minimum variable value * 1.0001) to the raw values. Results of both PGLSs for each variable allowed us to assess how the distribution of residuals in the two conditions deviate from normality (visually and through Shapiro-Wilk tests). Consequently, for each parameter, we took raw or log-transformed values according to the condition that yielded residuals distribution closer to normality. It resulted in PGLSs with raw values for PC scores and PGLSs with log-transformed values for CSP and trabecular parameters. Variables and body mass proxy were averaged for each species. Individuals of extant taxa only catalogued at the genus level (one humerus and one femur) were excluded. For each PGLS, we pruned the time-tree tips according to the represented species (which varied depending on specimen preservation) using the 'drop.tips' R function in the package 'ape' (Paradis et al. 2004). All PGLSs were performed through the "gls" function ("nlme" R package; Pinheiro et al. 2020) and phylogenetically informed with Pagel's lambda (λ) ("corPagel" function, "ape" R package; Paradis et al. 2004). By default, the maximum likelihood ('ML') method was set. When the fitted linear model failed to reach convergence, the restricted maximum likelihood ('REML') approach was chosen (Table 1). Mean taxa results for variables significantly discriminating slow arboreal xenarthrans (Table 1) were visualized using boxplots (Figs. S7, S10 and Online Resource 9). Distribution of the variables significantly correlated with the size proxy are additionally discussed based on boxplots made with the residuals of a linear regression of log-transformed values against log-transformed centroid sizes (Figs. S8, S11; Online Resource 9). This was not done for PC scores, for which such operation would be partly redundant with the Procrustes superimposition. Boxplots were generated through the "ggplot" function ("ggplot2" R package; Wickham 2016). Significant differences in humeral/femoral PC scores were visualized with bivariate scatterplots if more than one PC score was significantly correlated with locomotor ecology (Figs. S9 and S12; Online Resource 9). Morphological variation captured by 3D GM analysis was visualized warping meshes to minimum and maximum PC scores ('warpRefMesh' function, 'geomorph' R package; Adams et al. 2021) (Figs. S9, S12, S13 and S16; Online Resource 9). Also, differences between PC score extreme values were highlighted in terms of single landmarks coordinates by plotting extreme configurations and drawing the corresponding vectors between them ('deformGrid3d' function, 'Morpho' R package; Schlager et al. 2020) (Figs. S14, S15, S17 and S18; Online Resource 9).

Convergence Analyses
To test if the variables significantly correlated with locomotor ecology evolved through convergence and to provide evidence for the expected models of convergence, we used the pattern-based approach of the 'convevol' R package (Stayton 2015a). The method computes several indices measuring degree of convergence by comparing phenotypic distances between convergent taxa with past distances estimated with ancestral state reconstruction assuming Brownian motion (BM). The method is based on the a priori definition of taxa that are expected to converge. Since we were interested in assessing the convergence resulting in morphological similarity between (and not within) the two lineages of 'tree sloths', data for the species of Bradypus and Choloepus were averaged for each genus, and the time-tree was modified accordingly, obtaining a single tip per genus.
To estimate how similar putatively convergent taxa are compared to maximum distance in their lineage, C1 is computed. It compares the maximum phenotypic distance in the evolutionary history of two converging taxa with their current phenotypic distance in order to establish how much of the maximum distance has been decreased by convergent evolution. C1 is expressed as a proportion (i.e., the ratio between the current distance and the maximum distance is subtracted from 1), ranging from 0 (no convergence at all or divergence) to 1 (lineages that have perfectly converged, resulting in indistinguishable states of the analyzed traits). This property makes C1 a powerful tool, allowing comparisons of convergence strength among different traits (i.e., with different units). On the other hand, it implies that C1 does not account for the absolute amount of evolution due to convergence. As exemplified by Stayton (2015a), two taxa might converge by halving their maximum past phenotypic distance through small or large phenotypic changes yet have the same value of C1 (= 0.5). For this reason, C2 is introduced, in which the current phenotypic distance is subtracted from the maximum phenotypic distance. C2 is thus expressed in the same units employed for distance calculation and directly reflects absolute evolutionary changes. In order to measure the amount of total evolutionary change that occurred in the lineages leading to the convergent taxa can be explained by the studied convergence, C3 is introduced. C3 is computed by dividing C2 by the sum of all the phenotypic changes that occurred along the lineages leading Table 1 Summary of univariate PGLSs and phylogenetic ANCO-VAs [variable ~ locomotor ecology + body mass proxy]. Locomotor ecology and body mass proxy (abbreviated as Loc Ecol and BMP, respectively) p-values inform on the significant correlation of each trait with these two factors. P-values for testing normal distribution of each PGLS residuals are shown in Shapiro (resid). The method used in the gls() function (either maximum likelihood, ML, or restricted maximum likelihood, REML) and the fitted λ value for each PGLS are also shown. Parameters significantly affected by locomotor ecology are highlighted with grey rows. The locomotor ecology p-values allowing significant correlation to be identified are in bold. Abbreviations: Av. Br. Len: Average Branch Length; BS/TV: bone surface to total volume; BV/TV: bone volume to total volume; Conn.D: con-nectivity density; CSA: cross-sectional area; CSS: cross-sectional shape; DA: degree of anisotropy; Imax: maximum second moment of area; Imin: minimum second moment of area; PC: principal component; Tb.Th: trabecular thickness. Subscripts: PC1-10 indicates which principal component from those extracted from humeral and femoral external shape data (only the first ten are taken); CSP 50 and CSP Aver , respectively, indicate CSP extracted from a single 2D slice corresponding to the 50% level of the entire bone length or averaged from those extracted on the entire diaphysis with a slice-by-slice approach; Parameter prox , Parameter dist , Parameter lat.con and Parameter med.con , respectively, indicate trabecular variables quantified in the proximal epiphyses or distal epiphyses / lateral or medial condyle in the femur. Further details are given in Materials and Methods from the common ancestor of the converging taxa to them. Finally, to establish how important the identified convergence process is compared to the total amount of evolution that occurred in the clade to which all convergent taxa belong, C4 is computed. C4 is computed by dividing C2 by the sum of all the phenotypic changes in the clade defined by the common ancestor of the converging taxa and not only along the lineages leading to converging tips (as for C3). See Stayton (2015a) for further details on these indices. For both the humerus and femur, a multivariate convergence analysis was run separately for each of the anatomical levels (3D GM, CSP 50 , CSP Aver , proximal epiphyseal and distal epiphyseal trabecular data) including only variables found as significantly correlated with locomotor ecology (with PGLS, see above). In femoral medial condyle trabecular data, only one variable showed significant relationship with locomotor ecology (i.e., DA; see Table 1 and below). In this case, a univariate convergence analysis was run. Furthermore, we pooled all variables significantly correlated with locomotion in two datasets, one for the humerus (Pheno.sig hum ) and one for the femur (Pheno. sig fem ). Multivariate convergence analyses on sets of variables with different units were preceded by the standardization and centering of each variable ('scale' R function). C1-C4 indices and each associated p-value (C1 p-C4 p, based on evolutionary simulations under BM with 1000 simulations) were computed for each convergence analysis. For the multivariate convergence analyses, we used the modified version of 'convevol' functions developed by Zelditch et al. (2017) ('convrat' and 'convratsig' functions). For the univariate dataset of DA med.con , we adapted the 'convevol' functions used for multivariate analyses to univariate ones ('calcConv1dZelditch' and 'convSig1D' functions; see Online Resource 1).

Humerus
We performed the complete humeral and femoral convergence analysis (steps detailed above) four times. First, we quantified morphological convergence defining slow arboreal xenarthrans as putatively convergent taxa (i.e., Choloepus, Bradypus, and Cyclopes). Then, to isolate the relative contributions to the general pattern, convergence strength was quantified only in suspensory slow arboreal species (i.e., Bradypus and Choloepus). Finally, we performed two additional convergence analyses: one to quantify how strongly Cyclopes converges with Bradypus and another to measure convergence between Cyclopes and Choloepus.
To visualize convergence patterns, we ran PCAs on Pheno. sig hum and Pheno.sig fem . and generated bivariate phylomorphospaces of pairs of the first three PCs (cumulatively explaining around the 90% of variability for both Pheno.sig hum and Pheno.sig fem ; Tables S9 and S10, Online Resource 9). To avoid confusion with the PCs extracted for external shape analyses (see above), PCs obtained through PCA on Pheno.sig hum and Pheno.sig fem will be referred to as PC pheno . Phylomorphospaces were generated with the 'phylomorphospace' R function ('phytools' package; Revell 2012). Statistical comparisons of morphological disparity, which may help to identify complete or incomplete overlap in morphospaces (Stayton 2006;Grossnickle et al. 2020), were not performed due to the low number of putatively convergent taxa, possibly resulting in low statistical power.

3D Geometric Morphometrics
As shown in Table 1, only PC1 and PC3, within the PCs representing humeral shape, are significantly correlated with locomotor ecology. PC1 explains 44.7% of the humeral shape variation (Table S3 and Online Resource 9) and describes variability mainly related to the general elongation, smoothness, and relative development of processes (Figs. S9, S13, S14 and Online Resource 9). Humeri that score low on PC1 are slender and less robust, with reduced processes (such as tubercles), a low degree of torsion (shown by a trochlea-capitulum long axis that is approximately mediolateral), rounded articular surfaces (both head and capitulum), and relatively wide distal epiphysis (with shallower trochlea; e.g., Fig. 3a, b). PC1 values are significantly lower for taxa that have a slow arboreal locomotor ecology and are not correlated with body size (Table 1). Although the PGLS residuals are normally distributed (Table 1), not all slow arboreal xenarthrans contribute to this correlation. The distribution of PC1 values shows that low PC1 scores distinguish 'tree sloths' but not Cyclopes, which is conversely characterized by high PC1 values (Figs. S7, S9 and Online Resource 9).
PC3 accounts for 7.7% of humeral shape variation (Table S3 and Online Resource 9) and highlights subtler differences (Figs. S9, S13, S15 and Online Resource 9). Low PC3 scores characterize morphologies with a more proximally projecting humeral head, more distal greater tubercle, a proximodistally elongated capitulum, and a trochlea that is highly reduced in the mediolateral direction but quite extended in the anteroposterior direction (e.g., Fig. 3c). Higher PC3 values are associated with a more prominent medial epicondyle and a round capitulum (e.g., Fig. 3d-f).
Excluding the effects of body mass, which significantly correlates with PC3, slow arboreal xenarthrans significantly differ from non-slow arboreal taxa in their PC3 scores (Table 1). 'Tree sloths' exhibit moderately lower PC3 values, while extremely low values set apart Cy. didactlyus (Figs. S7, S8, S9 and Online Resource 9).

Cross-sectional Properties
Accounting for size effects, slow arboreal xenarthrans are discriminated by Imax 50 , Imax Aver , Imin 50 , Imin Aver , CSA 50 and CSA Aver (Fig. 4a-c and Table 1). For all these parameters, 'tree sloths' show clearly lower values compared to all other xenarthrans. Imax, Imin and CSA values for Cyclopes are not as low as in 'tree sloths' since they partially overlap with values in the non-slow arboreal range. However, they fall in the lower part of the non-slow arboreal range, always below the average of non-slow arboreal xenarthrans (Figs. S7, S8 and Online Resource 9).  (Table 1). Indeed, they show significantly more numerous (i.e., higher Conn. D prox ), thinner (i.e., lower Tb.Th prox ), and shorter (i.e., lower Av.Br.Len prox ) trabeculae in the humeral head. Moreover, slow arboreal taxa exhibit less uniformly oriented (i.e., lower DA dist ) and shorter trabeculae (i.e., lower Av.Br.Len dist ) in the capitulum (Figs. S7 and S8; Online Resource 9). There is a wide overlap between the Tb.Th prox values of slow and non-slow arboreal taxa. This is due to non-slow arboreal anteaters, which exhibit low values comparable to those of 'tree sloths' and Cyclopes (Online Resource 4). The 'tree sloth' Choloepus didactylus is an outlier of the slow arboreal range, showing a high Tb.Th prox value. For DA dist (after size correction), Cyclopes exhibits extremely low values. Strikingly, the non-slow arboreal anteater Tamandua mexicana falls within the slow arboreal range for size-corrected values of the five humeral trabecular parameters found as significantly correlated with locomotor ecology (Fig. S8, Online Resource 9).

3D Geometric Morphometrics
Among the PCs representing femoral shape, PC1 and PC3 show a significant correlation with locomotor ecology (Table 1). PC1 explains 39.2% of the femoral shape variation (Table S4, Online Resource 9) and relates to general robustness and diaphyseal curvature (Figs. S16, S17, Online Resource 9). Lower PC1 scores characterize slender and straight femora, with a lower degree of torsion (associated with an approximately mediolateral direction for the long axis of the distal epiphysis), large, sub-spherical and proximally protruding head (i.e., smaller angle of femoral neck), reduced trochanters and a short neck (e.g., Fig. 3g, h).
Slow arboreal xenarthrans exhibit significantly lower PC1 scores. Anteaters occupy the lower part of the nonslow arboreal range. Myrmecophaga tridactyla is within the slow arboreal PC1 range, and T. mexicana has the lowest value of the non-slow arboreal taxa (Online Resource 5, Figs. S10, S12 and Online Resource 9).
PC3 explains 9.8% of the femoral shape variation (Table S4, Online Resource 9) and it accounts for minor differences, mainly regarding the distal epiphysis. Higher PC3 scores identify wide and anteroposteriorly flat distal epiphyses, with a shallow patellar groove and condyles of approximately the same size (Figs. S10, S12, S16, S18 and Online Resource 9).

Cross-sectional Properties
Slow arboreal xenarthrans are significantly discriminated for Imax 50 , Imax Aver , CSA 50 , CSS 50 , CSS Aver , including when accounting for the effect of body mass (significantly correlated to Imax 50 , Imax Aver and CSA 50 ) ( Table 1). Slow arboreal xenarthrans show lower values for these parameters (Figs. S10, S11; Online Resource 9). Size-corrected Imax 50 , Imax Aver and CSA 50 values for M. tridactyla and Tamandua -especially T. mexicana -are intermediate between slow arboreal xenarthrans (with some overlap) and armadillos. Tolypeutes matacus has the lowest value among armadillos for size-corrected Imax 50, Imax Aver , and CSA 50 (even within the slow arboreal range in the case of the latter). Regarding CSS 50 , a striking outlier of the non-slow arboreal range is the armadillo Euphractus sexcinctus. This species has a value of 1.62, lower than Bradypus. However, looking at CSS Aver , E. sexcinctus lies within the range of armadillos. T. mexicana is within the slow arboreal range for both CSS 50 and CSS Aver (Online Resource 5, Figs. S10, S11; Online Resource 9).

Trabecular Architecture
A significant relationship with locomotor ecology was found for DA prox , DA med.con and BV/TV prox (Table 1). Slow arboreal species have significantly less uniformly oriented trabeculae in the femoral head and in the medial condyle (i.e., lower DA prox and DA med.con ). Moreover, their femoral head VOI also comprises less bone (i.e., lower BV/TV prox ). 'Tree sloths' are further distinguished by low values for these three parameters, except for C. hoffmanni, which features higher DA med.con values. Relatively high DA prox values are found in Cyclopes, which is in the non-slow arboreal range, while this species is within the range of 'tree sloths' for BV/TV prox and DA med.con (for the latter, Cyclopes has the lowest values in the sample). For T. mexicana, DA prox is within the slow arboreal range, while its BV/TV prox value is the lowest of the non-slow arboreal distribution (Fig. S10, Online Resource 9).

Slow Arboreal Xenarthrans
Focusing on the variables correlated with the slow arboreal locomotor ecology (see above), we found significant convergence for several sets of traits (Table 2). C1 revealed significant convergence for humeral and femoral CSP (both CSP 50 and CSP Aver ), femoral shape, and distal humerus trabecular parameters. CSP exhibits the strongest degree of convergence in both bones (from 72% in femoral CSP Aver up to 90% in humeral CSP 50 ). C2 results confirmed that the strongest convergence is present in CSP (C2 = 1.28-2.5). Apart from CSP, the strongest significant convergence in slow arboreal xenarthrans is found in their femoral shape (66%), although implying a quite low absolute amount of phenotypic change (C2 = 0.068). Yet, the C3 score and related significance level suggest that 35.5% of the femoral shape evolution that occurred between the common ancestor of pilosans and the three slow arboreal pilosan species can be explained by convergence. For humeral 3D GM data and for proximal humerus trabecular properties, the C2 p-value is significant, but it was not confirmed for any of the other C indices obtained for this set of data. Despite the fact that only some anatomical levels showed significant convergence, the analysis of both humeral and femoral phenotypes of slow arboreal xenarthrans -collectively testing variables associated to locomotor ecology (Pheno.sig hum and Phenosig. fem ) -yielded significant convergence (Table 2).

Suspensory Slow Arboreal Xenarthrans
Quantifying convergence only between Bradypus and Choloepus, some additional anatomical levels provided significant C indices p-values (Table 2). Overall, CSPs are confirmed as the most strongly converging traits in Bradypus and Choloepus (but see femoral CSP Aver , as detailed below). Humeral CSP 50 even reaches the outstanding C1 value of 0.99, implying that the two genera of 'tree sloths' evolved a nearly identical phenotype for this trait. C1 and C2 p-values identify the additional significant convergence of humeral external shape (with C1 and C2 indices yielding a three times stronger convergence compared to the analysis in slow arboreal xenarthrans) and proximal femur trabecular parameters. The overall humeral and femoral morphological data (Pheno.sig hum and Pheno.sig fem ) of Bradypus and Choloepus exhibit significant convergence, with higher C1 and C2 values compared to the former convergence analysis. Although it clearly suggests that the average convergence is stronger if Cyclopes is not considered as a putatively convergent taxon, different patterns are present for the humerus and the femur across different anatomical levels. For the humerus, there is a consistently stronger convergence at every level. The only possible exception to this trend is represented by proximal humerus trabecular parameters, which yielded significant convergence for C2 in the analysis of slow arboreal xenarthrans and lost this significance in the analysis of 'tree sloths'. For the femur, mixed results were obtained. Some anatomical levels show an increased convergence in suspensory slow arboreal xenarthrans (i.e., CSP 50 and proximal femur trabecular parameters, for the latter with significant C1-C3 indices p-values only in the convergence analysis of suspensory slow arboreal xenarthrans). On the other hand, femoral CSP Aver and external shape -which were among the most strongly converging traits in the previous analysis -are associated with lower C1 values. For CSP Aver, only a slightly weaker convergence seems to be involved (C2 is increased), but the decrease of convergence strength for the femoral external shape is more evident since this level even loses significance for C1, C3 and C4 scores, compared with the convergence analysis in slow arboreal xenarthrans.

Contribution of Cyclopes to Convergence in Slow Arboreal Xenarthrans
The silky anteater provides a weaker yet significant contribution to both humeral and femoral morphological convergence in slow arboreal xenarthrans (Table 2). It is corroborated by the generally lower convergence strength for Pheno.sig hum and Pheno.sig fem . Indeed, significant p-values were yielded only for C2 in the analysis of convergence of Pheno.sig hum between Bradypus and Cyclopes and for C1-C2 in the analysis of convergence of Pheno.sig fem between Bradypus and Cyclopes. Moreover, this trend of weaker but significant contribution of Cyclopes to the convergence of slow arboreal xenarthrans is supported by Table 2 Results for the four quantitative convergence analyses performed through the 'convevol' R package on the traits found as significantly correlated with locomotor ecology (see Table 1). Taxa putatively considered as convergent include the slow arboreal xenarthrans (upper left), suspensory slow arboreal xenarthrans (bottom left), Cyclopes and Bradypus (upper right), Cyclopes and Choloepus (bottom right). Traits were grouped by anatomical levels: external shape (3D GM), diaphyseal anatomy (CSP) and epiphyseal internal bone (trabecular structure). Then these traits were grouped in humeral (Pheno.sig hum ) and femoral (Pheno.sig fem ) datasets. Significant p-values are shown in bold. * indicates anatomical levels for which only one trait was significantly correlated with locomotor ecology (i.e., DA for medial condyle trabecular parameters), which were analyzed using univariate convergence analyses partial convergence patterns in some anatomical levels.
For the humerus, the overall contribution of the silky anteater is more clearly weaker relative to the one provided by 'tree sloths'. CSP and distal epiphyseal trabecular properties show consistently lower C1 values. The humeral external shape even yields all C indices exactly equal to 0 (with p-value = 0.999), an outcome diagnostic of divergence. For the femur, though, Cyclopes does not show an evidently weaker contribution to convergence along several levels. Indeed, for some traits (e.g., CSP, external shape), Cyclopes tends to converge in a similar or even stronger way compared to 'tree sloths'. The convergences in femoral shape and CSP Aver between Bradypus and Cyclopes are the strongest of those identified for these anatomical levels (94% and 97%, respectively).

Visualisation of Convergence Patterns On Phylomorphospaces
Phenotypic convergence patterns in the humerus and the femur of Bradypus, Choloepus, and Cyclopes are shown on bivariate phylomorphospaces. Humeral and femoral variables found as influenced by locomotor ecology are overall represented by the first three PCs pheno of the humeral and femoral PCAs (performed on Pheno.sig hum and Pheno.sig fem , respectively). The biplots of humeral PC1 pheno and PC3 pheno (explaining 70.1% and 9.6% of variation, respectively; Table S9, Online Resource 9) and femoral PC1 pheno and PC2 pheno (explaining 49.3% and 24% of variability, respectively; Table S10, Online Resource 9) are shown in Fig. 5. For the humeral data, only the PC1 pheno -PC3 pheno biplot reveals a clustering of slow arboreal species (see below), which is interest for this study. The PC1 pheno -PC2 pheno biplot for humeral data and the other biplots not presented in Fig. 5 are provided in Figs. S19, S20, S21 and S22 (Online Resource 9). In both the humeral and the femoral phylomorphospace (Fig. 5)

Discussion
Tackling humeral and femoral bone morphology in xenarthrans through a multipronged integrative approach, we identified a set of expected traits that distinctively characterize taxa adapted to a slow arboreal locomotor ecology (Bradypus, Choloepus, and Cyclopes). Moreover, we quantified convergence for these traits at different levels of analytical detail and assessed convergence trends and evolutionary trajectories in phylomorphospaces. These results substantiated our hypotheses of incomplete convergence of slow arboreal xenarthrans (Bradypus, Choloepus, and Cyclopes).
The subset of suspensory slow arboreal xenarthrans (Bradypus and Choloepus), on the other hand, is more strongly converging and is consistent with a complete convergence pattern. Before discussing the evidence supporting the proposed convergence patterns, anatomical features for which slow arboreal xenarthrans are set apart are interpreted from a functional perspective.

Distinctive Traits of Slow Arboreal Xenarthrans and Functional Interpretations
Humeral External Shape PC1 from humeral shape data discriminates suspensory slow arboreal xenarthrans. It is mostly associated with traits for which Cyclopes resembles other anteaters. The suspensory locomotor style is associated with a humerus with reduced robustness, few muscle attachments along the diaphysis, as well as a rounded and relatively featureless proximal epiphysis (increasing rotatory movements ;Hamrick 1996;Carter and Beaupré 2007). The presence of a relatively larger and symmetrical humeral head with poorly developed tubercles increases shoulder mobility and extensive rotation at this joint (Mendel 1985a;Rose 1989;Pujos et al. 2007;Nyakatura and Fischer 2010). The relatively wide distal epiphysis with shallower trochlea in the humerus of 'tree sloths' probably relates to high degree of elbow extension (Hubbe 2008) and large supination/pronation of the forearm (Figueirido et al. 2016) (Fig. 3a, b). Humeral torsion, consistently low in 'tree sloths' and high in armadillos (especially C. truncatus), is challenging to interpret. Previous work has related humeral torsion to both phylogenetic and adaptive factors (Evans and Vernon 1945;Pieper 1998). Low humeral torsion is not exclusive to 'tree sloths' among xenarthrans. Through a visual assessment of oriented stacks, we found little humeral torsion in anteaters (especially the most terrestrial species, M. tridactyla) and some extinct sloths (e.g., Hapalops). Thus, we hypothesize that factors other than locomotion explain this feature in xenarthrans. The humerus of Cyclopes recalls a more heavilybuilt morphology with more accentuated processes (tubercles, deltopectoral shelf, epicondyles; Fig. 3c), indicating strong lever arms for muscles that pull (Rietveld et al. 1988;Polly 2007). The relative size of processes along PC1 shows an increasing gradient that directly follows the extent to which digging is prevalent in xenarthran habits (Toledo 2016;Amson et al. 2017). This variation goes from nearly absent in 'tree sloths' (low PC1 scores; e.g., Fig. 3a, b), to increasingly larger in anteaters and extinct sloths (e.g., Fig. 3c, d, f) and, especially, armadillos (high PC1 scores, e.g., Fig. 3e). There is an evident overlap between anteaters and armadillos for PC1 values, but higher PC1 values are reached by the armadillo Chlamyphorus truncatus (Fig.  S9, Online Resource 9). Digging requires accentuated muscle attachments for powerful forelimbs (Hildebrand 1985;Kley and Kearney 2007). Thus, since the silky anteater occasionally exhibits a peculiar hook-and-pull arboreal digging behavior (Montgomery 1983;Hayssen et al. 2012), we can hypothesize that the humeral shape variability captured by the PC1 score of Cyclopes is mainly driven by the presence of digging in locomotor habits of the species. Cyclopes may retain some of the morphological traits typical of anteaters (those related to a higher general robustness) because of phylogenetic inertia, which limited adaptations exclusive to slow arboreality. The study of fossil anteater postcrania should allow this hypothesis to be substantiated, resolving the ancestral state in vermilinguans. PC3 scores mainly identify characteristics that distinguish Cyclopes. This species exhibits a more proximally projecting humeral head and a more distally projecting greater tubercle, features more clearly related to arboreality (Pujos et al. 2007). In the trochlea-capitulum complex, Cyclopes shows unique traits within the xenarthrans here analyzed. The silky anteater exhibits a capitulum that greatly exceeds the trochlea in size. Moreover, the capitulum shows a distinctive proximodistal elongation (Fig. 3c). Shape and size of the capitulum and trochlea are possibly informative of forearm range of motion (Andersson 2004;Figueirido et al. 2015). This is wide in arboreal taxa, implying great extents of pronation and supination (Figueirido et al. 2016), and associated with a capitulum and trochlea of roughly similar size (Figueirido et al. 2016) and a sub-spherical capitulum (Toledo et al. 2013). Such features are exhibited by 'tree sloths' but not Cyclopes. Strikingly, the pattern we found contrasts with the findings of Figueirido et al. (2016), in which distal humerus shape in the silky anteater (analyzed with 2D GM) was associated with arboreal habits. Our result could be explained by a highly stable elbow joint more adapted to directionally stereotyped digging than multidirectional arboreal loadings. This is also supported by the distal humerus features highlighted by PC1 (see above). Interestingly, the proximal ulna generally shows a similar shape in arboreal climbers and diggers, since both the locomotor behaviors necessitate powerful flexion/extension of the elbow (White 1993). This has also been argued for Cyclopes by White (1993) but related to arboreality. Hence, the influence of digging on the shape of the elbow joint of Cyclopes might previously have been underestimated, although more recently considered with regard to ulnar morphology (Toledo et al. 2021). Our findings about the distal humeral morphology point to a more complex adaptive scenario that deserves further investigation. It is evident that, although generally recalling other anteaters in some aspects of humeral shape (features highlighted by PC1), Cyclopes shows some unique features of humeral external shape (evidenced by PC3), possibly indicating divergence.

Humeral Cross-sectional Properties
Imax, Imin, and CSA (both at midshaft and averaged for the whole diaphysis) are lower in slow arboreal xenarthrans (Fig. 4a, b, c), including when accounting for body size. CSA, Imax, and Imin provide general information different set of forces acting on the diaphysis, such as resistance to axial and bending loading (Crowder and Stout 2011). The distinctively lower values for the three parameters in 'tree sloths' possibly reflect their slow arboreal locomotor ecology. Indeed, this result suggests that the humerus of 'tree sloth' is weaker overall, as expected for a species engaging more cautious locomotor behaviors (Demes and Jungers 1993;Marchi et al. 2016). The sum of Imax and Imin (not analyzed in this work) positively relates to bone robustness (Crowder and Stout 2011 and references). Since both Imax and Imin are significantly lower in 'tree sloths', their humeri are less robust overall than those of other xenarthrans. In Cyclopes, the increased robustness and diaphyseal processes highlighted for the external shape (see above) might similarly affect humeral CSP. Indeed, Imax, Imin, and CSA values in Cyclopes are between those of 'tree sloths' and the rest of the dataset. Hence, digging could explain the relatively stronger humerus of Cyclopes within the slow arboreal range.

Humeral Trabecular Architecture
Thicker trabeculae have already been reported in the humeral head of armadillos and associated with intense activity related to terrestrial digging (Amson et al. 2017). Another trabecular parameter often yielding functional signal is DA. This variable is generally related to stereotypy of biomechanical loadings acting on a joint (e.g., Ryan and Ketcham 2002;Su et al. 2013;references in Kivell 2016), with multidirectional sets of stimuli expected to produce lower DA (i.e., more isotropic organisation). The lower DA dist of 'tree sloths' and Cyclopes is consistent with this explanation, since the arboreal environment provides more diverse loading directions (Patel et al. 2013 and references) compared to the terrestrial environment, which is partially (e.g., Myrmecophaga and Tamandua) or fully (armadillos) exploited by non-slow arboreal xenarthrans. As for the other trabecular parameters found here to significantly discriminate slow arboreal xenarthrans (higher Conn.D prox and lower Av. Br. Len prox and Av.Br.Len dist ), a functional interpretation is less clear. Amson et al. (2017) argued that these traits may be related to xenarthran locomotor ecology and found results that generally mirror ours (i.e., armadillos showing lower Conn.D prox and higher Av. Br. Len prox and Av.Br.Len dist compared to 'tree sloths' and anteaters). Importantly, Amson et al. (2017) found that anteaters showed the most extreme ranges for Conn.D prox , Av. Br. Len prox , and Av.Br.Len dist , mirroring our findings for slow arboreal xenarthans. However, vermilinguans showed wide overlap with folivorans, and Cyclopes was grouped with anteaters in the study of Amson et al. (2017). Pending further experimental work on trabecular bone ecological adaptation, the overall agreement between our results and those of Amson et al. (2017) allows us to consider higher Conn.D prox and lower Av. Br. Len prox and Av.Br.Len dist as additional traits discriminating slow arboreal xenarthrans.

Femoral External Shape
PC1 from femoral shape data discriminate slow arboreal xenarthrans, that yield significantly lower scores and show a gracile and straight femur with reduced trochanters (e.g., Fig. 3g, h). This morphology implies hindlimb adaptations for vertical climbing and suspensory locomotion in 'tree sloths' (Mendel 1981(Mendel , 1985bGranatosky et al. 2018;Gorvet et al. 2020) and some of these functions may be associated with the femoral morphology of the arboreal vertical climber Cyclopes. However, also the femur of semi-arboreal anteaters (i.e., Tamandua spp.) have these general morphological traits. Lower PC1 scores of nonslow arboreal xenarthrans correspond to anteaters, with M. tridactyla within the slow arboreal range. This unexpected result may be due to the generally long femur of M. tridactyla, with a third trochanter reduced in the mediolateral direction but extended along most of the diaphysis. It can be hypothesized that lower PC1 values more generally discriminate species showing some degree of arboreal behavior rather than a strictly slow arboreal locomotor ecology. However, this interpretation does not apply to M. tridactyla, the most terrestrial anteater, whose degree of elongation and processes reduction are apparently more pronounced compared to Cyclopes. Armadillos have a broad, robust, short and curved femur with pronounced trochanters (the greater projecting more proximally than the head and the third reaching mid-diaphysis; e.g., Fig. 3k). All these features can be associated with their strenuous digging style Milne and O'Higgins 2012;Marshall et al. 2021). The proximal epiphysis of taxa with low PC1 scores is characterized by reduced trochanters (with the greater extending far laterally), a large, curved and proximally-directed head (i.e., smaller neck angle), and a shorter neck. All these features arguably relate to arboreality. The reduced greater trochanter, together with head morphology and neck angle, allow higher hip abduction and multidirectional femoral rotation (Walker 1974;Fleagle and Meldrum 1988;White 1993;Marchi et al. 2016). Although not captured by PC1, anteaters differ concerning the neck angle. Cyclopes exhibits a smaller angle (i.e., femoral head directed more proximally, e.g., Fig. 3i), a feature shared with 'tree sloths' (e.g., Fig. g, h) that is indicative of arboreal habits (White 1993). The small third trochanter associated with low PC1 values (e.g., Fig. 3g, h, i) indicates reduced action of the hip extensor/flexor muscles attaching on the process. It is also a condition associated with arboreal habits (Sargis 2002). Strong hip flexors may nevertheless be needed for suspensory locomotion in 'tree sloths' (Mendel 1985a;Spainhower et al. 2021), which can be achieved through the modification of muscles not attaching on the third trochanter (Marshall et al. 2021). Similar to what was found in the humerus (see above), the degree of torsion in the femur is challenging to interpret. Through a qualitative evaluation of oriented bones, we identified low degree of femoral torsion in 'tree sloths', extinct sloths, M. tridactyla, and Tamandua species. Strikingly, Cyclopes is the pilosan showing the highest degree of femoral torsion. At present, we exclude a direct relationship with locomotor ecology. Further functional investigations on this feature are warranted.
PC3 from the femoral external morphology data generally identifies features related to the distal epiphysis shape. Slow arboreal xenarthrans, characterized by higher PC3 scores, show a mediolaterally wide but anteroposteriorly flat distal epiphysis. These characteristics can be ascribed to low values of the index of the distal femur shape (White 1993) and potentially set slow arboreal xenarthrans apart. Interestingly, Cylopes was found in the analysis of White (1993) as the taxon with the lowest index of distal femur shape (even lower than 'tree sloths'). In this work, only B. tridactylus yields higher PC3 scores than Cyclopes. Moreover, femora of taxa with high PC3 scores have both condyles that are roughly isometric, with their morphology defining a shallow patellar groove. All of these distal femoral features can be related to arboreality, since this disposition of the condyles reflects a frequently flexed knee joint, the typical condition of suspensory and highly arboreal species (White 1993).

Femoral Cross-sectional Properties
Lower femoral Imax and CSA in 'tree sloths' and Cyclopes ( Fig. 4e-g) are coherent with lower bending and axial strengths hypothesized to characterize long bones in slow and cautiously moving species (Demes and Jungers 1993;Marchi et al. 2016). Although not reaching significance, Imin 50 , Imin Aver , and CSA Aver value distributions are lower for slow arboreal xenarthrans. Conversely, the higher femoral strength (Imax and CSA) exhibited by armadillos is consistent with high stresses deriving from fossoriality and digging Milne and O'Higgins 2012;Marshall et al. 2021). Thus, the hindlimb strength required by locomotor ecology drives the distribution of Imax and CSA, as suggested by intermediate values in non-slow arboreal anteaters. An interesting result to note are the quite low Imax and CSA values shown by the southern three-banded armadillo, T. matacus. Since this species has been reported to be the least fossorial armadillo (Amson et al. 2017 and references therein), we can propose that degree of fossoriality is the main determining factor for femoral resistance to bending and axial loadings. Lower CSS (both 50 and Aver ) in slow arboreal xenarthrans again seems to reflect locomotor aspects. This parameter relates to the degree to which the diaphyseal cross-section is circular. A more circular cross-section (CSS close to 1) is expected in species undergoing more directionally variable loading regimes, while species using their limbs in more stereotyped ways (i.e., one main direction) should have a more elliptical cross-section (e.g., Ruff and Runestad 1992;Carlson 2005;Patel et al. 2013). Lower CSS values in slow arboreal xenarthrans (Online Resource 5) are consistent with more variable regimes shown in suspensory/arboreal quadrupedal species (Carlson 2005;Demes and Carlson 2009;Schmidt and Fischer 2010;Patel et al. 2013). This result mirrors those of Patel et al. (2013) for the humerus of 'tree sloths', which they found was more circular. Our humeral CSS results did not yield significant differences between slow arboreal and non-slow arboreal taxa (Table 1). It is probably due to the quite high CSS values found in Cyclopes or to the quite low CSS values found in some non-slow arboreal extinct sloths (Online Resource 4). The contrastingly low femoral CSS 50 and high femoral CSS Aver retrieved for the six-banded armadillo E. sexcinctus can be explained by the position of the third trochanter along the diaphysis in armadillos. Data from a single cross-section possibly missed the third trochanter, which strongly affects the average degree of circularity. This finding emphasizes the advantages of wholebone CSP analysis (e.g., Amson 2019; Profico et al. 2021), in addition to traditional methods based on single crosssections (e.g., Ruff 2002;Marchi et al. 2016).

Femoral Trabecular Architecture
Of the trabecular parameters computed from the femoral internal epiphyseal structure, DA prox , BV/TV prox , and DA med.con showed a significant relationship with locomotor ecology (Table 1). DA possibly relates to the stereotypy of loadings, and BV/TV probably relates to the magnitude of stresses acting on joints (Kivell 2016 and references therein).
Lower DA prox and DA med.con can be consistent with more variable mechanical requirements in arboreal environment (Patel et al. 2013), resulting in a lower degree of orientation of trabeculae. Moreover, lower BV/TV prox can relate to a cautious and less active locomotor ecology. Accordingly, distinctly lower DA prox , DA med.con , and BV/TV prox in 'tree sloths' and Cyclopes (but see exceptions below) may reflect locomotion. No significant behavioral differences have been reported between the Choloepus species, thus the unexpectedly higher DA med.con in C. hoffmanni is difficult to justify in an ecological context. Another striking result is the higher DA prox of Cyclopes (found in all three individuals analyzed; Online Resource 5). It might mean that the proximal femur is responding more to stereotyped loadings (such as those deriving from digging) in the silky anteater. However, it is difficult to relate higher DA in the femoral head of the silky anteater to digging, since this is performed with forelimbs in Cyclopes (Amson et al. 2017), though the hindlimbs must anchor the body to the arboreal substrate in use. Typical femoral morphological adaptations in xenarthrans more specialized for digging (i.e. armadillos) are mainly related to fossoriality, with hindlimbs playing a role in anchoring the posterior portion of the body to the soil, and the fossorial behaviors are obviously not present in arboreal environment (Marshall et al. 2021). Although hindlimbs of Cyclopes supposedly need to be adapted to provide support for arboreal digging, pending additional experimental data at this stage we cannot frame the higher DA prox of Cyclopes in a locomotor ecological explanation.

Convergence Patterns in Slow Arboreal Xenarthrans
The quantitative analysis fully confirms our expectations in indicating that slow arboreal xenarthrans show convergence, that the convergence between Bradypus and Choloepus is strongest, and that Cyclopes converges to a lesser but significant degree (Table 2). Indeed, the test for convergence is significant in slow arboreal xenarthrans when all variables are combined. In the convergence analysis of suspensory slow arboreal xenarthrans (i.e., Bradypus and Choloepus), both C indices and p-values clearly indicate that the overall convergence is stronger compared to the analysis in slow arboreal xenarthrans, for both the humerus and femur. Separate convergence analyses of Cyclopes and either Bradypus or Choloepus with all variables combined yield significant convergence only for few C indices and only in the convergence analysis between Cyclopes and Bradypus. Nevertheless, significant C index p-values at several anatomical levels clearly show that Cyclopes tends to converge morphologically with both 'tree sloths'. Remarkably, for some femoral anatomical levels, Cyclopes even contributes more to the general convergence pattern in slow arboreal xenarthrans than the convergence between Bradypus and Choloepus. This is especially evident for femoral external shape and CSP Aver , for which the convergence between Cyclopes and Bradypus is stronger than the convergence between Bradypus and Choloepus.
The two bivariate phylomorphospaces (PC1 pheno -PC3 pheno from the PCA on Pheno.sig hum and PC1 pheno -PC2 pheno from the PCA on Pheno.sig fem ; Fig. 5) demonstrate evolutionary trajectories of convergent evolution. Pilosans, to which slow arboreal xenarthrans belong, cluster separately for both the humerus and the femur. Yet, the putatively convergent taxa migrated toward a similar distinct area of the morphospace. In the sub-region defined by Bradypus, Choloepus, and Cyclopes within the morphospace of pilosans, it is evident that Cyclopes shows a different degree of similarity, compared to how Bradypus and Choloepus are similar each other (Fig. 5a, c). Indeed, while 'tree sloths' are closer to each other, Cyclopes occupies a more distant position (distinct from all other xenarthrans in the humeral plot) causing the non-complete overlap of slow arboreal xenarthrans in the morphospace. Hence, Cyclopes is more similar to 'tree sloths' than other anteaters (apart from T. mexicana; see below) but still shows an overall distinctive morphology. The non-complete overlap of the silky anteater with 'tree sloths' can be interpreted as indicative of incomplete convergence (Herrel et al. 2004) in the humeral and femoral phenotype of slow arboreal xenarthrans. This is supported by the fact that convergence explains a minor percentage of the overall evolution in the lineages leading to slow arboreal xenarthrans (C3: 12-23.5%). On the other hand, the more extensive clustering of Bradypus and Choloepus in phylomorphospaces allows us to propose a model of complete convergence for 'tree sloths'.
Complete convergence only occurs when few restricted niches are defined and within them a single optimum is specified by functional constrains. In such cases, there is negligible influence of phylogenetic factors (Zelditch et al. 2017). The highly specialized locomotion of 'tree sloths' and the peculiar constraints that it imposes potentially generate a suitable context for complete convergence and, accordingly, for outstandingly similar morphotypes to converge in two lineages that diverged around 35 Mya (Delsuc et al. 2019). On the other hand, incomplete convergence may be explained by a wide and homogeneous slope instead of a distinct peak in the adaptive landscape, a condition yielding a selective plane that is directional and inclined (Phillips and Arnold 1989;Polly 2004;Grossnickle et al. 2020). Alternatively, a single peak could be present, but it may be quite large, and convergent groups, being far from the peak, might stay in a restricted region of the morphospace (Grossnickle et al. 2020). The latter scenario might account for the convergence pattern observed among slow arboreal xenarthrans. Indeed, adding the silky anteater to 'tree sloths' defines a broader ecological niche and thus a larger adaptive peak. The adaptive landscape, when exhibiting the conditions for incomplete convergence, might be better described by a uniform plane (e.g., Polly et al. 2016). This scenario results in similar/parallel shifts in phylomorphospaces, suggesting that ancestral differences are responsible for the incomplete convergence pattern (Grossnickle et al. 2020). Slow arboreal taxa show evolutionary trajectories from ancestors which are directionally similar (roughly parallel in the femur; Fig. 5c). To some extent, incomplete convergence of Cy. didactlyus may reflect the retention of an ancestral phenotype. Hence, we explain the incomplete convergence in slow arboreal xenarthrans with the locomotor ecology of Cyclopes that shows both independently evolved similarities with 'tree sloths' and plesiomorphic traits (i.e., shared with other anteaters). This intermediate degree of locomotor similarity would have resulted in humeral and femoral morphotypes of Cyclopes that only partially converge with 'tree sloths', yielding an overall unique morphology. Future studies could address morphological evolution in Vermilingua (and, specifically, the lineage leading to Cyclopes) to better understand how the distinctive morphotype of the silky anteater arose.
It is important to notice that the evolutionary trajectory for the overall humeral phenotype of Cyclopes, despite directionally similar and pointing toward the same morphospace region of 'tree sloths', is slightly divergent. As exemplified in Grossnickle et al. (2020), one or few divergent evolutionary trajectories may be at the base of the relatively large morphospace occupied by incompletely convergent clades. We highlighted divergence for one humeral anatomical level of Cyclopes (i.e., external shape). Low 3D GM PC3 values mainly identified traits for which the silky anteater shows a unique (hence, divergent) morphology. Nevertheless, the significant convergence with 'tree sloths' for several levels in the humerus of Cyclopes brings us to believe that a model of incomplete convergence better describes our outcomes.
Future studies should substantiate the models of convergence proposed here by, for example, more thoroughly defining the degree of morphospace overlap using quantitative analyses of disparity or model-fitting analyses Mahler et al. 2013;Mahler and Ingram 2014;Wölfer and Nyakatura 2019;Grossnickle et al. 2020). For the latter, it is noteworthy that the pattern-based approach used here ('convevol', Stayton 2015a) is based on past phenotypic distances which, in turn, are estimated through ancient state reconstructions assuming BM. This model of past phenotype estimation may be overly simplistic in some cases. Thus, further analyses should test alternative models of evolution.

Unexpected Results for T. mexicana
The non-slow arboreal anteater T. mexicana showed unexpected patterns for several traits. These traits, more reminiscent of slow arboreal taxa, all pertain to the internal structural variables (both diaphyseal and epiphyseal). That is the case for all the humeral trabecular parameters found as significantly related to locomotor ecology. Although other non-slow arboreal xenarthrans occasionally showed values compatible with slow arboreal variability for isolated traits (see above), T. mexicana is the only species to show this pattern for an extensive set of traits, which caused T. mexicana to cluster with slow arboreal xenarthrans in phylomorphospaces of Pheno.sig hum and, especially, Phenosig fem (Fig. 5). This striking result is hard to interpret. No major ecological differences were reported between T. mexicana and the other Tamandua species studied, T. tetradactyla (Brown 2011;Hayssen 2011;Navarrete and Ortega 2011;Chiapero et al. 2021). Among their few behavioral differences, one can cite the larger home range of T. tetradactyla (Montgomery 1985;Rodrigues et al. 2001;Brown 2011) and an activity budget of T. mexicana that is characterized by a greater percentage time spent feeding compared to locomotion (Brown 2011). These differences are probably too subtle to justify the high degree of internal structural similarity shared by T. mexicana, 'tree sloths', and Cyclopes. Unexpectedly, T. mexicana was previously found to share a trait of potential functional significance with slow arboreal xenarthrans, humeral and femoral cortical compactness (Alfieri et al. 2021). These findings come from the humerus and femur of a single individual (FMNH 123994), which was included both in our study and that of Alfieri et al. (2021). Further studies including more than one specimen of T. mexicana will be crucial for clarifying this puzzling finding.

Evolvability of Internal and External Structure in the Humerus and Femur
The convergence analyses highlighted different patterns depending on the anatomical level studied. Internal bone structure (i.e., CSP and trabecular architecture) more strongly converges than external shape. Indeed, the latter is one of the levels showing weak/non-significant convergence between slow arboreal xenarthrans (in the humerus) and between suspensory slow arboreal xenarthrans (in the femur). Even when showing strong convergence according to C1 (i.e., femoral shape among Bradypus, Choloepus, and Cyclopes and, especially, between Bradypus and Cyclopes; humeral shape between Bradypus and Choloepus), C2 clearly suggests that high levels of similarity detected by C1 involve very low amounts of absolute morphological evolution. Moreover, external shape (humerus of Cyclopes) is the only level for which we identified divergence (Table 2). These results are in agreement with the expected lower functional and/or evolutionary plasticity of external morphology due to stronger phylogenetic inertia and/or interdependence with other articulating elements (Ruff and Runestad 1992;Lieberman et al. 2001;Currey 2002;Kivell 2016). Concerning internal bone structure, diaphyseal structure is the most strongly converging level, while epiphyseal parameters show a relatively weaker convergence (Table 2). This difference in convergence strength between internal structural levels could relate to the different ways the diaphysis and epiphysis of a long bone are expected to reflect biomechanical loadings and ecological constraints. The two long bone levels relate to bending/torsional loadings (assumed to be higher in the diaphysis) versus axial loadings (arguably higher in the epiphysis) (Biewener et al. 1996;Pontzer et al. 2006;Carter and Beaupré 2007;Barak et al. 2011). Since one of the main features shared by slow arboreal xenarthrans is having their forelimbs and hindlimbs mainly loaded in tension (Preuschoft 2002;Nyakatura and Andrada 2013;Schmitt 2017, 2019;Hanna et al. 2017), one might expect epiphyseal parameters (more responding to axial loadings, i.e., tension/ compression) to be more strongly affected by ecomorphological convergence. Moreover, trabecular bone is reportedly more responsive to biomechanical stimuli and potentially more prone to functional adaptation to ecology and locomotion than diaphyseal cortical bone (Kivell 2016 and references therein).
This study confirms that the use of CSP and trabecular parameters as traits potentially correlated with ecology should be strongly encouraged. Importantly, further analyses of diaphyseal and epiphyseal variables are also needed in order to elucidate their relative contributions to ecomorphological convergence at these related but not necessarily co-varying levels. Expecting different relative contributions of internal diaphysis and epiphysis to convergence only based on prevalent regimes of loadings (i.e., bending vs. axial) may be oversimplicistic. Several questions should still be answered in detail on the factors affecting prevalently trabecular bone functional adaptation (e.g., contractile forces vs. substrate reaction forces, Judex and Carlson 2009;Robling 2009; high-magnitude loadings vs. high-frequency loadings; Bertram andSwartz 1991, Kivell 2016 and references therein for details and other issues). CSP investigation started earlier than trabecular bone investigation, and several studies have established a straightforward relationship with the biomechanical environment (e.g., Jones et al. 1977;van Der Meulen et al. 1993, 1996Ruff et al. 1994Ruff et al. , 2006van Der Meulen and Carter 1995). However, some in vivo works have revealed that diaphyseal internal bone functional adaptation can be more complex than modeled (e.g., Lanyon and Rubin 1985;Judex et al. 1997;Demes et al. 1998;Wallace et al. 2014). Recent analyses of how CSP and trabecular structure have evolved according to ecology using a phylogenetic framework have revealed complex/unexpected adaptative scenarios that can probably be explained by heterogeneous and disparate bone-dependent or trait-dependent models of morphological diversification (e.g., Amson and Kilbourne 2019;Scheidt et al. 2019;Parsi-Pour and Kilbourne 2020). It is another element in support of the need to use model-fitting analysis in future studies, as suggested above. Hopefully further studies may provide additional experimental and analytical evidence regarding the relative responses of CSP and trabecular bone to biomechanical forces and ecological constraints, both through ontogeny and evolutionary history.

Potential Limitations
We recognize that some aspects of this work may represent limitations. To capture humeral and femoral external shape, we used a high-density approach that combines anatomical landmarks with semi-landmarks sliding on curves and surfaces (Gunz et al. 2005;Gunz and Mitteroecker 2013;Bardua et al. 2019). It is not unusual that studies applying high-dimension landmarking (e.g., Cornette et al. 2013;Fabre et al. 2013a;2013b;Aristide et al. 2016;Dumont et al. 2016;Segall et al. 2016;Pavoni et al. 2017) involve a large number of variables (p, proportional to the number of landmarks) compared to observations (n, i.e., analyzed specimens). In these cases, PCA may be problematic; indeed, a larger p/n ratio potentially causes larger first relative eigenvalues and more unstable corresponding eigenvectors (Mitteroecker and Bookstein 2011;Bookstein 2017). This can return spurious results, as artificially created/inflated between-group differences due to random noise that can obscure signal (Bookstein 2017;Cardini et al. 2019). We used 720 and 809 landmarks (including anatomical and curve + surface sliding semi-landmarks) to analyze 39 humeri and 40 femora, respectively. Hence, the possibility that the differences found in this work may result from a biased analysis should be discussed. The recommendations of Bookstein (2017) of controlling p/n in geometric morphometric studies are fundamental to avoid the detection of non-natural between-group differences. On the other hand, it is common that the usage of only anatomical landmarks prevents to capture the morphology of important structures, even within the same clade. This occurs because such morphologies can be challenging or impossible to represent in all taxa (i.e., ambiguous homology; Goswami et al. 2019), especially if the clade is highly diverse morphologically. The problem is exacerbated when analysing elements with scarce homologous and easily detectable landmarks, such as limb bones . In this study, we faced similar issues, since xenarthrans are well known for their high evolutionary morphological plasticity (Billet et al. 2012;Amson et al. 2017;Amson and Nyakatura 2018a, b). For instance, the diaphyseal morphology would hardly be represented by homologuous anatomical landmarks alone (e.g., highly reduced/absent deltopectoral shelf of the humerus and third trochanter of the femur in 'tree sloths' ; Fig. 3). Landmarking protocols that do not reliably capture highly variable morphologies also fall short of one of the main principles of morphometry (Conaway 2021). Hence, high-density landmarking including curve and surface semi-landmarks is often essential to provide complementary information on shape that is not otherwise captured by landmarks alone (Goswami et al. 2019). Using high dimension landmarking and decreasing p/n by increasing sample size would be the most desirable solution. However, reaching sample sizes that meet the criteria suggested by Bookstein (2017) is often not possible given the availablity of specimens in collection-based research, especially if fossils are also studied. To reach a compromise, we believe that a case-by-case evaluation should be performed, also considering that effects of high landmarks dimensionality (compared to standard landmarking) may depend on the analyzed taxon and skeletal element (Conaway et al. 2019).
To evaluate the effect of the large number of landmarks on our results, both humeral and femoral 3D GM analyses were re-run only using anatomical landmarks (n = 21 for the humerus and n = 22 for the femur). For both the humerus and the femur, biplots of the first three PCs (PC1-PC2, PC1-PC3 and PC2-PC3) extracted from the analysis using anatomical + sliding semi-landmarks are compared with the corresponding scatterplots of PCs resulting from the analysis using only anatomical landmarks (Figs. S23, S24, S25, S26, S27, S28 and Online Resource 9). It is evident that major groups and overall discrimination of taxa are preserved. Moreover, we re-ran univariate PGLSs to identify a potential relationship with locomotor ecology (as detailed in Materials and Methods). We used the first ten PC scores for both humeral and femoral data resulting from only anatomical landmarks coordinates (results in Tables S11, S12 and Online Resource 9). Then, we compared these results to those obtained using PC scores resulting from anatomical + sliding semi-landmarks coordinates (Table 1). PC1 for the femur and PC3 for the humerus and the femur are confirmed as significantly related to locomotor ecology even if only anatomical landmarks are used. The only difference concerns PC1 for humeral data, for which the correlation is not significant when using only anatomical landmarks (p-value = 0.25; Table S11, Online Resource 9). Using anatomical landmarks + sliding semi-landmarks, we found humeral PC1 to overall discriminate 'tree sloths' from all other xenarthrans. Some of the related anatomical features distinguishing the humerus of 'tree sloths' identified here (e.g., highly reduced/absent processes along the diaphysis; Fig. 3) are undoubtadbly better represented by including information from curves and surfaces. Thus, as mentioned above, using only anatomical landmarks potentially entails loosing fundamental information and, in this case, a lower degree of discrimination shown by 'tree sloths' for humeral PC1 (Figs. S23, S24 and Online Resource 9). Given the consistent distribution of species on PC scatterplots and PC score correlations with locomotor ecology, we view the biasing effect of noise deriving from high p/n in this work as negligible.
Another aspect of the methodology herein employed that can be viewed as a limitation is our choice to prune variables to only include the ecologically affected ones in the convergence analyses (i.e., those that significantly discriminate slow arboreal xenarthrans through PGLS). Indeed, one might object that by doing so, the probability to show strong convergence through the pattern-based approach of 'convevol' (Stayton 2015a) is artificially increased. We can frame the present study as a 'morphometry-driven approach to ecomorphology' (sensu Feilich and López-Fernández 2019) in that morphological variation is quantified, the variation among studied taxa is described, and we focused on hypotheses concerning the correlation of form and function. Key elements of this approach are represented by those aspects of morphology that should be investigated and how to use them in subsequent analyses to best answer the experimental question (Feilich and López-Fernández 2019). As theorized by Feilich and López-Fernández (2019), the selection of measurements to be analyzed, in turn, reflects different approaches, within the 'morphometry-driven' approach. Choosing only strongly ecologically-related functional variables, as we did when identifying ecological relatedness through PGLS, means following the 'functionally-informed approach'. Using variables directly related to functional modelling (as we discussed above in '29'), this generally allows one to discuss the functional implications of the identified patterns (Feilich and López-Fernández 2019). Several past studies have followed our same approach of excluding variables preliminarily identified as independent of ecology from downstream multivariate and/or phylogenetic comparative analyses (e.g., Chen and Wilson 2015; Amson and Nyakatura 2018a;Schwab et al. 2020;Sosiak and Barden 2021).
We ran another convergence analysis including all variables (Table S13, Online Resource 9). Most of the convergence patterns previously identified are still significant but weaker: slow arboreal xenarthrans show convergent humeral and femoral CSP, femoral shape, and humeral overall phenotype (i.e., all anatomical levels pooled together); Bradypus and Choloepus show convergent humeral and femoral CSP, humeral shape, and humeral overall phenotype; Cyclopes and Bradypus show convergent femoral shape. Some levels show weaker convergence, with p-values for their C-scores loosing significance (i.e., distal humerus trabecular parameters and femoral overall phenotype among slow arboreal xenarthrans; distal humerus trabecular parameters and femoral overall phenotype between Bradypus and Choloepus; humeral CSP between Cyclopes and 'tree sloths'). Yet, several levels show strikingly stronger convergence (i.e., humeral shape, proximal humerus and proximal femur trabecular parameters in slow arboreal xenarthrans; proximal humerus trabecular parameters and femoral shape between Bradypus and Choloepus). Overall, one can identify the general patterns which substantiated the expected modes of convergence. Slow arboreal xenarthrans show significant convergence at several levels. The convergence strength increases on average for the analysis of suspensory slow arboreal xenarthrans, supporting their pattern of complete convergence. The separate contribution of Cyclopes to the convergence in slow arboreal xenarthrans is clearly minor although significant, supporting the pattern of incomplete convergence with Choloepus and Bradypus (Table S13, Online Resource 9). The increase of convergence strength for some levels when all variables are analyzed suggests that, in the present study, including only the subset of variables showing significant relationship with ecology does not merely correspond to an artificial increase in the statistical significance of the downstream convergence analysis. It rather reflects that including all variables or only those strongly influenced by ecology in ecomorphological studies (as we did with the main convergence analysis) represent two distinct theoretical approaches (as defined by Feilich and López-Fernández, 2019).

Conclusions
In this work, we integratively examined data collected from external shape and internal structure (diaphyseal and epiphyseal) of the humerus and the femur of several species of xenarthrans. Through phylogenetic comparative methods, we highlighted previously unrecognized patterns of ecomorphological convergence in slow arboreal species: 'tree sloths' (Bradypus and Choloepus) and the silky anteater (Cyclopes). A set of traits shared by all or some slow arboreal xenarthrans were identified. For most of these characteristics, it is possible to derive a functional interpretation. We then assessed the patterns of convergence in the features recognized as correlated with locomotor ecology. The presence of ecomorphological convergence in slow arboreal taxa was quantitatively shown, with Bradypus, Choloepus, and Cyclopes occupying a distinctive area of humeral and femoral phylomorphospaces. Slow arboreal xenarthrans mainly converge in their humeral and femoral internal structure, especially at the diaphyseal level. Bradypus and Choloepus exhibit a higher degree of morphological similarity, while Cyclopes falls between these taxa and anteaters. This trend results in non-complete overlap of slow arboreal xenarthrans in the morphospace defined by the studied parameters. This led us to propose a model of incomplete convergence for slow arboreal xenarthrans. In contrast, our results for 'tree sloths' are interpreted as evidence for complete convergence between suspensory slow arboreal xenarthrans. In future works, the integrative and quantitative analysis of morphological convergence could be extended to non-xenarthran slow arboreal species or to other ecological adaptations. A wider methodological (e.g., including quantitative disparity comparison and/or evolutionary modelfitting analysis) and taxonomic framework may further elucidate key mechanisms underlying ecomorphological convergence. distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.