Abstract

INTRODUCTION

In tetrapods, limb bones play a dual role in supporting the body and ensuring efficient locomotion in a given environment, while being intimately related to the muscles attached on them. Therefore, the shape of the limb bones of quadrupeds is strongly related to body size, body mass and locomotor constraints (Hildebrand, 1974;Polly, 2007;Biewener & Patek, 2018). The strong tendency of many quadruped lineages to converge towards high body mass across their evolutionary history (Cope, 1887;Depéret, 1907;Raiaet al., 2012;Bakeret al., 2015;Bokmaet al., 2016) highlights repeated patterns of musculoskeletal changes related to increase in size and weight for diverse morphologies. The condition of heavy animals is often called ‘graviportality’ and classically opposed to ‘cursoriality’, which is encountered in running animals like horses and many other ungulate taxa (Hildebrand, 1974;Carrano, 1999). Many works have investigated the skeletal features often found in ‘graviportal’ tetrapods, such as more vertical and thicker limbs, a re-orientation of girdle bones, changes in limb segment proportions or an internal micro-anatomical restructuration (Gregory, 1912;Osborn, 1929;Hildebrand, 1974;Coombs, 1978;Eisenmann & Guérin, 1984;Biewener, 1989a,b;Bertram & Biewener, 1990;Houssayeet al., 2016). These changes in body proportion are also linked with changes in locomotor capacities. All these modifications can lead to a high diversity of body plans for a single given mass (Hildebrand, 1974;Polly, 2007) and, consequently, to modifications of the bone shape. However, bone shape changes associated with body mass increase are poorly documented among quadrupeds.

Although only five species of modern rhinos survive today (Dinerstein, 2011), the Rhinocerotoidea was a flourishing superfamily during most of the Cenozoic. A rich and well-preserved fossil record led to the description of more than a hundred species distributed in Eurasia, North America and Africa, occupying a huge diversity of ecological niches and showing a great variety of locomotor conditions (Prothero & Schoch, 1989;Cerdeño, 1998;Prothero, 2005;Biasattiet al., 2018). Rhinocerotoidea ranged from less than 100 kg inHyrachyus, the most ancient representative of the superfamily (Antoine, 2002;Baiet al., 2017), to between 10 and 17 tons in giant Paraceratheriidae (Fortelius & Kappelman, 1993;Prothero, 1998,2013;Qiu & Wang, 2007;Larramendi, 2016) (Table 1). Between these two extremes, numerous lineages showed convergent increases in body mass, with many species exceeding 1 ton or more (Cerdeño, 1998). In addition to this variation in body mass, the evolutionary history of rhinocerotoids exhibits fluctuations in their general body plan (from cursorial to graviportal), their degree of brachypody (or gracility, i.e. reduction of their relative limb length), their ecological affinities (from open environments to presumed semi-aquatic lifestyles), their number of forelimb digits (tetradactyl or tridactyl manus), the presence, position and number of horns and the size of their head, all of which may also have co-varied with the shape of long bones (Guérin, 1989;Prothero & Schoch, 1989;Cerdeño, 1998;Prothero, 1998,2005,2013;Antoine, 2002;Becker, 2003;Beckeret al., 2009;Baiet al., 2017).

Consequently, members of the superfamily represent a rich diversity of body mass, size and proportion and constitute a great example for exploring how the evolution of long bone shape in the group could be associated with these parameters. A few studies previously investigated shape variation of the limb bones in either modern or fossil rhinocerotoids, but rarely in regard to mass, size or degree of brachypody/gracility (Guérin, 1980;Prothero & Sereno, 1982;Becker, 2003;Malletet al., 2019,2020;Etienneet al., 2020). To date, no comprehensive morphofunctional analysis has explored covariation patterns between the shape of the long bones and each of these parameters at the scale of the entire superfamily.

Here, we investigate the shape variation of the forelimb bones among the superfamily Rhinocerotoidea in relation to bone size, body mass and degree of gracility. We performed phylogenetically-informed shape analyses of the three forelimb bones (humerus, radius, ulna) in a 3D geometric morphometric context. We chose to focus on forelimb bones, because they play a crucial role in supporting body weight and in braking during locomotion in quadrupeds (Hildebrand, 1974;Duttoet al., 2006;Henderson, 2006). Previous works on modern rhinos indicate a greater association of both mass and size with the shape of the forelimb bones over that of the hind limb ones (Malletet al., 2019,2020). In accordance with the literature, we hypothesize: (1) a strong association of bone size, body mass and degree of gracility with bone shape; (2) different expression of this association on the stylopodium and zeugopodium respectively (Alexanderet al., 1979;Prothero & Sereno, 1982;Biewener, 1989b;Bertram & Biewener, 1992;Malletet al., 2019,2020); and (3) a strong phylogenetic signal in bone shape variation, with differences depending on the considered bone. Testing these hypotheses will enable us to precisely highlight whether and how body mass could have played a role in shaping the evolution of forelimb bones among rhinocerotoids.

MATERIAL AND METHODS

Studied specimens

We selected 283 modern and fossil specimens housed at 15 institutions representing a total of 94 humeri, 105 radii and 84 ulnae (seeSupporting Information, Table S1, for the complete list of studied specimens). The data set included 69 taxa (five modern and 64 fossil species) belonging to almost all families of the superfamily Rhinocerotoidea (no representatives of the recently defined family Eggysodontidae were included) (Fig. 1). Taxa were selected to include as much body shape and mass diversity as possible and to cover the largest temporal range, however, selection also depended on what material was available. Taxonomic attributions were verified or updated using recent literature, directly with specimen numbers when available, or using taxonomic lists and institution databases for each locality. We retained the most recent binomial names considered as correct following the International Commission on Zoological Nomenclature rules (seeSupporting Information, Table S1).

Figure 1.

Composite cladogram of the studied species among Rhinocerotoidea (Mammalia: Perissodactyla). Families, subfamilies, tribes and subtribes are defined by a colour code following the cladistic framework ofAntoineet al. (2003) andBeckeret al. (2013). All silhouettes representing a member of each group are to scale (provided bywww.phylopic.org under a Creative Commons license).

Open in new tabDownload slide

We only considered adult individuals with fully fused epiphyses. We chose complete bones displaying no or negligible taphonomic effects (e.g., shallow surface cracks not altering the global shape), rejecting massively crushed specimens or those restored with plaster. We also considered incomplete bones in partial shape analyses (see below), as long as they were not crushed or distorted. Almost no information regarding sex was available for fossil specimens: even if sexual dimorphism is known for some species and may slightly affect the shape of long bones (Guérin, 1980;Dinerstein, 1991;Mead, 2000;Zschokke & Baur, 2002;Mihlbachler, 2007;Chenet al., 2010), we assumed that this intraspecific variation was largely exceeded by interspecific shape changes (according toMalletet al., 2019). For each species, we selected up to three specimens for each bone. All anatomical terms follow classic veterinary terminology and anatomical works on Perissodactyla and rhinoceroses (Guérin, 1980;Federative Committee on Anatomical Terminology, 1998;Antoine, 2002;Prothero, 2005;Barone, 2010a;Heissig, 2012;Baiet al., 2017). These terms are illustrated in theSupporting Information (Fig. S1). Locations of muscle insertion followEtienneet al. (2021).

3D models

Most bones were digitized with a structured-light three-dimensional scanner (Artec Eva) and reconstructed with Artec Studio Professional software (v.12.1.1.12;Artec 3D, 2018). This software was also used to reconstruct bones broken in two or more pieces (without any parts lacking) into a single complete mesh. Some specimens were digitized with a photogrammetric approach, followingMallison & Wings (2014) andFauet al. (2016). Sets of photos were used to reconstruct 3D models using Agisoft Photoscan software (v.1.4.2;Agisoft, 2018). Two specimens were digitized using medical computed tomography scanners at the Royal Veterinary College, London (Equine Hospital) and at the University of California, San Francisco (Department of Radiology & Biomedical Imaging). For these specimens, bone surfaces were extracted as meshes using Avizo software (v.9.5.0;Thermo Fisher Scientific, 2018). Because a few specimens displayed small parts lacking on the shaft, we used Geomagic Studio (v.2014.3.0.1781;3D Systems Corporation, 2014) to fill holes. We used the ‘curvature filling’ tool to ensure that the added polygons matched the curvature of the surrounding mesh. Finally, each mesh was decimated to reach 250 000 vertices and 500 000 faces using MeshLab software (v.2016.12;Cignoniet al., 2008). We performed our analyses on left bones: when left bones were not available in some specimens, we used mirrored right bones instead.

3D geometric morphometrics

The shape variation of our sample was analysed using a 3D geometric morphometrics approach, a widely-used methodology to quantify and visualize the morphological differences between objects by comparing the spatial coordinates of points called landmarks (Adamset al., 2004;Zelditchet al., 2012). We quantified bone shape by placing a set of anatomical landmarks and curve and surface sliding semi-landmarks on the meshes, followingGunz & Mitteroecker (2013) andBotton-Divetet al. (2016). Anatomical landmarks and curves were placed on meshes using the IDAV Landmark software (v.3.0;Wileyet al., 2005). We created a template to place surface sliding semi-landmarks for each bone. The geometric location of landmarks and sliding semi-landmarks is derived from previous morphometric works on rhinoceros long bones (Malletet al., 2019,2020) to cover the shape diversity of the sample (seeSupporting Information, Data S1, for details on landmark numbers and locations). Two specimens [Dicerorhinus sumatrensis (Fischer, 1814) NHMUK ZE 1948.12.20.1, for the humerus and the ulna, andCeratotherium simum (Burchell, 1817) RMCA 1985.32-M-0001, for the radius] were chosen to be the initial specimens on which all anatomical landmarks, curve and surface sliding semi-landmarks were placed. We selected these two individuals for their average shape and size ensuring that all points were correctly projected on other bones despite the great shape and size ranges of the sample. These specimens were then used as templates for the projection of surface sliding semi-landmarks on the surface of all other specimens. Projection was followed by a relaxation step to ensure that projected points matched the actual surface of the meshes. Curve and surface sliding semi-landmarks were then slid to minimize the bending energy of a thin plate spline (TPS) between each specimen and the template at first, and then four times between the result of the previous step and the Procrustes consensus of the complete data set. Therefore, all landmarks can be treated at the end as geometrically homologous (Gunzet al., 2005;Gunz & Mitteroecker, 2013).

Because we chose to work at the species level, we then computed and analysed species mean shapes (Botton-Divetet al., 2017;Serioet al., 2020). After the sliding step, we computed a first Generalized Procrustes Analysis (GPA) with all specimens to remove the effect of size, location and orientation of the different landmark conformations (Gower, 1975;Rohlf & Slice, 1990). Then we computed the Procrustes consensus (or mean shape) of each species in the same geometric space. These Procrustes consensuses were superimposed in a second GPA in order to pool all species means into a single morphospace. This process was repeated for each bone separately. Because our data set contained more variables than observations, we computed a Principal Component Analysis (PCA) to reduce dimensionality (Baylac & Frieß, 2005;Gunz & Mitteroecker, 2013) and visualize the distribution of the species in the morphospace. We also computed theoretical shapes associated with both the minimum and maximum of the first two components of PCAs using a TPS deformation of a template mesh. Phylogenetic relationships between taxa (see below) were then plotted in the morphospace and compared to Neighbour Joining (NJ) trees computed on PC scores. Projection, relaxation, sliding processes, GPAs, PCAs and theoretical shape computation were conducted using the ‘Morpho’ package (v.2.8;Schlager, 2017) in the R environment (v.3.5.3;R Core Team, 2014). Phylogeny was plotted on the morphospace using the ‘geomorph’ package (v.3.2.1;Adams & Otárola-Castillo, 2013). NJ trees were computed with the ‘ape’ package (v.5.3;Paradis & Schliep, 2019).

Analyses of partial bones

Fossil long bones of rhinoceros can show redundant breakage patterns due to various taphonomic agents throughout the diagenesis process, for example high sedimentary pressure on fragile anatomical areas, trampling by heavy animals after burial (Hullot & Antoine, 2020) or scavenger action on parts containing marrow (Guérin, 1980). This is notably the case of the proximal part of the humerus or the olecranon process of the ulna, which was frequently damaged and prevented us from using some specimens in whole bone shape analyses. In order to include a higher number of relevant taxa in our sample despite these alterations, we performed analyses of partial bones presenting important lacking parts. We also included complete bones in the analyses of partial bones. FollowingBarduaet al. (2019), we used curve-sliding semi-landmarks to define artificial lines acting as a limit for the sliding of surface semi-landmarks and virtually remove damaged or lacking parts from analyses. These limit lines involved at least one anatomical landmark to ensure that they were geometrically homologous on all specimens. They were also placed on complete bones, which were all included in the analyses of partial bones. Limit lines were finally removed after the sliding process and before the GPA to consider only true biological shape information in our analyses. Three data sets were used: distal half of the humerus, ulna without olecranon tubercle and proximal half of the ulna (seeSupporting Information, Data S1, for details on landmarks and sliding semi-landmarks in templates of partial bones).

Phylogenetic framework

Although recent publications refined the phylogenetic relationships in Rhinocerotoidea (Wanget al., 2016;Tissieret al., 2018) and in Ceratomorpha (Baiet al., 2020), these studies only include a small part of all genera of rhinocerotoids known worldwide. To date, no comprehensive and consensual phylogeny of the whole superfamily Rhinocerotoidea exists. To assess the effect of phylogenetic relationships on shape variation, we constructed a composite cladogram using trees previously computed on craniodental and postcranial characters or molecular data. Branch relations, lengths and occurrence dates were reconstructed after the works ofCerdeño (1995),Antoine (2002),Antoineet al. (2003,2010,2022),Prothero (2005),Boada-Saña (2008),Piraset al. (2010),Beckeret al. (2013),Lu (2013),Wanget al. (2016),Averianovet al. (2017),Tissieret al. (2018,2020) andBaiet al. (2020). We used the cladistic framework ofAntoineet al. (2003) andBeckeret al. (2013) to define families, subfamilies, tribes and subtribes (Fig. 1). The relationships between the five modern taxa remain controversial, especially regarding the position of the Sumatran rhinoceros (Dicerorhinus sumatrensis) and its extinct relatives (e.g.Tougardet al., 2001;Orlandoet al., 2003;Fernandoet al., 2006;Price & Bininda-Emonds, 2009;Steiner & Ryder, 2011;Yuanet al., 2014;Welkeret al., 2017;Cappelliniet al., 2019;Margaryanet al., 2020). It is likely that these uncertainties may be due to a hard polytomy at the base of the crown-group containing the five modern species (Willerslevet al., 2009;Gaudry, 2017). We therefore chose to consider a hard polytomy in our analyses and to address phylogenetic uncertainties using an Nearest Neighbour Interchange (NNI) procedure (see below).

To address the effect of phylogenetic relationships on shape data for each bone, we evaluated their phylogenetic signal by computing a multivariate K statistic (K_mult) on PC scores (Adams, 2014). This index allows the comparison between the rate of observed morphological change and that expected under a Brownian motion on a given phylogeny (Blomberget al., 2003;Adams, 2014). Because the K_mult computation requires fully bifurcating trees, we removed polytomies using the functionmulti2di in the ‘ape’ package (Paradis & Schliep, 2019). This function resolves polytomies by randomly creating a new branch with a null length from one branch of the polytomous node (Swenson, 2014;Paradis & Schliep, 2019). K_mult was then computed using the functionK.mult in the ‘phylocurve’ package (Goolsby, 2015).

Body mass, centroid size and gracility index

We explored the association of three variables related to body proportions and size (body mass, centroid size of the bone and gracility index) with the shape of each long bone of the forelimb within Rhinocerotoidea. Mean body mass (BM) of each species was retrieved from the literature, compiling up to three estimations per species to compute mean BMs (seeTable 1;Supporting Information, Table S2). However, BM estimations are highly heterogenous and can vary by a factor of three for a single species depending on the considered method and morphological proxy (dental, cranial or postcranial measurements), the specimen developmental stage and the geological formation. Moreover, regression equations for BM estimation were rarely developed for Perissodactyla or rhinoceroses only, resulting in potentially biased results for fossil Rhinocerotoidea (Prothero & Sereno, 1982). We managed to collect BM estimation for only 40 of the 69 taxa constituting our sample. Consequently, we chose to also consider the centroid size (CS) of each bone, which is classically used to address allometric variation, i.e. shape variation linked to size (Zelditchet al., 2012;Mitteroeckeret al., 2013;Klingenberg, 2016;Hallgrímssonet al., 2019). Centroid size, defined as the square root of the sum of the square of the distance of each point to the centroid of the landmark set (Zelditchet al., 2012), is known to be a good proxy of the mass of the animal (Ercoli & Prevosti, 2011;Cassiniet al., 2012), especially for limb bones of rhinoceros (Malletet al., 2019;Etienneet al., 2020). Given the large range of body shapes in Rhinocerotoidea (Fig. 1) and the fact that the same mass can be associated with both a slender or a robust body condition, we used the mean gracility index (GI-MC3) as an estimator of the degree of brachypody (seeTable 1;Supporting Information, Table S2). This index is computed dividing the transverse width of the third metacarpal by its maximal length and has been much used for rhinocerotoids (Colbert, 1938;Arambourg, 1959;Guérin, 1980;Cerdeño, 1998;Becker, 2003;Beckeret al., 2009;Scherleret al., 2013). The higher the GI-MC3 value, the shorter the limb length: species with a high GI-MC3 value are considered as more brachypodial (or less gracile) than species with low values. We computed this index by measuring third metacarpals, when available in collections, or compiling up to three GI-MC3 values in the literature to compute mean GI-MC3. These metacarpals were mostly associated with long bones for modern species, and mostly associated with a similar locality for fossil species (Supporting Information, Table S2). We addressed the effect of phylogeny on log-transformed CS, log-transformed cubic root of the mean BM and log-transformed mean GI-MC3 using the univariate K statistic to quantify the strength of the phylogenetic signal (Blomberget al., 2003). We tested for correlation between these three variables respectively using a linear regression on Phylogenetic Independent Contrasts (Felsenstein, 1985). We used the functioncontMap of the ‘phytools’ package (Revell, 2012) to plot these three variables along the phylogeny.

Variation patterns, and thus co-variation, can be expressed and analysed at different levels: across species (evolutionary variation), within a species at a single developmental stage (static variation) or within a species across developmental stages (ontogenetic variation) (Klingenberg, 2014). Here we explored the evolutionary co-variation of bone shape with each of the three variables (BM, CS, GI-MC3) considering a multivariate approach using Phylogenetic Generalized Least Squares (PGLS), a regression model taking into account the phylogenetic framework and computed here on Procrustes coordinates to quantify the shape variation related to CS, BM and GI-MC3 (Martins & Hansen, 1997;Rohlf, 2001;Klingenberg & Marugán-Lobón, 2013;Adams & Collyer, 2018). This was done using the functionprocD.pgls of the ‘geomorph’ package (v.3.2.1;Adams & Otárola-Castillo, 2013), suited for 3D geometric morphometric data. However, the functionprocD.pgls uses a Brownian Motion of evolution to compute PGLS, which assumes non-directional trait changes, while other models might assume a different computational hypothesis. To account for these changes depending on the considered model, we also computed PGLS under a Phylogenetic Ridge Regression (RR) model of evolution (Castiglioneet al., 2018). Phylogenetic RR takes into account variations of evolutionary rates along the different branches of a phylogenetic tree, accounting for potential accelerations and decelerations of the phenotypic changes among groups in a more accurate way than Brownian Motion. Therefore, we used the functionPGLS_fossil of the ‘RRphylo’ package (v.2.5.0;Castiglioneet al., 2018) to compute PGLS with a Ridge Regression model and compare it to the results obtained under a Brownian Motion model in order to see whether our results were robust to model variations.

Because the phylogeny of Rhinocerotoidea remains debated for both extant and extinct taxa (see above), we assessed the effect of potential uncertainty in taxa position in the phylogeny on PGLS by using a NNI procedure. The NNI algorithm generates new trees by swapping two adjacent branches of a specified tree (Felsenstein, 2004). We generated new trees using thenni function of the package ‘phangorn’ (Schliep, 2011) and computed PGLS with these rearranged trees to estimate the ranges of R² andP values.

All statistic tests have been considered as significant forP values ≤ 0.01. However, given that recent statistical works call for a continuous approach of theP value (Hoet al., 2019;Wassersteinet al., 2019), we chose to mention results having aP value up to 0.05 as well.

RESULTS

Correlation between BM and GI-MC3

The evolutionary variations of mean BM and mean GI-MC3 both show a significant phylogenetic signal (K_BM = 1.75,P < 0.001; K_GI-MC3 = 1.70,P < 0.001) and they are significantly correlated with one another when phylogeny is taken into account (correlation coefficient [r] = 0.44,P < 0.001). The mapping of mean BM and GI-MC3 along the phylogeny (Fig. 2) clearly indicates that, despite this significant correlation, there is not a strict correspondence between high BM and high GI-MC3 values. This is particularly visible for Paraceratheriidae, large Elasmotheriinae and Teleoceratina.

Figure 2.

Evolution of body mass (BM) and gracility index (GI-MC3) along the phylogeny for the studied species in Rhinocerotoidea (Mammalia: Perissodactyla). Left: mean BM; Right: mean GI-MC3. Computations were made on log-transformed cubic root of mean BM and log-transformed GI-MC3. Values at nodes and along branches were reconstructed based on a Brownian motion model of evolution (Revell, 2012). Colour coding for taxa followsFigure 1. Dashed lines indicate missing data. Evolution of the third metacarpal shape depending on the GI-MC3 value is illustrated by specimensHyrachyus modestus AMNH FM 17436 (minimum) andTeleoceras fossiger AMNH FM 2636 (maximum).

Open in new tabDownload slide

Differences in PGLS between Brownian Motion and phylogenetic RR

PGLS computed under a Brownian Motion model (using the geomorph functions) and under a Phylogenetic RR model (using the RRphylo functions) show similar results (seeSupporting Information, Table S3, for a detailed comparison between both models). Significant regressions under a Brownian Motion model remain significant under a RR model, in addition, non-significant results under a Brownian Motion model remain non-significant under a RR model. R²,P-values and shape deformations are extremely close in both cases. Only regression plots differ, those obtained under a RR model showing a much higher spread of specimens, making their interpretation more difficult. For all these reasons, we chose to present only results obtained under a Brownian Motion model in the following sections.

Humerus – complete bone

The species distributions in the NJ tree (Fig. 3A) and in the phylomorphospace (Fig. 4A) computed on the complete humeri are mostly congruent with phylogeny, which is not surprising since the phylogenetic signal carried by its shape variation is strong (K_mult = 1.16,P < 0.01). Along the NJ tree, small-sized and early-diverging Hyrachyidae and Hyracodontidae are followed by a cluster mixing Rhinocerotidae and Rhinocerotinaeincertae sedis withParamynodon Matthew, 1829 (Amynodontidae),Urtinotherium Chow & Chiu, 1963 (Paraceratheriidae),Menoceras Troxell, 1921 (Elasmotheriinae), some Aceratheriini and evenDicerorhinus sumatrensis. Other Aceratheriini are grouped close to Teleoceratina andDihoplus megarhinus de Christol, 1835 (Rhinocerotini), while almost all Rhinocerotina form a well-separated group (Fig. 3A). The phylomorphospace of the first two axes of the PCA, representing 63.7% of the global variance, is structured in a similar way (Fig. 4A). PC1 carries 54% of the variance. Along PC1, Hyrachyidae and Hyracodontidae plot towards negative values whileParamynodon is close to a central cluster groupingUrtinotherium,Menoceras, Aceratheriini and Teleoceratina, as well as Rhinocerotidae and Rhinocerotinaeincertae sedis. In this cluster,Aphelops Cope, 1874 shares a shape proximity with all Teleoceratina, whereas other Aceratheriini are closer to more ancient taxa (Amphicaenopus Wood, 1907,Menoceras,Plesiaceratherium Young, 1937,Protaceratherium Forster-Cooper, 1911 andTrigonias Lucas, 1900). All members of the subtribe Rhinocerotina group together towards positive values withDicerorhinus Gloger, 1841,Dihoplus Brandt, 1878,Rhinoceros Linnaeus, 1758 andStephanorhinus Kretzoi, 1942 overlapping theAphelops–Teleoceratina cluster. The highest PC1 values are associated with the modern African clade (Ceratotherium Gray, 1868–Diceros Gray, 1821) and their extinct relatives, and theCoelodonta Bronn, 1831 clade. PC2 represents 9.7% of the global variance. It is mainly driven by an opposition between Hyrachyidae, Hyracodontidae and Rhinocerotina towards negative values and Amynodontidae, Paraceratheriidae and all other Rhinocerotidae towards positive values.Urtinotherium is strongly isolated from all other species towards maximal positive values.

Figure 3.

Neighbour Joining trees computed on all PC scores obtained from the PCAs performed on shape data. Colour coding followsFigure 1 and abbreviations followTable 1. Point size is proportional to the mean log centroid size of each species. A, complete humerus; B, distal partial humerus; C, radius; D, complete ulna; E, ulna without olecranon tuberosity; F, proximal partial ulna.

Open in new tabDownload slide

Figure 4.

Results of the PCA performed on morphometric data of complete humerus (A) and distal partial humerus (B) and shape variation associated with the first two axes of the PCA (caudal view). Blue: negative side of the axis. Orange: positive side of the axis. Phylogenetic relationships are plotted in the morphospace. Colour coding followsFigure 1 and abbreviations followTable 1. Point size is proportional to the centroid size of each species.

Open in new tabDownload slide

The shape variation along PC1 is mainly related to bone slenderness (Fig. 3A;Supporting Information, Fig. S2A). Towards the minimal values, the humerus is thin and straight, with a greater trochanter developed craniomedially; an asymmetrical bicipital groove; a rounded humeral head oriented proximocaudally; a poorly developed deltoid tuberosity; a poorly developed supracondylar crest; a narrow olecranon fossa; and a symmetrical trochlea with a developed capitulum. The shape associated with maximal values is highly robust and thick, with a strong development of the lesser tubercle over the greater one; a large symmetrical bicipital groove with an intermediate tubercle; a deltoid tuberosity highly developed laterally; a strong development of the lateral epicondyle and the epicondylar crest; a large and rectangular olecranon fossa; and an asymmetrical trochlea with a reduced capitulum. Along PC2, shape variation mainly concerns epiphyseal elements. Towards the positive maximum, the humerus displays a greater tubercle developed cranially; a rounded head oriented proximally; a strong deltoid tuberosity situated at the middle of the shaft; a larger shaft diameter; a strong proximolateral development of the epicondylar crest; and a trochlea flattened proximodistally. The shape associated with minimal values exhibits a deltoid tuberosity situated above the midshaft; a poorly developed epicondylar crest with a lateral epicondyle directed laterodistally; and an asymmetrical trochlea medially developed.

The evolutionary variation of the centroid size of complete humeri bears a significant phylogenetic signal (K_CS = 1.28,P < 0.001) and is highly correlated with BM (r = 0.64) and marginally correlated with GI-MC3 (r = 0.37,P = 0.03) (Table 2). PGLS results indicate that CS, BM and GI-MC3 are all significantly correlated with humerus shape (Table 3). The NNI procedure indicates that the correlation with BM is more strongly affected by phylogenetic uncertainties than that with CS (Table 3). This may be related to a smaller and less diverse sample for BM values. Regression of shape against CS shows a good fit to the regression line, most of the species following a marked common trend with little divergences far away from the line (Fig. 5A). Most of the Rhinocerotini (i.e., Rhinocerotina and Teleoceratina) are situated below the regression line, while the other species are situated above.Urtinotherium appears as slightly shifted from the general trend. In the absence of taxa such as Hyracodontidae, Amynodontidae and Paraceratheriidae, regression of shape against BM shows a good fit to the regression line. The trend is strongly driven byHyrachyus Leidy, 1871, which potentially constitutes a bias. However, a clear separation exists between Aceratheriini, all situated below the regression line, and Rhinocerotini, mainly situated above the line (Fig. 5B). Results for GI-MC3 indicate a good fit to the regression line as well. Almost all Rhinocerotina group together above the line whereas Teleoceratina are situated below the line. All other species are mixed close to the common trend. Hyrachyidae and Hyracodontidae are isolated towards the minimal values (Fig. 5C). If shape variation related to these three variables mainly concerns an increase of robustness towards maximal values (Fig. 5;Supporting Information, Fig. S3A–C), that related to BM (that lacks heavy Paraceratheriidae) is slightly different from those related to CS and GI-MC3, with a stronger mediolateral development of both epiphyses relatively to the shaft (Fig. 5B). Most of the shape variation occurs on the medial face of the bone and on strong muscular insertions like the deltoid tuberosity and the epicondylar crest for the three variables. In addition, BM variation affects the bicipital groove, while variation of GI-MC3 implies shape changes located distally and caudally to the humeral head, from the deltoid tuberosity and tricipital line to the lesser tubercle convexity (Fig. 5C).

Figure 5.

Significant PGLS regression plots for complete humerus performed on shape data and log-transformed centroid size (CS) (A), log-transformed cubic root of mean body mass (BM) (B), log-transformed mean gracility index (GI-MC3) (C). Point colour coding followsFigure 1. Point size is proportional to mean log CS of each species. On the right, shapes associated with minimum and maximum fitted values (top row) and colour maps of the location and intensity of the shape deformation (bottom row). Blue: minimum value of the regression. Orange: maximum value of the regression. For each bone, the shape associated with the minimum was coloured depending on its distance to the shape associated with the maximum (blue indicates a low deformation intensity and red indicates a high deformation intensity). Orientation from left to right in each case: caudal, lateral, cranial and medial.

Open in new tabDownload slide

Humerus – distal part

The phylogenetic signal carried by the shape variation of the distal humeri is strong (K_mult = 1.22,P < 0.01). The species distributions in the NJ tree (Fig. 3B) and in the phylomorphospace are highly similar to those observed for the complete humeri (Fig. 4B). On the NJ tree, all Amynodontidae are grouped together withJuxia Chow & Chiu, 1964 (small Paraceratheriidae), while the giant Paraceratheriidae group together are close to some Aceratheriini (Aphelops,Chilotherium Ringström, 1924). Other Aceratheriini are mixed with Teleoceratina and more basal taxa, while Rhinocerotina form a homogeneous cluster all together. A similar organization is observable in the phylomorphospace, where the first two axes represent 70.2% of the global variance. PC1 carries 55% of the global variance and PC2 carries 15.2%. The species distribution along both axes is largely similar to that observed for the complete humerus (Fig. 4B). Small and large Amynodontidae group together with the light paraceratheriidJuxia, while heavier Paraceratheriidae form an isolated cluster along PC2. Within Rhinocerotina, species seem distributed from the smallest to the largest along PC1 despite some exceptions (e.g.,Dihoplus megarhinus,Rhinoceros unicornis Linnaeus, 1758).Chilotherium shows the highest positive value on PC1.

The shape variation along PC1 is similar to that observed on complete bones (Fig. 4B;Supporting Information, Fig. S2B). Towards positive maximal values, PC1 is mainly associated with an increase of thickness, with a strong development of the epicondylar crest; a broad olecranon fossa; and an asymmetrical trochlea with a reduced capitulum. Along PC2, the shape variation is also almost identical to that observed on complete bones.

The evolutionary variation of the centroid size of partial humeri carries a significant phylogenetic signal (K_CS = 1.39,P < 0.001). The correlation between CS and BM is higher than for the complete humeri (r = 0.72) and correlation between CS and GI-MC3 is significant (r = 0.50) (Table 2). As for complete bones, PGLS results indicate a significant correlation between humerus shape and CS, and BM and GI-MC3, respectively. The NNI procedure indicates that phylogenetic uncertainties do not highly affect the relation between shape and the three variables (Table 3). The regression plot of shape against CS indicates an excellent fit to the regression line with a tendency similar to that observed on complete bones, but with Amynodontidae, Hyrachyidae, Hyracodontidae and Paraceratheriidae slightly shifted towards less robust shapes for a given CS than Rhinocerotidae, which have a generalized ‘large-and-heavy-head’ plan, with respect to other families among Rhinocerotoidea. (Fig. 6A). The presence of Amynodontidae and Paraceratheriidae in the regression of shape against BM highlights a similar tendency and a strong fit to the regression line (Fig. 6B). The regression plot of shape against GI-MC3 is almost identical to that obtained on complete bones with a good fit to the regression line as well (Fig. 6C). Similarly, shape variation is like that of complete bones for the three variables, mainly affecting the general robustness and muscular insertions such as the epicondylar crest that is broadened (Fig. 6;Supporting Information, Fig. S3D–F). Only shape variation associated with BM slightly differs with an epicondylar crest less developed than for complete bones towards maximum values.

Figure 6.

Significant PGLS regression plots for distal partial humerus performed on shape data and log-transformed centroid size (CS) (A), log-transformed cubic root of mean body mass (BM) (B), log-transformed mean gracility index (GI-MC3) (C). Point colour coding followsFigure 1. Point size is proportional to mean log CS of each species. On the right, shapes associated with minimum and maximum fitted values (top row) and colour maps of the location and intensity of the shape deformation (bottom row). Blue: minimum value of the regression. Orange: maximum value of the regression. For each bone, the shape associated with the maximum was coloured depending on its distance to the shape associated with the minimum (blue indicates a low deformation intensity and red indicates a high deformation intensity). Orientation from left to right in each case: caudal, lateral, cranial and medial.

Open in new tabDownload slide

Radius

As for the humerus, the phylogenetic signal carried by shape data of the radii is strong (K_mult = 1.15,P < 0.01). However, the species distributions in the NJ tree (Fig. 3C) and in the phylomorphospace (Fig. 7) are less reminiscent of the phylogeny and seem likely related to the degree of brachypody. Along the NJ tree, Hyrachyidae group with Hyracodontidae, Paraceratheriidae and small Elasmotheriinae. Aceratheriini, Rhinocerotina and Teleoceratina are mixed together with larger Elasmotheriinae, most of the species being sorted by their gracility rather than mass or size. This pattern is highly similar to that seen on the PCA, with the first two axes representing 75% of the global variance (Fig. 7). PC1 gathers 70.7% of the global variance. Along this axis,Triplopus Cope, 1880 (Hyracodontidae) constitutes the positive maximum. Contrary to the morphospace obtained for the humerus, two of the biggest species of the sample,Juxia andUrtinotherium (Paraceratheriidae), plot together with the smallest and lightest species.Paraceratherium groups with small Elasmotheriinae and Rhinocerotidaeincertae sedis, as well asAmynodon Marsh, 1877 andParamynodon (Amynodontidae). Towards negative values, Aceratheriini and Rhinocerotini (Teleoceratina and Rhinocerotina) are grouped together with larger Elasmotheriinae (Hispanotherium Crusafont & Villalte, 1947 andElasmotherium Fischer, 1808). Within this cluster,Stephanorhinus,Dicerorhinus and someDihoplus plot withAphelops,Peraceras Cope, 1880 andHoploaceratherium Ginsburg & Heissig, 1989, whereas larger Rhinocerotina (Ceratotherium,Rhinoceros,Diceros,Coelodonta) are closer toBrachypotherium Roger, 1904 andDiaceratherium Dietrich, 1931 (Teleoceratina). OnlyTeleoceras Hatcher, 1894 andCoelodonta antiquitatis (Bronn, 1831) plot outside the main cluster towards the maximal negative values. PC2 represents only 4.3% of the variance and no obvious organization of the specimens is visible along this axis.

Figure 7.

Results of the PCA performed on morphometric data of the radius and shape variation associated with the first axis of the PCA (cranial view). Blue: negative side of the axis. Orange: positive side of the axis. Phylogenetic relationships are plotted in the morphospace. Colour coding followsFigure 1 and abbreviations followTable 1. Point size is proportional to the mean log centroid size of each species.

Open in new tabDownload slide

As for the humerus, shape variation of the radius along PC1 is mainly related to bone slenderness (Fig. 7;Supporting Information, Fig. S2C). The shape associated with maximal values is thin and slender, with slight craniocaudal and mediolateral bends; a rectangular glenoid cavity with a lateral expansion for the capitulum; a shaft mediolaterally as large as the two epiphyses; a rectangular and shallow distal articular surface; and a poorly developed radial styloid process. Conversely, the shape associated with the minimal values is massive with a large asymmetrical glenoid cavity; almost no lateral development of the cavity for the capitulum; both epiphyses mediolaterally larger than the diaphysis; a radial styloid process developed distally; and a rectangular and deep distal articular surface.

As for the humerus, the evolutionary variation of the centroid size of the radius carries a significant phylogenetic signal (K_CS = 0.82,P < 0.001). The correlation between CS and BM is significant and high (r = 0.80), whereas CS and GI-MC3 are not correlated (Table 2). PGLS results indicate that the radius shape is significantly correlated with BM and GI-MC3, the latter correlation being stronger than the former (Table 3). PGLS computed on NNI trees indicate that correlation with BM is affected by phylogenetic uncertainties and may be non-significant depending on the tree configuration. Conversely, correlation with CS appears always as non-significant and GI-MC3 always as significant for whatever the tree configuration is (Table 3). The regression plot of radius shape against BM indicates a rather good fit to the regression line (Fig. 8A). Paraceratheriidae andTeleoceras deviate strongly from the common regression trend, whereas Rhinocerotina and Elasmotheriinae strongly follow it. The shape variation associated with maximal values of BM is mainly related to a mediolateral development of both epiphyses, notably on the lateral part of the proximal epiphysis, where the m. biceps brachii inserts (Etienneet al., 2021). These changes are also associated with a slight increase of robustness towards high BM values (Fig. 8A;Supporting Information, Fig. S3G). The regression plot of shape against GI-MC3 indicates an excellent fit to the regression line, with a strong common trend shared by all members of the superfamily. Although most Rhinocerotina are situated above the regression line, they are mixed together with Aceratheriini, Teleoceratina and large Elasmotheriinae. Giant Paraceratheriidae plot together with small Elasmotheriinae and almost all Amynodontidae, while the small paraceratheriidJuxia is close toHyrachyus andHyracodon Leidy, 1856.Triplopus plots towards minimal values (Fig. 8B). GI-MC3 variation is correlated with a mediolateral development of the bone appearing stronger on the lateral side of both epiphyses than on the medial one, although less marked than for BM (Fig. 8B;Supporting Information, Fig. S3H) and with an overall increase in robustness.

Figure 8.

Significant PGLS regression plots for the radius performed on shape data and log-transformed cubic root of mean body mass (BM) (A) and log-transformed mean gracility index (GI-MC3) (B). Point colour coding followsFigure 1. Point size is proportional to mean log CS of each species. On the right, shapes associated with minimum and maximum fitted values (top row) and colour maps of the location and intensity of the shape deformation (bottom row). Blue: minimum value of the regression. Orange: maximum value of the regression. For each bone, the shape associated with the minimum was coloured depending on its distance to the shape associated with the maximum (blue indicates a low deformation intensity and red indicates a high deformation intensity). Orientation from left to right in each case: caudal, lateral, cranial and medial.

Open in new tabDownload slide

Ulna – complete bone

The shape variation of the complete ulnae carries a strong phylogenetic signal (K_mult = 0.93,P < 0.01). The NJ tree (Fig. 3D) shows a grouping of Hyrachyidae, Hyracodontidae, small Elasmotheriinae andJuxia (which slightly isolates from this cluster). Aceratheriini group together with someDiaceratherium but also Amynodontidae, while all other Teleoceratina are grouped together and slightly isolate from other species.Metamynodon Scott & Osborn, 1887 (Amynodontidae) is placed between Aceratheriini and Teleoceratina, while all Rhinocerotina group together (also withDiaceratherium lamilloquense Michel, 1983). A similar structure is observed on the phylomorphospace, with the first two axes representing 70.2% of the global variance (Fig. 9A). The first axis carries 54.7% of the variance.Juxia (Paraceratheriidae) plots towards minimal values. Small Elasmotheriinae group together withTrigonias,Protaceratherium (Rhinocerotidaeincertae sedis) andParamynodon (Amynodontidae) towards minimal values. However,Amphicaenopus andMetamynodon group with a cluster containing Aceratheriini and Rhinocerotini, as well as someDiaceratherium. Within Rhinocerotina, larger taxa such asCeratotherium simum,Rhinoceros unicornis,Coelodonta antiquitatis orDihoplus pikermiensis Toula, 1906 group towards slightly higher values.Prosantorhinus,Brachypotherium,Diaceratherium aurelianense Nouel, 1866 andTeleoceras constitute the highest positive values. The second axis accounts for 15.5% of the global variance.Hyrachyus and Rhinocerotina group together in the negative part of the axis with almost no overlapping of the other species. Rhinocerotidaeincertae sedis plot around null values together withHyracodon, Amynodontidae, small Elasmotheriinae,Hoploaceratherium andDiaceratherium lamilloquense. All other Teleoceratina group with Aceratheriini andJuxia towards the highest positive values.

Figure 9.

Results of the PCA performed on morphometric data of complete ulna (A), ulna without olecranon tuberosity (B) and proximal partial ulna (C) and shape variation associated with the first two axes of the PCA (caudal view). Blue: negative side of the axis. Orange: positive side of the axis. Phylogenetic relationships are plotted in the morphospace. Colour coding followsFigure 1 and abbreviations followTable 1. Point size is proportional to the mean log centroid size of each species.

Open in new tabDownload slide

As for the humerus and the radius, the shape variation of the ulna along PC1 is mainly related to the bone slenderness (Fig. 9A;Supporting Information, Fig. S2D). The shape associated with minimal values is very thin and slender with an olecranon tuberosity developed proximally; a symmetrical and mediolaterally flattened articular surface for the humerus; a shaft bent in the craniocaudal direction and highly compressed mediolaterally; a narrow and shallow distal articular surface; and an articular surface for the pisiform developed proximally. Conversely, the shape associated with maximal values is robust and massive with a strong olecranon tuberosity developed proximocaudally; a large and asymmetrical articular surface for the humerus; a massive and straight shaft with a triangular section; a distal epiphysis developed mediolaterally; a wide and deep distal articular surface; and a reduced articular surface for the pisiform. Along PC2, the shape associated with minimal values displays an olecranon tuberosity developed proximodistally; an anconeus process developed cranially; a shaft bent craniocaudally; and a narrow distal articular surface. The shape associated with maximal values displays an olecranon tuberosity developed mainly caudally; an anconeus process poorly developed cranially; a shaft with a curved caudal border and a straight cranial border; and a wide and medially tilted distal articular surface.

The evolutionary variation of the centroid size of the complete ulna carries a significant phylogenetic signal (K_CS = 0.84,P = 0.002). Centroid size is not significantly correlated with GI-MC3 but is marginally significantly correlated with BM (r = 0.44,P = 0.04) (Table 2). PGLS results highlight a strong and significant correlation between ulna shape and GI-MC3 only (Table 3). PGLS computed on NNI trees confirm that neither BM nor CS are significantly correlated with shape whatever the phylogenetic configuration (Table 3). As for the radius, the regression plot of shape against GI-MC3 shows a good fit to the regression line and highlights a strong common trend with few outliers. However, groups are more clearly separated than for the radius, with almost all Rhinocerotina plotting above the regression line, while other species plot below the line. Teleoceratina and Aceratheriini form well-separated groups with few overlapping with other species. Small Elasmotheriinae plot with Amynodontidae, whileHyrachyus (Hyrachyidae) andJuxia (Paraceratheriidae) plot towards minimal values (Fig. 10A). A higher GI-MC3 is associated with a more robust and straighter ulna, showing a craniocaudal and mediolateral broadening and a strong development of the olecranon tubercle, as well as a development of the lateral insertion area for digit extensors along the shaft (Fig. 10A;Supporting Information, Fig. S3I).

Figure 10.

Significant PGLS regression plots for complete ulna (A) and ulna without olecranon tuberosity (B) performed on shape data and log-transformed mean gracility index (GI-MC3). Point colour coding followsFigure 1. Point size is proportional to mean log CS of each species. On the right, shapes associated with minimum and maximum fitted values (top row) and colour maps of the location and intensity of the shape deformation (bottom row). Blue: minimum value of the regression. Orange: maximum value of the regression. For each bone, the shape associated with the minimum was coloured depending on its distance to the shape associated with the maximum (blue indicates a low deformation intensity and red indicates a high deformation intensity). Orientation from left to right in each case: caudal, lateral, cranial and medial.

Open in new tabDownload slide

Ulna – without the olecranon tuberosity

Shape data of the ulna without the olecranon tuberosity carry a strong phylogenetic signal (K_mult = 0.81,P < 0.01). The NJ tree (Fig. 3E) and phylomorphospace (Fig. 9B) are similar to those obtained for the complete ulnae. One of the main differences with the complete ulna is the position of the heavyElasmotherium (Elasmotheriinae): the NJ tree highlights that this genus shares shape similarity with poorly related taxa likeAmphicaenopus (Rhinocerotidaeincertae sedis) andMetamynodon (Amynodontidae). In the phylomorphospace, the two first axes of the PCA account for 68.9% of the global variance. PC1 represents 51.5% while PC2 accounts for 17.4%. Again,Elasmotherium plots far away from smaller Elasmotheriinae likeSubhyracodon Brandt, 1878 andMenoceras along PC1, and closer toAmphicaenopus,Aphelops andMetamynodon (Fig. 9B). The shape variation associated with both axes is largely equivalent to that observed for the complete ulna (Fig. 9B;Supporting Information, Fig. S2E). PC1 is mainly driven by a change of slenderness and proportion of both epiphyses relatively to the shaft, with a highly massive and robust bone towards the positive maximum. PC2 is mainly driven by changes of both orientation of the olecranon development and straightness of the shaft. Towards minimal values, the olecranon is oriented almost completely caudally and the cranial border of the shaft is fully straight.

As for the complete ulna, the evolutionary variation of the CS of the ulna without the olecranon tuberosity carries a significant phylogenetic signal (K_CS = 0.78,P = 0.003). Conversely, CS is significantly and strongly correlated with BM (r = 0.52), but not with GI-MC3 (Table 2). Results of the PGLS indicate only a significant correlation between shape and GI-MC3, which is not affected by phylogenetic uncertainties. Conversely, the correlation between shape and CS remains non-significant regardless of phylogenetic uncertainties (Table 3). If the regression plot displays a trend relatively similar to that observed on complete ulnae, the fit to the regression line is poorer. Rhinocerotini (Rhinocerotina and Teleoceratina) are much more distant from the common regression slope, contrary to what it is observed for the radius and complete ulna.Elasmotherium andAmphicaenopus plot close to Rhinocerotina, which form a well-isolated cluster above the regression line (Fig. 10B). Shape variation related to GI-MC3 is highly similar to that observed along PC1, with a much more pronounced variation along the lateral side of the shaft (Fig. 10B;Supporting Information, Fig. S3J). As observed on the radius, PGLS computed on BM display marginally non-significant results and NNI trees lead to significant or non-significant correlations between shape and BM depending on the considered phylogeny (Table 3). The regression plot shows a rather good fit to the regression line, despite some clear outliers. Most Rhinocerotina plot below the regression line, together with some Teleoceratina, while Aceratheriini form a central cluster.Elasmotherium plots towards maximal values whileMenoceras plots towards negative values. This poorly significant regression can be related to the isolation ofJuxia away from the common trend (seeSupporting Information, Fig. S4, for regression plot). The shape variation related to BM mainly concerns the caudal border of the ulna, particularly the area placed distally to the olecranon (seeSupporting Information, Fig. S4).

Ulna – proximal part

Shape data of the proximal parts of the ulnae carry a strong phylogenetic signal (K_mult = 0.72,P < 0.01). The NJ tree (Fig. 3F) and the phylomorphospace (Fig. 9C) show marked differences with previous analyses on the complete bones or on the ulna without the olecranon tubercle. The NJ tree is more congruent with phylogenetic groupings than is the phylomorphospace (Figs 3F,9C). Rhinocerotina form a homogeneous cluster (except forLartetotherium Ginsburg, 1974) close to a group containing small Elasmotheriinae,Protaceratherium,Hyracodon andHyrachyus. Aceratheriini and Teleoceratina are mixed together. Paraceratheriidae and Amynodontidae plot withAmphicaenopus (Rhinocerotidaeincertae sedis) among the Aceratheriini–Teleoceratina group. On the phylomorphospace, the two first axes of the PCA carry 55.2% of the global variance (Fig. 9C). PC1 represent 31.1% of the global variance. Along this axis,Hyrachyus is isolated towards positive values.Hyracodon (Hyracodontidae) and Amynodontidae plot in a cluster grouping Rhinocerotina, Elasmotheriinae and Rhinocerotidaeincertae sedis. Aceratheriini and Teleoceratina isolate towards negative values. Paraceratheriidae are placed between the Rhinocerotina cluster and the Aceratheriini–Teleoceratina one, together with other taxa likeAmphicaenopus,Lartetotherium andMetamynodon. The second axis, representing 24.1% of the variance, is mainly driven by the isolation of Paraceratheriidae from all other species, especially the two big forms of the genusParaceratherium, towards minimal values. Almost all other species form a single and mixed cluster from null to positive values without any clear organization.

As for the complete ulna, the shape variation of the proximal part of the ulna along PC1 mainly relates to the slenderness of the bone (Fig. 9C;Supporting Information, Fig. S2F). The shape associated with maximal values is thin and slender with a high olecranon tuberosity, developed in proximal direction and mediolaterally flattened; an anconeus process developed cranially; a symmetrical articular surface for the humerus flattened mediolaterally; and a long synostosis surface for the radius. Conversely, the shape associated with minimal values is thick and massive, with a strong olecranon tubercle developed proximocaudally and enlarged mediolaterally; an anconeus process poorly developed; a wide and asymmetrical articular surface for the humerus; and a short synostosis surface for the radius. Along PC2, variation is mainly driven by the proportion, shape and orientation of the olecranon. Towards minimal values, the proximal part of the ulna displays a massive and short olecranon, mediolaterally compressed and poorly caudally developed; a large and trapezoid articular surface for the humerus; a poorly-developed anconeus process; and a long synostosis surface for the radius developed medially. The shape associated with maximal values displays a thinner and squared olecranon developed proximocaudally; a more triangular articular surface for the humerus; and a short synostosis surface for the radius.

The evolutionary variation of the centroid size of the proximal part of the ulna carries a significant phylogenetic signal (K_CS = 1.91,P < 0.001). As for the ulna without the olecranon tuberosity, the centroid size is significantly and strongly correlated with BM (r = 0.85) but not with GI-MC3 (Table 2). PGLS indicate a significant correlation between shape and each of the three variables, similar to what is observed on the complete and partial humerus (Table 3). However, both regression plots of shape against CS and BM must be considered with caution, as the dispersion of specimens poorly fits the regression line. For CS, Aceratheriini and Teleoceratina form a cluster situated below the regression line, together withProtaceratherium and small Elasmotheriinae, while Rhinocerotina plot near the line. Amynodontidae, Paraceratheriidae and Rhinocerotidaeincertae sedis plot above the line, whileHyrachyus plots towards minimal CS values (Fig. 11A). Similarly, for BM, Hyrachyidae and Paraceratheriidae plot far away from the common regression slope, whereas among Rhinocerotidae, some Aceratheriini and Teleoceratina are grouped together below the line (Fig. 11B). Conversely, the regression plot for the GI-MC3 is close to those obtained on the humerus and radius, with an excellent fit to the regression line. All species are close to the common regression line, with a marked overlap between the different groups (Fig. 11C). Shape variation related to both CS, BM and GI-MC3 is highly similar and mainly concerns a mediolateral broadening towards high values, as well as a caudal development of the caudal border of the ulna (Fig. 11;Supporting Information, Fig. S3K–M). This broadening is more marked for shape variation correlated with GI-MC3.

Figure 11.

Significant PGLS regression plots for proximal partial ulna performed on shape data and log-transformed centroid size (CS) (A), log-transformed cubic root of mean body mass (BM) (B), log-transformed mean gracility index (GI-MC3) (C). Point colour coding followsFigure 1. Point size is proportional to mean log CS of each species. On the right, shapes associated with minimum and maximum fitted values (top row) and colour maps of the location and intensity of the shape deformation (bottom row). Blue: minimum value of the regression. Orange: maximum value of the regression. For each bone, the shape associated with the minimum was coloured depending on its distance to the shape associated with the maximum (blue indicates a low deformation intensity and red indicates a high deformation intensity). Orientation from left to right in each case: caudal, lateral, cranial and medial.

Open in new tabDownload slide

Evolution of CS values along the phylogeny

The evolution of CS values along the phylogeny for the distal part of the humerus, complete radius and proximal part of the ulna (these three samples being the largest) is relatively congruent between the different taxa (Fig. 12). Hyrachyidae-Hyracodontidae and giant Paraceratheriidae possess, respectively, the lowest and highest values for each bone. However, the CS of the radius shows a greater variation along the phylogeny than that of the humerus and ulna. Many taxa among Elasmotheriinae, Aceratheriini and Teleoceratina display low values relative to those observed on the humerus and ulna, these two bones displaying similar patterns of CS variation.

Figure 12.

Evolution of centroid size (CS) along the phylogeny for the studied species. A, distal partial humerus; B, radius C, proximal partial ulna. Computations were made on log-transformed CS. Values at nodes and along branches were reconstructed based on a Brownian motion model of evolution (Revell, 2012). Colour coding for taxa followsFigure 1.

Open in new tabDownload slide

Discussion

Relations between bone shape and mass, size and gracility

Our results highlight the strong relations existing between the shape variation of the forelimb bones and the changes in bone size, body mass and degree of gracility in Rhinocerotoidea, confirming the first hypothesis. However, these relations appear complex and variable depending on the bone, the anatomical area and the considered parameter, resulting in congruent and non-congruent changes along the limb.

Congruent shape variation associated with all variables

Centroid size appears almost always significantly correlated with body mass, despite missing data, heterogeneous weight estimations and marginally non-significant results for the complete ulna. This suggests that the CS of the long bones is relevant to approximate the weight of a species (Ercoli & Prevosti, 2011;Cassiniet al., 2012;Botton-Divetet al., 2017), at least on rhinocerotoids despite their diversity of body size and shape. However, beyond this general strong correlation, the variation of CS along the phylogeny for the radius differs from that observed for the humerus and ulna. Some groups may also strongly differ from the general trend shown by the whole superfamily because of specific morphological changes (i.e., Teleoceratina) (see below). Conversely, while BM correlates with GI-MC3, the latter is poorly related to CS except for the distal part of the humerus (and marginally for the complete humerus and proximal part of the ulna). This highlights that, beyond the significant correlation between bone size and body mass, these parameters do not vary conjointly with the degree of brachypody among the superfamily.

The complete humerus, distal humerus and proximal ulna share strong similarities in having their shape variation always correlated with CS, BM and GI-MC3. An increase of these variables is always associated with an increase of the bone robustness, confirming previous observations on modern (Malletet al., 2019,2020) and fossil rhinocerotoids (Prothero & Sereno, 1982;Etienneet al., 2020). Other areas mainly impacted by shape modification across the superfamily are epiphyses, which mainly extend in the mediolateral direction in heavy species. These global changes tend to indicate the existence of a common trend within the entire superfamily Rhinocerotoidea for these bones, where shape varies relatively congruously with size, mass and gracility despite the morphological diversity of these species. Convergent shape tendencies may be observed among both tetradactyl and tridactyl species (seeSupporting Information, Fig. S5, for visualization of shape variation depending on the digit number). However, as digit number is invariable within each family, we were unable to evaluate the extant of convergence.

The shape changes linked to size, mass or gracility are particularly congruent on the humerus, affecting mainly the medial side of the bone, notably the lesser tubercle tuberosity where the m. subscapularis inserts and the midshaft where the m. teres major inserts, and the m. latissimus dorsi; these muscles acting as adductors and extensors of the arm (Etienneet al., 2021). On the lateral side, most shape changes are located on the deltoid tuberosity, where the m. deltoideus inserts (Etienneet al., 2021), being more laterally developed and more distally situated on the shaft with high values of body mass, centroid size and gracility index (with a maximum for GI-MC3; see below). This distal displacement of the mm. deltoideus and of the teres major is coherent with an increase in strength of the lever arm for arm flexion and extension required to move a heavier body and limbs (Hildebrand, 1974;Polly, 2007). Such a distal displacement observed simultaneously among taxa opposed in gracility and body mass (like Teleoceratina and Paraceratheriidae) can appear paradoxical. A longer and stronger lever arm in large Paraceratheriidae is likely related to longer and heavier limbs requiring more strength to be moved. Similarly, this condition in highly brachypodial taxa may result from a difference in mass repartition: a lower centre of gravity associated with a relatively high body mass and small limbs requires powerful muscles with strong insertions to move efficiently (Hildebrand, 1974;Coughlin & Fish, 2009;Biewener & Patek, 2018). Similar observations can be done for the distal epiphysis, where most of the changes are located on the medial and lateral epicondyles and the epicondylar crest when mass, size and brachypody increase. These changes are likely associated with the development of powerful muscles for the extension movements of carpals and digits (Fisheret al., 2007;Barone, 2010a;Etienneet al., 2021) and can relate to changes in mass repartition and position of the centre of gravity as well.

Contrary to what is observed for the humerus and the proximal ulna, the shape variation of the radius and ulna (complete and without the olecranon) are only significantly correlated with body mass and gracility index (although marginally for the former parameter for the ulna). Both the radius and ulna show a reduction of the craniocaudal curvature and a straightening of the shaft with increasing body mass and brachypody. These changes are coherent with modifications observed on the humerus, highlighting the necessity to resist both higher pressure forces and stronger bending in brachypodial species (Bertram & Biewener, 1992;Milne, 2016;Hendersonet al., 2017). On the ulna, the congruent changes observed along the caudal edge of the bone towards high body mass and degree of brachypody are likely linked to a modification of the orientation of olecranon tuberosity (see below).

Non-congruent shape variation associated with variables

Beyond congruent shape variations between bones or body proportions (size, mass and gracility), some anatomical areas appear to vary more in association with one particular variable. On the humerus, this is likely the case of the bicipital groove, which is reoriented cranially and becomes more symmetrical with the apparition of an intermediate tubercle for high body mass only. This conformation is likely to play a role as a ‘passive stay-apparatus’, a feature convergently also present in modern and some fossil horses, reducing the muscular energy needed to stand for long periods (Hermanson & MacFadden, 1992). A relatively developed intermediate tubercle is observed in many groups showing high body mass (Paraceratheriidae, Aceratheriini, Rhinocerotina, Teleoceratina and, to a lesser extent, Amynodontidae), indicating the presence of a partially or fully functional passive stay-apparatus in these heavy species. Although this feature in horses is associated with a cursorial condition, equids spending long periods of time in a standing pose, its development among Rhinocerotoidea appears mainly related to their body mass.

A pronounced development of the radial tuberosity, where the m. biceps brachii inserts (Etienneet al., 2021), is observable on the radius. This development is only associated with body mass increase. This may be related to the strong flexion forces exerted by this muscle on the radius, likely related to the strength needed to move heavier limbs in large taxa (or a relatively short limb in species with a low centre of gravity). Moreover, the m. biceps brachii is also a relevant muscle involved in the passive-stay apparatus of the shoulder joint (Hermanson & MacFadden, 1992). The development of the radial tuberosity in association with body mass only is therefore coherent with changes observed on the humeral bicipital groove for the same variable.

On the ulna, the lateral border of the shaft shows a marked variation associated only with a high degree of brachypody. This area corresponds to the insertion of the lateral digit extensors (Etienneet al., 2021) and its development is coherent with that observed on the epicondylar crest of the humerus (see above). As for other extensors previously described, the marked development of these insertions along the ulna in brachypodial species may relate to the lowering of the centre of gravity and the higher power needed to move a short-limbed body efficiently.

Congruent variations between bones

Congruent shape variations are also observed between bones, which partially weakens the second hypothesis. The tricipital line running from the deltoid process to the humeral head on the humerus is particularly affected by changes in the degree of brachypody. This area corresponds to the insertion of the lateral head of the m. triceps brachii (Etienneet al., 2021). On the proximal ulna, an increase of size and mass, but above all of brachypody, involves morphological changes of the olecranon tuberosity, where the tendon of the m. triceps brachii also inserts (Etienneet al., 2021), one of the most powerful extensors of the forelimb ensuring stance of the body and opposing to gravity (Watson & Wilson, 2007;Barone, 2010b;Etienneet al., 2021). Furthermore, the development of its insertion is associated with a reorientation of the whole olecranon towards high body mass and degree of brachypody. These changes indicate a wider angle for elbow opening and a modification of the angulation of the olecranon process relatively to the shaft, known to strongly change with body mass among quadrupeds (Jenkins, 1973;Fujiwara, 2009;Fujiwara & Hutchinson, 2012;Milne, 2016;Hendersonet al., 2017).

Similarly, the distal trochlea of the humerus undergoes strong changes linked simultaneously to increases in mass, size and brachypody, becoming asymmetrical, wider and flattened, with a drastic reduction of the capitulum and a huge development of the medial lip. This conformation responds to changes observed on the radius and ulna when mass and brachypody increase. The proximal articular surfaces of the radius and ulna, forming the trochlear notch, lose their asymmetry and concavity in brachypodial taxa. Such coherent changes of the elbow region confer more degrees of freedom in the mediolateral direction, contrary to the structure encountered in light and cursorial rhinos only allowing craniocaudally constrained movements. This likely allows the elbow joint to support stronger constraints in multiple directions due to heavy weight (Polly, 2007). Such changes are coherent with similar modifications observed on the ankle joint of Rhinocerotoidea (Etienneet al., 2020), but also with observations made on modern rhinocerotoids (Malletet al., 2019), indicating a development of the medial parts of the limb bones over the lateral ones for heavier species. All these morphological modifications in the elbow region, directly linked to a higher mass in heavy taxa, may relate to a lowering of the centre of gravity of the animal in brachypodial species, involving more muscle power and longer lever arms when associated with shorter limb segments for a given mass (Hildebrand, 1974).

Differences between the stylopodium and the zeugopodium

Beyond these coherent changes located on precise anatomical areas, the patterns of shape variations appear different between the stylopodial and the zeugopodial elements. While the variations of the humeral shape follow a trend common to the whole superfamily and are simultaneously related to size, mass and gracility, those of the radius and the ulna are highly related to the degree of brachypody, coupled with an effect of the body mass being more marked on the radius than on the ulna. This relation between shape and brachypody is strikingly high for the radius. All these results likely indicate a deep functional breakdown between the stylopodium and the zeugopodium. This is coherent with an increase of the variation of limb elements along a proximodistal gradient, as hypothesized by previous authors (Hallgrímssonet al., 2002;Young & Hallgrímsson, 2005). Thanks to its obliquus orientation in the limb, the humerus ensures weight support by allowing the dissipation of stresses, while also being the support for muscles linked both to the pectoral girdle and the carpals. Therefore, it ensures the flexion and extension of the whole forelimb (Polly, 2007). Conversely, the radius and the shaft of the ulna, oriented vertically, are strongly aligned with pressure constraints due to gravity. The proximal articular surface of the radius supports the entirety of the humerus and, consequently, a significant part of the body weight—the forelimb itself supporting a larger proportion of the total weight than the hind limb (Henderson, 2006;Regnaultet al., 2013;Stilsonet al., 2016;Panagiotopoulouet al., 2019). However, our results highlight that the zeugopodial shape is more highly related to variations of brachypody than of weight, underlining the importance of the repartition of mass in the body and the position of the centre of gravity, rather than to the absolute body mass itself. The influence of the body mass value itself is more visible at lower taxonomic levels (i.e. within families or subfamilies), as has been observed among modern rhinos (Malletet al., 2019,2020).

Modularity of the elbow joint

Beyond the congruences previously described between the humerus and the ulna, the exploration of the shape of both complete and partial bones, driven at first by taphonomic constraints, led to unexpected functional observations. Whereas the complete and distal humeri show similar results, strong differences occur between the whole ulna and its proximal part (in their relations between shape, size, mass and gracility), while the shaft and the distal part seem to follow the same pattern as the radius. The proximal part of the ulna displays similar patterns of variation as the humerus (complete and distal), its shape being not only linked to gracility as in the complete ulna, but also to mass and size. This is particularly visible in Paraceratheriidae, whose complete ulna is close to the plesiomorphic condition, but whose proximal part of the ulna shows a derived morphology coherent with that of the humerus. Additional analyses on the isolated proximal part of the radius do not show this morphological shift and led to results highly similar to those obtained on the complete radius (C.M., pers. obs.). The elbow is known as a simple yet crucial hinge joint among quadrupeds, involved both in locomotion and stability of the body (Jenkins, 1973;Fujiwara, 2009;Fujiwara & Hutchinson, 2012). The humerus and ulna share complementary articular surfaces and are connected by numerous muscles (m. anconeus and flexor and extensor muscles of the carpals and digits) and a strong joint cap (Barone, 2010a;Etienneet al., 2021). Consequently, the humerus and ulna are strongly integrated among quadrupeds, i.e. they show a noticeable shape covariation (Fabreet al., 2014;Martín-Serraet al., 2015;Hanotet al., 2017;Botton-Divetet al., 2018), notably among modern rhinos (Malletet al., 2020). Our results indicate that this covariation is likely to mainly concern the distal part of the humerus and the proximal part of the ulna, leading to consider the elbow as a probable modular structure among Rhinocerotoidea, i.e. an anatomical unit covarying more in itself than with other units (Klingenberg, 2008). Beyond purely functional requirements, this potential modularity can also be related to an evolutionary covariation of size, mass and gracility among Rhinocerotoidea. Similar observations have been highlighted in small carnivorans (Fabreet al., 2014) and this assertion remains yet to be tested on modern and fossil rhinocerotoids through modularity tests (Goswami & Polly, 2010).

Bone shape and phylogenetic relationships

In addition to functional requirements, the evolutionary legacy between species has a strong but unequally distributed influence on the shape variation of the forelimb. Shape, size, mass and the degree of brachypody all carry a strong phylogenetic signal underlining that their variation is constrained by historical factors (Cubo, 2004). This influence is particularly visible on the humerus: most of the considered groups display a marked shape homogeneity despite variation in body proportions. This is not the case for the radius and the ulna, where the different groups are split depending more on their mass or degree of brachypody rather than their phylogenetic affinities. This is coherent with previous results on modern rhinos indicating that the shape of the stylopodium is more related to phylogeny than that of the zeugopodium (Malletet al., 2019,2020). This pattern seems to occur at the level of the whole superfamily, in accordance with the hypothesis of an increase of variation of the limb elements along a proximodistal gradient (Hallgrímssonet al., 2002;Young & Hallgrímsson, 2005).

However, one particular group does not seem to follow this general trend. While being closely related to stem clades like Hyracodontidae, giant Paraceratheriidae exhibit a humeral shape close to that of more derived groups like Aceratheriini. Marked shape changes relative to the shape displayed by Hyracodontidae or Hyrachyidae are observable on the humerus. Conversely, the shapes of the radius and ulna (except for the proximal part of the latter) appear to retain a plesiomorphic condition close to that of small Hyrachyidae and Hyracodontidae, these bones displaying few morphological changes except their striking relative size. These observations underline the particularity of this group among Rhinocerotoidea, whose unique body shape has puzzled biologists since their discovery (Granger & Gregory, 1936;Fortelius & Kappelman, 1993;Prothero, 2013). These considerations appear contradictory with our previous findings indicating that the radius shape is strongly related to the degree of brachypody and poorly related to phylogeny (and conversely for the humerus). It is possible that Paraceratheriidae underwent particular developmental processes constraining zeugopodium shape, while the stylopodium was subject to marked morphological changes to ensure its role in body support and propulsion, constituting a unique pattern in the superfamily. Ecological factors may also have a role in shaping the forelimb of Paraceratheriidae but this question remains to be addressed in a dedicated study.

Two other groups show marked differences with the common trend of shape variation among Rhinocerotoidea: Rhinocerotini. Species belonging to Teleoceratina likeTeleoceras show a high degree of brachypody and their forelimb bones often display an extreme shape relatively to the whole superfamily, particularly on the zeugopodium. Their extreme brachypody had sometimes been associated with a semi-aquatic ecology, although this hypothesis is now considered unlikely (MacFadden, 1998;Mead, 2000;Mihlbachler, 2003;Prothero, 2005;Clementzet al., 2008;Wang & Secord, 2020). Nevertheless, the extreme limb bone variation observed inTeleoceras,Prosantorhinus Heissig, 1973 andBrachypotherium, and their high abundance in Miocene fluviolacustrine settings with respect to coeval rhinos is compatible with that of semi-aquatic hippos (P.-O. Antoine, pers. comm.) (Harrison & Manning, 1983;Wermelinger, 1998;Antoine, 2002;Hullot & Antoine, 2020). Despite the unique limb morphology of Teleoceratina, our results highlight many shape resemblances with fully terrestrial Aceratheriini (Aphelops,Peraceras) and Rhinocerotina (Coelodonta) and do not support the hypothesis of a semi-aquatic ecology either. A morphofunctional analysis focused on this subtribe could help to understand the factors driving this particular limb construction.

Finally, Rhinocerotina (i.e., the subtribe encompassing all living rhinoceroses) display a high shape homogeneity, particularly on the humerus and the ulna, despite a broad range of body mass and body proportions. The range of shape variation within this subtribe appears thus highly constrained by the evolutionary history. The ecological preferences encountered in Rhinocerotina do not seem to strongly impact the shape variation (Guérin, 1980;Cerdeño, 1998). However, this relative homogeneity in the whole superfamily likely encompasses different trends of shape variation between genera that remain to be explored in detail.

CONCLUSION

The relations between shape variation of the forelimb bones, body proportions and phylogeny among Rhinocerotoidea vary, but general trends are clearly observed despite this complexity. A common trend in the superfamily is the increase of bone robustness towards a higher body mass and higher degree of brachypody. The reinforcement of the insertions for the extensor muscles enables the animals to counteract the gravitational constraints when body mass increases. However, strong differences in shape variation exist between the stylopodium and the zeugopodium. The shape of the humerus is modified following size, mass and brachypody in a similar way in the whole superfamily, while being also strongly constrained by the evolutionary history. Conversely, the shape of the zeugopodium appears mostly associated with the degree of brachypody, namely the distribution of mass in the body (centre of gravity), rather than with the absolute mass itself. Surprisingly, the shape variation of bones in the elbow caudal region shows striking similarities, suggesting a likely modular organization of the humerus and ulna. Beyond these general trends, groups like Paraceratheriidae, Teleoceratina and Rhinocerotina display divergent patterns that remain to be fully understood. Consequently, this exploration of the forelimb shape among Rhinocerotoidea encourages the application of the same morphofunctional approach on the hind limb to highlight how shape patterns converge or diverge between limbs under a similar weight constraint.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article at the publisher’s web-site.

Data S1. Designation and location of the anatomical landmarks placed on each bone.

Figure S1. Summary of the anatomical areas of the rhino long bone. Bones figured here belong toCeratotherium simum. A, humerus. Abbreviations: B.g.: bicipital groove; C.: capitulum; D.t.: deltoid tuberosity; E.c.: epicondylar crest; G.t.: greater tubercle; G.t.c.: greater tubercle convexity; h.: Head; I.t.: intermediate tubercle; L.e.: lateral epicondyle; L.l.b.: lateral lip border; L.t.: lesser tubercle; L.t.c.: lesser tubercle convexity; M.e.: medial epicondyle;M.i.i.: m. infraspinatus insertion; M.l.b.: medial lip border;M.t.m.t.: m. teres major tuberosity; N.: neck; O.f.: olecranon fossa; T.: trochlea; T.g.: trochlear groove. B, radius. Abbreviations: A.s.s.: articular surface for the scaphoid; A.s.sl.: articular surface for the semilunar; C.p.: coronoid process; D.a.s.u.: distal articular surface for the ulna; I.c.: interosseous crest; I.s.: interosseous space; L.g.c.: lateral glenoid cavity; L.i.r.: lateral insertion relief; L.s.a.s.: lateral synovial articular surface; M.g.c.: medial glenoid cavity; M.s.a.s.: medial synovial articular surface; P.a.s.u.: proximal articular surface for the ulna; P.p.: palmar process; R.s.p.: radial styloid process; R.t.: radial tuberosity. C, ulna. Abbreviations: A.p.: anconeal process; A.s.h.: articular surface for the humerus; A.s.p.: articular surface for the pisiform; A.s.sl.: articular surface for the semilunar; A.s.t.: articular surface for the triquetrum; D.a.s.r.: distal articular surface for the radius; I.c.: interosseous crest; I.s.: interosseous space; M.t.o.: medial tuberosity of the olecranon; O.t.: olecranon tuberosity; P.b.: palmar border; U.s.p.: ulnar styloid process.

Figure S2. Shape deformations associated with the first two axes of the PCA for each bone. Blue: minimal values. Orange: maximal values. Orientation from left to right: caudal, lateral, cranial, medial, proximal and distal views. (A) complete humerus; (B) distal partial humerus; (C) radius; (D) complete ulna; (E) ulna without olecranon tuberosity; (F) proximal partial ulna.

Figure S3. Shape deformations associated with minimum and maximum values of the centroid size (CS), body mass (BM) and gracility index (GI-MC3) for significant regressions with shape. Blue: minimal values. Orange: maximal values. Orientation from left to right: caudal, lateral, cranial, medial, proximal and distal views. (A, B, C) complete humerus; (D, E, F) distal partial humerus; (G, H) radius; (I) complete ulna; (J) ulna without olecranon tuberosity; (K, L, M) proximal partial ulna.

Figure S4. Significant PGLS regression plots of ulna without olecranon tuberosity performed on shape data and log-transformed cubic root of mean body mass (BM). Point colour coding followsFigure 1. Point size is proportional to mean log CS of each species. On the right, shapes associated with minimum and maximum fitted values (top row) and colour maps of the location and intensity of the shape deformation (bottom row). Blue: minimal values. Orange: maximal values. For each bone, the shape associated with the minimum was coloured depending on its distance to the shape associated with the maximum (blue indicates a low deformation intensity and red indicates a high deformation intensity). Orientation from left to right in each case: caudal, lateral, cranial and medial views.

Figure S5. Number of digits for each species plotted on the Neighbour Joining trees computed on all PC scores as inFigure 3. Colour coding of species names followsFigure 1 and abbreviations followTable 1. Colour coding for number of digits as indicated on the bottom of the figure. Point size is proportional to the mean log centroid size of each species. (A) complete humerus; (B) distal partial humerus; (C) radius; (D) complete ulna; (E) ulna without olecranon tuberosity; (F) proximal partial ulna.

Table S1. Complete list of all the studied specimens.

Table S2. Complete list of data compiled from the literature for gracility index on third metacarpal and mean body mass.

Table S3. Summary of the differences in P and R² values between the PGLS computed under a Brownian Motion (BM) model (geomorph) and a Ridge Regression (RR) model (RRphylo). Only variables with significant results are presented here.

ACKNOWLEDGEMENTS

The authors warmly thank all curators of the visited institutions for granting us access to the studied specimens, particularly E. Hoeger, S. Ketelsen, R. O’Leary and J. Meng (American Museum of Natural History, New York, USA), J.-M. Pouillon and C. Bouix (Association Rhinopolis, Gannat, France), G. Rößner (Bayerische Staatssammlung für Paläontologie und Geologie, Munich, Germany), D. Berthet (Centre de Conservation et d’Étude des Collections, Musée des Confluences, Lyon, France), E. Robert (Collections de Géologie de Lyon, Université Lyon 1 Claude Bernard, Lyon, France), Y. Laurent (Muséum d’Histoire Naturelle de Toulouse, Toulouse, France), J. Lesur, A. Verguin (Muséum National d’Histoire Naturelle, Paris, France), R. Portela-Miguez, P. Brewer and R. Pappa (Natural History Museum, London, UK), L. Costeur and F. Dammeyer (Naturhistorisches Museum Basel, Basel, Switzerland), A. Folie, C. Cousin, O. Pauwels and S. Bruaux (Royal Belgian Institute of Natural Sciences, Brussels, Belgium), E. Gilissen (Royal Museum for Central Africa, Tervuren, Belgium) and D. Brinkman (Yale Peabody Museum, New Haven, CT, USA). We also thank W. Liu from the Institute of Vertebrate Paleontology and Paleoanthropology (Beijing, China), for providing a 3D model ofJuxia, M.C. Reyes from the National Museum of the Philippines (Manila, Philippines) and T. Ingicco from the MNHN (Paris, France) for providing the 3D models ofN. philippinensis and J. Hutchinson from the Royal Veterinary College (London, UK) for providing us with CT-scan data from the University of California Museum of Paleontology (Berkeley, USA). We would like to warmly thank P.-O. Antoine and an anonymous reviewer for their constructive comments that allowed us to greatly improve the quality of the final manuscript. We are grateful to S. Castiglione and P. Raia (University of Naples Federico II, Naples, Italy) for their help in using the RRphylo package. Many thanks to K. Gaignebet and C. Bouquet for their help in reconstructing many 3D models. C.M. acknowledges C. Étienne, R. Lefebvre and R. Pintore (MNHN, Paris, France) for constructive discussions and advice on R programming, data analyses and interpretation. The authors declare that there is no conflict of interest.

FUNDING

This work was funded by the European Research Council and is part of the GRAVIBONE project (ERC-2016-STG-715300).

AUTHOR CONTRIBUTIONS

C.M. designed the study with significant inputs from A.H., R.C. and G.B. C.M. acquired the data with inputs from A.H. C.M. performed the analyses with the help of R.C and G.B. and all authors interpreted the results. C.M. drafted the manuscript. All authors reviewed and contributed to the final version of the manuscript, read it and approved it.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. Most of the 3D models will be or have been deposited on the 3D online repository MorphoSource (https://www.morphosource.org/projects/000366286?locale=en).

References