Interspecific variation in the limb long bones among modern rhinoceroses—extent and drivers

Landmark digitization

Following the procedure described by Gunz, Mitteroecker & Bookstein (2005), Gunz & Mitteroecker (2013) and Botton-Divet et al. (2016), we defined the bones’ shape using anatomical landmarks and curve and surface sliding semi-landmarks. Each curve is bordered by anatomical landmarks as recommended by Gunz & Mitteroecker (2013). We placed all landmarks and curves using the IDAV Landmark software (v3.0—Wiley et al., 2005). We used 35 anatomical landmarks on the humerus, 23 on the radius, 21 on the ulna, 27 on the femur, 24 on the tibia and 12 on the fibula. Details of landmark numbers and locations used for each bone are given in Table S1 and Fig. S1.

Following the procedure detailed by Botton-Divet et al. (2016), we created a template to place surface semi-landmarks for each bone: a specimen (C. simum RMCA 1985.32-M-0001) was randomly chosen on which all anatomical landmarks, curve and surface sliding semi-landmarks were placed. We then used this template for the projection of surface sliding semi-landmarks on the surface of the 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 two times between the result of the preceding step and the Procrustes consensus of the complete dataset. Therefore, all landmarks can be treated at the end as geometrically homologous (Gunz, Mitteroecker & Bookstein, 2005) and analyzed with classic procedure such as generalized Procrustes analysis (GPA; see below). Projection, relaxation and sliding processes were conducted using the Morpho package in the R environment (R Development Core Team, 2014). Details of the process are provided in the documentation of the package (Schlager, 2018).

Repeatability tests

For each bone, we tested the repeatability of the anatomical landmark digitization taking measurements ten times on three specimens of the same species, C. simum, chosen to display the closest morphology and size. We superimposed these measurements using a GPA and visualized the results using a principal component analysis (PCA). Results showed a variation within specimens clearly smaller than the variation between specimens (see Fig. S2) and allowed us to consider our anatomical landmarks as precise enough to describe shape variation.

Generalized Procrustes Analyses

After the sliding of all semi-landmarks, we performed GPA (Gower, 1975; Rohlf & Slice, 1990) to remove the effects of size and of the relative position of the points and to isolate only the shape information. As our dataset contained more variables than observations, we used a PCA to reduce dimensionality as recommended by Gunz & Mitteroecker (2013) and visualize the specimen distribution in the morphospace. We computed theoretical consensus shape of our sample and used it to calculate a TPS deformation of the template mesh. We then used this newly created consensus mesh to compute theoretical shapes associated with the maximum and minimum of both sides of each PCA, as well as mean shapes of each bone for each species. GPA, PCA and shape computations were done using the “Morpho” and “geomorph” packages (Adams & Otárola-Castillo, 2013; Adams, Collyer & Kaliontzopoulou, 2018; Schlager, 2018) in the R environment (R Development Core Team, 2014). Neighbor Joining method was used to construct trees based on relative Euclidian distances between individuals based on all principal component scores obtained with the PCA, allowing a global visualization of the relationships between all the specimens. Trees were computed with the “ape” package (Paradis, Claude & Strimmer, 2004).

Allometry effect

We tested the effect of allometry, defined as “the size-related changes of morphological traits” (Klingenberg, 2016). Pearson’s correlation tests were performed to look for correlation between the principal components and the centroid size (log₁₀) for each bone. We also used the function procD.allometry of the “geomorph” package to perform a Procrustes ANOVA (a linear regression model using Procrustes distances between species instead of covariance matrices—see Goodall, 1991) to quantify the shape variation related to the centroid size, and to visualize theoretical shapes associated with minimal and maximal sizes of our sample (Adams & Otárola-Castillo, 2013; Adams, Collyer & Kaliontzopoulou, 2018). This test was performed taking into account group affiliation (e.g., species) to highlight respective roles of centroid size and species determination on the shape variation. In the absence of individual body mass for the majority of our sample, we also performed a Procrustes ANOVA with the cube root of the mean mass attributed to each species (Table 1), each specimen being associated with the mean mass of its species. As for the centroid size, theoretical shapes associated with minimal and maximal mean mass were computed using the predicted Procrustes residuals (details on the procedure are given in the “geomorph” documentation). Plots of the multivariate regressions of shape scores (i.e., regression of shape on size; see Drake & Klingenberg, 2008) against log-transformed centroid size were also computed.

Humerus

The first two axes of the PCA computed on the humerus represent 60.6% of the total variance (Fig. 1A). The first axis represents more than half of the global variance (53%) and the five species appear clearly sorted along it, opposing Dicerorhinus sumatrensis on the positive side to C. simum on the negative one, i.e., the lightest and heaviest species, respectively. Diceros bicornis is grouped with C. simum on the negative part of the axis, whereas R. sondaicus is on the positive part. R. unicornis occupies the center of the axis, between Diceros bicornis and R. sondaicus. Points distribution in the morphospace and Neighbor Joining trees indicate a clear separation between African and Asiatic rhinos (Fig. S4A). The theoretical shape at the PC1 minimum (Figs. 1B, 1D, 1F and 1H) shows a massive morphology, with mediolaterally and craniocaudally broad epiphyses and shaft; a wide humeral head, with very little overhanging of the diaphysis in the caudal direction; a lesser tubercle more strongly developed than the greater tubercle, with an intermediate tubercle separating a widely open bicipital groove into unequal parts, the lateral one being the largest; a lesser tubercle convexity medially extended whereas the greater tubercle one is quite reduced in this direction; a broad and diamond-shaped m. infraspinatus imprint on the lateral side; a broad deltoid tuberosity not extending beyond the lateral border of the bone; a shaft with its maximal width situated between the humeral head and the deltoid tuberosity; a distinct but very smooth and flat m. teres major tuberosity; a very large distal epiphysis because of the development of the lateral epicondyle; a smooth epicondylar crest; a mediolaterally wide and craniocaudally compressed medial epicondyle; shallow and proximodistally compressed olecranon fossa and trochlea; a wide trochlea displaying a main axis tilted in the dorsoventral direction; and a capitulum with a relatively small surface area. At the opposite, the theoretical shape at the PC1 maximum (Figs. 1C, 1E, 1G and 1I) shows a slender and thin aspect; a more rounded humeral head overhanging the diaphysis caudally; a greater tubercle more strongly developed than the lesser one and extending medially, conferring a more closed aspect to the bicipital groove, where the intermediate tubercle is almost absent; a poorly developed lesser tubercle convexity whereas the greater tubercle one is massive; a rounded and reduced m. infraspinatus insertion; a deltoid tuberosity strongly protruding laterally; a straight and thin shaft; no visible m. teres major tuberosity; a narrow distal epiphysis, with a small development of the lateral epicondyle; a sharp epicondylar crest; a craniocaudally developed medial epicondyle overhanging the olecranon fossa; a deep and wide olecranon fossa; a far less compressed trochlea, with a less dorsoventrally tilted axis; and a very reduced capitulum.

Along the second axis (7.6%), we observe this time that C. simum and Dicerorhinus sumatrensis are grouped together on the negative part of the axis, with the three other species on the positive part, whereas they are opposed along the first axis. This second axis expresses the separation between the lightest and the heaviest rhino species on the one hand and the three other species on the other hand. The theoretical shape at the PC2 minimum displays a humeral head stretched in the caudal direction; a lesser tubercle more developed than the greater one, delimiting an open bicipital groove; a proximodistally extended distal epiphysis, with an epicondylar crest starting almost on the middle of the shaft; a rounded and wide olecranon fossa. At the opposite, the theoretical shape at the PC2 maximum shows a rounded humeral head; a strong development of both tubercles and a more closed bicipital groove; a mediolaterally stretched distal epiphysis, with the epicondylar crest starting at the distal third of the shaft; an olecranon fossa proximodistally compressed and more rectangular; and a well-developed lateral epicondyle.

Radius

The first two axes of the PCA performed on the radius express 52.3% of the total variance (Fig. 2A). The first axis (36.4%) opposes Dicerorhinus sumatrensis and Diceros bicornis to R. unicornis and C. simum. R. sondaicus overlaps both R. unicornis and Diceros bicornis clusters. The specimens of Dicerorhinus sumatrensis are split in two discrete clusters along the first axis, but no clear explanation linked to age, sex or geographic origin was associated with this distribution. Point dispersion along this axis indicates an important intraspecific variation for Dicerorhinus sumatrensis, and to a lesser extent for Diceros bicornis and R. sondaicus. Unlike for the humerus, phylogenetically related species are not grouped together on PCA and Neighbor Joining trees (Fig. S4B). The theoretical shape at the PC1 minimum (Figs. 2B, 2D, 2F and 2H) shows a massive morphology with large shaft and epiphyses; an asymmetrical proximal articular surface (constituting the ulnar notch), with a medial portion appearing nearly twice as large as the lateral one; a protruding lateral insertion relief (distally to the lateral coronoid process sensu Budras, Sack & Röck, 2009) (i.e., insertion area of the m. extensor digitorum communis) whereas the radial tuberosity is little prominent; a mediolaterally reduced lateral synovial articulation surface for the ulna; a rectangular and thin medial synovial articulation surface for the ulna; a triangular proximal articular surface for the ulna as wide mediolaterally as proximodistally; a thick shaft with an interosseous space opening close to the proximal epiphysis: consequently, the interosseous crest runs along the diaphysis to the distal articular surface for the ulna; a broad distal epiphysis in the mediolateral direction, with a strong medial tubercle developed on the dorsal face; a distal articular surface compressed in the dorsoventral direction; an articular surface for the scaphoid little extended proximally; a trapezoidal and wide articular surface for the semilunar (i.e., lunate bone or lunatum); a well-developed radial styloid process. The theoretical shape at the PC1 maximum (Figs. 2C, 2E, 2G and 2I) displays a more slender morphology; a less asymmetrical proximal articular surface, despite the development of the medial part; an almost absent lateral insertion relief; a completely flat radial tuberosity; a mediolaterally reduced lateral synovial articulation surface for the ulna; a rectangular and thin medial synovial articulation for the ulna; a triangular proximal articular surface for the ulna, mediolaterally short and proximodistally stretched; a thin and slender shaft, with an interosseous space opening at the proximal third of the total length; a poorly visible interosseous crest; a less dorsoventrally compressed distal epiphysis and a poorly developed lateral tubercle (on the dorsal side); a dorsoventrally wide distal articular surface, with a proximally extended articular surface corresponding to the scaphoid; a trapezoidal and reduced articular surface for the semilunar; a less developed radial styloid process with a rounded border.

The second axis (15.9%) discriminates mainly R. sondaicus from the four other species. R. unicornis displays little extension along this axis; neither does Diceros bicornis, only driven on the negative side by a single individual. Extension of R. unicornis morphospace occupation along the second axis is very limited, contrary to that of C. simum and Dicerorhinus sumatrensis clusters. As on the first axis, Dicerorhinus sumatrensis is split in two clusters, one in the negative part and the other around null values. The theoretical shape at the PC2 minimum displays a slender morphology, with a strongly asymmetrical proximal articular surface; a proximally reduced palmar process, opposed to the coronoid process; a distal epiphysis dorsoventrally broad, with a developed lateral prominence; a little developed radial styloid process; a slight proximal extension of the articular surface for the scaphoid. The theoretical shape at PC2 maximum displays a more massive shape; a deeper and more symmetrical proximal articular surface with a well-developed palmar process; a dorsoventrally compressed distal epiphysis with a more developed styloid process.

Ulna

The first two axes of the PCA performed on the ulna express 41.5% of the total variance (Fig. 3A). The first axis (22.1%) separates Dicerorhinus sumatrensis and Diceros bicornis on the positive part and R. sondaicus, R. unicornis and C. simum on the negative part. However, the clusters of C. simum and R. unicornis overlap along this axis. The general pattern on both PCA and Neighbor Joining trees is close to the one observed for the radius (Fig. S4C). The theoretical shape at the PC1 minimum (Figs. 3B, 3D, 3F and 3H) displays a thick morphology with large epiphyses; a massive olecranon tuberosity with a medial tubercle—where inserts the medial head of the m. triceps brachii, as well as the mm. flexor carpi ulnaris and flexor digitorum superficialis—oriented dorsally; an anconeal process poorly developed dorsally and mediolaterally wide, as is the articular surface constituting the trochlear notch (receiving the humeral trochlea); a medially stretched medial part of the articular surface for the humerus; a short interosseous crest ending at the shaft half, with the interosseous space; a broad shaft with a triangular cross-section; a straight palmar border whereas the shaft is medially curved; a massive distal epiphysis with a wide insertion surface for the radius; a mediolaterally wide and little concave articular surface for the triquetrum (i.e., triquetral, pyramidal or cuneiform bone), while the one responding to the pisiform is crescent-shaped and little extended proximally. The theoretical shape for the PC1 maximum (Figs. 3C, 3E, 3G and 3I) displays a more gracile morphology; a slender olecranon tuberosity with a medial tubercle where inserts the medial head of the m. triceps brachii oriented in the palmar direction; a dorsally developed and mediolaterally narrow anconeal process, as is the articular surface of the trochlear notch; a slightly medially stretched medial part of the articular surface; a sharp interosseous crest; a thin and straight shaft; a mediolaterally compressed and little concave distal epiphysis; a mediolaterally narrow articular surface for the triquetrum; a triangular and proximally well-developed articular surface for the pisiform.

The second axis (19.4%) separates quite clearly the three Asian species from the African ones. The theoretical shape at the PC2 minimum displays a slender and straight morphology with a high square-shaped olecranon process, mediolaterally flattened, more stretched in the palmar direction; a wide and squared anconeal process; a straight and regular shaft; a mediolaterally compressed distal epiphysis with a concave articular surface for the triquetrum and a distally developed styloid process; a proximally extended articular facet for the pisiform. The theoretical shape at the PC2 maximum displays a more massive and medially concave shape with an olecranon process mediolaterally inflated and rounded in the palmar direction; an anconeal process little developed dorsally and laterally tilted; a proximodistally compressed and extending medially articular surface constituting the trochlear notch; a mediolaterally wide articular surface for the triquetrum; a little developed styloid process; a poorly extended proximally and square-shaped articular surface for the pisiform.

Femur

The first two axes of the PCA performed on the femur express 45.0% of the global variance (Fig. 4A). The first principal component (36.1%) clearly isolates Dicerorhinus sumatrensis on the positive part from the other species. The clusters of Diceros bicornis, R. sondaicus and R. unicornis overlap on the negative part of the axis. Diceros bicornis and R. unicornis specimens overlap a substantial part of the cluster of C. simum too. The general pattern observed on the Neighbor Joining tree is closer to the humerus one, with African and Asiatic species grouped together, respectively (Fig. S4D). The theoretical shape at the PC1 minimum (Figs. 4B, 4D, 4F and 4H) shows a massive morphology with large epiphyses and a curved medial border, conferring a concave aspect to the diaphysis; a large femoral head, off-centered relatively to the shaft main axis, supported by a very large neck; a small and shallow fovea capitis mediocaudally oriented; a greater trochanter convexity expanding strongly laterodistally; the absence of trochanteric notch between the convexity and the top of the trochanter (Fig. 4F); a proximodistally reduced trochanteric fossa; a sharp lesser trochanter running along the medial edge, which is craniocaudally flattened below the humeral head; a third trochanter extending strongly laterally, cranially and proximally toward the greater trochanter convexity, and much curved toward the medial direction; a quite irregular shaft section along the bone—flattened below the proximal epiphysis and more trapezoidal toward the distal epiphysis; a broad distal epiphysis with developed medial and lateral epicondyles; a shallow supracondylar fossa; a wide trochlea, with a main rotation axis aligned with the shaft axis; a large and cranially expanded medial ridge of the trochlea separated from the lateral one by a deep trochlear groove; a medial condyle surface area larger than the lateral condyle one, both being separated by a narrow intercondylar space. At the opposite, the theoretical shape at the PC1 maximum (Figs. 4C, 4E, 4G and 4I) is more slender with a straight and regular shaft; a rounded femoral head aligned with the shaft main axis and supported by a thinner neck; a more pronounced and rounded fovea capitis oriented almost completely caudally; a greater trochanter convexity little developed laterodistally; a more pronounced trochanter top despite the absence of trochanteric notch; a thin lesser trochanter situated on the caudal border of the medial side; a rounded third trochanter more developed laterally than cranially; a quite regular and trapezoidal shaft section; a mediolaterally broader and medially oriented distal epiphysis; an almost absent supracondylar fossa; a less developed medial trochlear ridge separated from the lateral one by a shallow trochlear groove; a lateral condyle more oblique and divergent relatively to the medial one, increasing the intercondylar space; symmetrical medial and lateral condylar surfaces.

The second axis (8.9%) clearly opposes Dicerorhinus sumatrensis, C. simum and Diceros bicornis on the positive part to the two Rhinoceros species on the negative part, the cluster of Dicerorhinus sumatrensis being driven toward negative values by a single individual. The theoretical shape at the PC2 minimum is mainly characterized by a flattened femoral head with a strong neck; a rounded and large fovea capitis mediocaudally oriented; a laterodistally expanded greater trochanter convexity; a long and thin lesser trochanter; an extremely developed third trochanter in lateral, cranial and proximal directions; a straight and regular shaft; a broad distal epiphysis with important development of both epicondyles; a trochlea rotation axis aligned with the main axis of the shaft. The theoretical shape at the PC2 maximum displays a more rounded head, with a more stretched neck; no fovea capitis at all but a little groove on the head border; a greater trochanter convexity little expanded laterodistally; a short and more medially developed lesser trochanter; a rounded third trochanter little developed in cranial and proximal directions; a straight shaft; a distal epiphysis less mediolaterally broad; a narrower intercondylar space; a more inflated medial condyle.

Tibia

The first two axes of the PCA performed on the tibia express 50.0% of the global variance (Fig. 5A). The first axis (29.1%) separates roughly Diceros bicornis and Dicerorhinus sumatrensis on the positive part and C. simum, R. sondaicus and R. unicornis on the negative part. Diceros bicornis shows an important intraspecific variation along both axes. Neighbor Joining tree structure is less clear than for previous bones: both Rhinoceros species isolate from most of the other specimens, C. simum appears also separated from Diceros bicornis and Dicerorhinus sumatrensis. However, one C. simum and three Dicerorhinus sumatrensis specimens are closer from the Rhinoceros group than from their own respective species (Fig. S4E). The theoretical shape at the PC1 minimum (Figs. 5B, 5D, 5F and 5H) shows a massive morphology with broad shaft and epiphyses, both in craniocaudal and mediolateral directions; medial and lateral intercondylar tubercles having the same height and a reduced central intercondylar area; a broad cranial intercondylar area; a medial articular surface larger than the lateral one, with a caudally extended sliding surface for the m. popliteus tendon; a U-shaped popliteal notch; a rounded tibial tuberosity, laterally deflected and medially bordered by a shallow groove; a shallow extensor groove; a regularly triangular and distally extended proximal articular surface for the fibula; a thick tibial crest disappearing at the middle of the shaft, where the bone section is the smallest; a mediolaterally broad and rectangular in section distal epiphysis; a triangular-shaped distal articular surface for the fibula reduced in height, surmounted by a smooth interosseous crest running toward the middle of the shaft; a roughly rectangular distal articular surface for the talus, with a lateral groove larger and shallower than the medial one, separated by a prominent intermediate process without synovial fossa; an articular surface with a rotation axis aligned with the bone main axis; a prominent medial malleolus. The theoretical shape at the PC1 maximum (Figs. 5C, 5E, 5G and 5I) displays a relatively gracile morphology with a thin shaft; a lateral intercondylar tubercle more proximally extended than the medial one and a relatively large central intercondylar area; a cranially extended lateral condylar surface, reducing the cranial intercondylar area; roughly equal medial and lateral articular surface areas; a V-shaped popliteal notch; a tibial tuberosity slightly more laterally deflected; a deeper tuberosity groove; a nail-shaped proximal articular surface for the fibula; a sharper tibial crest disappearing just before the first half of the shaft; a distal epiphysis more compressed craniocaudally; a distal articular surface for the fibula displaying a large triangle synostosis area occupying a third of the shaft and prolonged by a sharp interosseous crest. There is no major difference in the distal articular shape between PC1 maximum and minimum, except that the caudal apophysis is less prominent in the distal direction.

The second axis (20.9%) clearly separates the two African species (C. simum and Diceros bicornis) on the positive part from the three Asian species (Dicerorhinus sumatrensis, R. sondaicus and R. unicornis) on the negative part. The theoretical shape at the PC2 minimum displays a slightly more slender morphology; a proximal plateau higher cranially than caudally and forming a closer angle with the diaphysis axis; a high intercondylar eminence; a lateral articular surface more caudally extended than the medial one; a tibial tuberosity well separated from the condyles by deep tuberosity and extensor grooves; a straight shaft ending with divergent borders forming a large and rectangular distal epiphysis; a distal articular surface for the fibula forming a regular triangle surmounted by a sharp interosseous crest; a medially extended medial malleolus, resulting in a rectangular articular surface with the talus, where the medial groove is narrow and deep, occupying a third of the area, whereas the lateral groove is shallow and broad. The theoretical shape at the PC2 maximum displays a more massive morphology, with a craniocaudal inflation of the epiphyses; a proximal plateau almost perpendicular to the diaphysis axis; a lower intercondylar eminence; a lateral condyle surface almost twice less large than the medial one, which is more developed caudally; a massive tibial tuberosity strongly deviated laterally, delimited by very shallow tuberosity and extensor grooves and resulting in a very large cranial intercondylar area; a straight shaft ending with almost parallel medial and lateral borders and a square-shaped distal epiphysis; a less medially deflated medial malleolus; a squared distal articular surface for the talus with medial and lateral grooves showing similar surface area and depth.

Fibula

The first two axes of the PCA performed on the fibula express 55.9% of the global variance (Fig. 6). Contrary to the five previous analyses, the first axis (40.7%) here seems particularly driven by a strong intraspecific variation. The clusters of C. simum and Dicerorhinus sumatrensis are stretched along the PC1 and overlap with almost every other specimens. The cluster of Diceros bicornis is quite stretched along the axis too and only the two Rhinoceros species display less intraspecific variation. This pattern does not seem linked to sex, age class or condition (wild or captive): despite the presence of slightly more females and subadults on the negative part of the component, we did not consider this observation as robust enough to state on this question. This cluster distribution along the PC1 seems linked to the presence of irregular crests along the shaft, associated with an important variation of the outline of the crests running along the diaphysis, and a slight rotation of the fibular head (see Fig. S3). Consequently, we chose to display and analyze the specimen distribution along the second and third components instead. Theoretical shapes associated with the PC1 are available in Fig. S3.

PC2 and PC3 express 22.9% of the global variance (Fig. 7A). The second component (15.2%) opposes C. simum on the negative side to Dicerorhinus sumatrensis on the positive side, whereas Diceros bicornis, R. sondaicus and R. unicornis have a more central disposition. As for the tibia, the Neighbor Joining tree structure appears less clearly sorted by species than for other bones. If Rhinoceros species group together and African ones as well, Dicerorhinus sumatrensis sample is split in two subgroups mixed with R. unicornis and African rhinos, respectively (Fig. S4F). The theoretical shape at the PC2 minimum (Figs. 7B, 7D, 7F, 7H and 7J) displays a broad morphology with large epiphyses and a straight shaft; a rounded head with a craniomedially oriented proximal articular surface for the tibia; a head width similar to the shaft one; a robust shaft with two strong craniolateral and caudolateral lines running down the distal epiphysis and enlarging craniocaudally toward the distal epiphysis; a sharp and irregular interosseous crest; a mediolaterally compressed distal epiphysis with little development of the two distal tubercles at the end of the lateral crests; a shallow lateral groove; a triangular distal articular surface for the tibia, occupying only the last distal quarter of the bone length; a short and ovoid articular surface for the talus with a sharp distal ridge. The theoretical shape at the PC2 maximum (Figs. 7C, 7E, 7G, 7I and 7K) displays a slender morphology with a strongly curved shaft; a mediolaterally flat head extending craniocaudally and overhanging strongly the diaphysis; a thin shaft with two sharp lateral crests running along it: these crests end with two developed tubercles surrounding a deep lateral groove; a distal articular surface for the tibia extending from the distal third of the shape and forming a stretched triangle; a wider and kidney-shaped articular surface for the talus, forming two distal tips corresponding to the two lateral tubercles: between them on the distal face, a large groove is visible, ending at the center of the face.

The third component (7.7%) mainly opposes Diceros bicornis on the positive part to R. sondaicus on the negative part. However, this opposition is mainly driven by a small number of individuals (two for Diceros bicornis and four for R. sondaicus). The specimens of R. sondaicus are divided into two clusters, with three individuals overlapping notably with Dicerorhinus sumatrensis. The theoretical shape at the PC3 minimum shows a massive morphology, with broad shaft and epiphyses; a craniocaudally broad head, overhanging the shaft laterally; a proximal articular surface for the tibia almost completely medially oriented; a straight shaft displaying a constant width along the bone; craniolateral and caudolateral crests running almost parallel toward the distal end of the bone, forming two developed tubercles surrounding a deep groove; an interosseous space covered by irregular reliefs and bordered by a sharp interosseous crest; a distal articular surface for the tibia forming a cranially deviated triangle; a kidney-shaped distal articular surface for the talus, with a distal border separated from the lateral tubercles by a groove stopping at the middle of the distal face. The theoretical shape at the PC3 maximum shows an extremely thin morphology with a flattened and poorly developed head; a proximal articular surface oriented almost completely in the cranial direction; a torsion of almost 90° between the orientation of the proximal and distal articular surfaces for the tibia; a very thin and flat shaft; craniolateral and caudolateral crests running along the diaphysis ending on the distal epiphysis with few developed tubercles; a distal articular surface for the tibia forming a slender triangle; a relatively small distal articular surface for the talus, with a less pronounced kidney-shape; a groove on the distal face mediolaterally compressed.

Forelimb bones

The difference between high and low size among extant rhinos is expressed on the humerus by a general enlargement in both craniocaudal and mediolateral directions, particularly for the proximal first half. This may be related to the constraints exerted both by weight bearing and braking role of the forelimb during locomotion (Dutto et al., 2006). The important development of the lesser tubercle at the expense of the greater tubercle in non-Dicerorhinus species allows both a greater stability of the shoulder articulation, preventing hyperextension and a larger insertion area for the medial end of the m. supraspinatus, also considered as a shoulder stabilizer (Fisher, Scott & Naples, 2007; Watson & Wilson, 2007). This muscle being one of the main extensors of the forelimb (Barone, 2010b), the developed lesser tubercle acts as a strong medial lever arm for extension movements. This configuration has been previously interpreted as mechanically advantageous for the muscles inserting on the shoulder joint, while the lateral reinforcement of the greater tubercle was supposed to help resisting the adduction of the arm (Hermanson & MacFadden, 1992). The development of the lesser tubercle may also help supporting the scapula (more elongated among African rhinos, J. MacLaren, 2019, personal communication) and be associated with a lengthening of the m. subscapularis tendons. In addition, the lesser tubercle also displays an important development in Diceros bicornis, more pronounced than in R. unicornis and R. sondaicus, though these species are heavier and taller. This indicates a possible effect of phylogenetic proximity or similar habitats between the African species (see above). The development of the intermediate tubercle for some rhinos may be related to the presence of a forelimb passive stay apparatus, as demonstrated in horses (Hermanson & MacFadden, 1992; Mihlbachler et al., 2014). Although less developed than in equids, the intermediate tubercle is present in all rhinos at different degrees (well visible in African taxa, less developed in Rhinoceros species and poorly developed in Dicerorhinus). This may indicate different degrees of development of passive stay mechanism possibly linked to phylogeny and ecology (Shockey, 2001). On the distal epiphysis of the humerus, the mediolateral enlargement observed toward high body mass ensures both a greater stability of the elbow articulation and larger insertion areas for the different flexor and extensor muscles for the digits (Barone, 2010a). The distal trochlea of the humerus is also subjected to a proximodistal compression and a mediolateral extension, increasing the articular surface area to dissipate compressive forces, important for maintaining posture at high body masses (Jenkins, 1973).

Forelimb paired zeugopodial bones seem to express complementary shape variations linked to body mass. Whereas the radius broadens mainly mediolaterally with increasing body mass, the ulna expands in the craniocaudal direction; they respond conjointly to the increase in body mass and bone size to form a structure reinforced in all directions, as it has been observed on the humerus. All rhinos have an ulnar proximal epiphysis situated caudally to the radius, while its shaft expands laterally, possibly allowing a mediolateral weight distribution. Moreover, almost all the weight is borne by the proximal articular surface of the radius (Bertram & Biewener, 1992), which expands medially and becomes asymmetrical for heavier rhinos. The concave radial tuberosity shows a deep m. biceps brachii insertion delivering a strong forearm flexion (Antoine, 2002) and the developed insertion lateral relief offers a greater surface for the m. extensor digitorum communis (Guérin, 1980). As this relief is more developed in African species than in Asiatic ones, this may suggest an effect of phylogeny or locomotion in different habitats or both. On the ulna, the developed olecranon process constitutes a powerful lever arm for forearm extensors such as the m. triceps brachii and the m. anconeus, also acting upon the bone for gravitational support. The medial development of the olecranon process is related to larger insertions for the mm. flexor carpi ulnaris, flexor digitorum profundus and flexor digitorum superficialis, all essential to resist hyperextension of the wrist. The cranially reduced anconeal process allows a greater extension of the forearm than in other taxa (e.g., bovids or equids) (Hildebrand, 1974) but prevents a complete verticality of the member as observed in elephants for example (Osborn, 1929). The distal epiphysis shows a reduction of both radial and ulnar styloid processes toward high body mass, adding a mediolateral degree of freedom to the wrist articulation. However, the proximally reduced articular surface for the scaphoid limits the craniocaudal wrist flexion (Yalden, 1971). These morphological traits allow the foot to bear the weight on different substrates while limiting the risk of wrist hyperflexion (Domning, 2002).

Hind limb bones

In the hind limb, the femur expands mainly in the mediolateral direction for rhinos with high body mass and bone size, tending to indicate a stronger resistance to constraints both linked to body propulsion and weight bearing (Lessertisseur & Saban, 1967), exerted in the mediolateral direction (Hildebrand, 1974). The mediolateral reinforcement of the femur is mainly located under the head and the neck, responding to a concomitant enlargement of the medial condyle and epicondyle on the distal epiphysis, both indicating an increase of the body load near the sagittal plane. The more distal location of the lesser trochanter improves the lever arm of the mm. psoas major and iliacus, developing slower but stronger hip flexions (Hildebrand, 1974; Polly, 2007). The same phenomenon is observed with the third trochanter, situated half way along shaft—contrary to in cursorial Perissodactyla like equids, where the third trochanter is more proximally situated (Hermanson & Macfadden, 1996; Holbrook, 2001; Barone, 2010a). However, it has been shown that the relative position of the third trochanter barely varies among extinct rhinoceroses considered as “cursorial” or “semi-cursorial” (Prothero, 2005). This position along the shaft may thus be influenced by both mechanical and phylogenetic constraints. The extreme development of the third trochanter associated with a distolateral development of the greater trochanter also creates a large lever arm for the fascia glutea, the mm. gluteus superficialis and gluteus medius allowing strong hip flexion and abduction. This association appears the greatest for R. unicornis, where the greater and third trochanters can be fused by a bony bridge. Conversely, the greater trochanter is less proximally developed than in related groups like horses and tapirs (Radinsky, 1965; Hermanson & Macfadden, 1996; Holbrook, 2001); as this trochanter is the insertion area for the m. gluteus medius, the main extensor of the hip, the extension in rhinos seems less powerful than in cursorial perissodactyls. On the distal epiphysis, the lateral torsion of the rotation axis of the trochlea in heavy rhinos also indicates a more laterally deviated position of the knee. This conformation may improve weight bearing, shifting the body mass laterally to the body, as previously observed on a study of pressure patterns of the feet in C. simum (Panagiotopoulou, Pataky & Hutchinson, 2019). No real difference in the bone curvature related to body proportion was noticed, confirming previous observations on the independence of femur curvature with regard to body mass increase in quadrupedal mammals (Bertram & Biewener, 1992).

On the hindlimb zeugopodial elements, when the proximal epiphysis of the tibia broadens craniocaudally, the proximal fibular epiphysis is reduced in this direction, despite an increased general robustness. The proximal epiphysis of the fibula is also oriented far more cranially than in lighter specimens. The enlargement of the tibial plateau thus seems to involve a relative reduction in size of the fibular head. The distal epiphyses of both bones covary too, with a broadening mainly expressed in the craniocaudal direction. The medial condyle of the tibial plateau enlarges strongly, resulting into an asymmetrical proximal epiphysis. Moreover, the broadening of the tibial tuberosity correlates with a stronger and larger patellar ligament, reinforcing the knee articulation and therefore the lever arm created by the patella (Hildebrand, 1974). On the distal epiphysis, the two malleoli are more mediolaterally inflated but less distally expanded, allowing the tarsal articulation to move more freely in heavier rhinos (Lessertisseur & Saban, 1967). This trait is associated with a slightly shallower distal articular surface, conferring more important degrees of freedom to the ankle articulation for high body mass (Polly, 2007). This observation is coherent with similar analyses conducted on rhino ankle bones (C. Etienne, 2019, personal communication, showing notably that the talus bone is flattened and has a shallower groove toward high body mass among rhinos).

In addition to the reduction of the proximal epiphysis, the fibula displays a straighter diaphysis for large rhinos as opposed to the greatly curved one for lighter rhinos (see Figs. S6 and S7). This is consistent with previous observations: although the fibula was not considered in their study, Bertram & Biewener (1992) noted a decrease of tibia curvature while body mass increases among terrestrial mammals. In our rhino sample, the tibia shows a very slight straightening of the diaphysis. However, this straightening, perhaps linked to load carrying capacity, appears to be more pronounced on the fibula.