Osteopathology in Rhinocerotidae from 50 Million Years to the Present

Individual elements of many extinct and extant North American rhinocerotids display osteopathologies, particularly exostoses, abnormal textures, and joint margin porosity, that are commonly associated with localized bone trauma. When we evaluated six extinct rhinocerotid species spanning 50 million years (Ma), we found the incidence of osteopathology increases from 28% of all elements of Eocene Hyrachyus eximius to 65–80% of all elements in more derived species. The only extant species in this study, Diceros bicornis, displayed less osteopathologies (50%) than the more derived extinct taxa. To get a finer-grained picture, we scored each fossil for seven pathological indicators on a scale of 1–4. We estimated the average mass of each taxon using M1-3 length and compared mass to average pathological score for each category. We found that with increasing mass, osteopathology also significantly increases. We then ran a phylogenetically-controlled regression analysis using a time-calibrated phylogeny of our study taxa. Mass estimates were found to significantly covary with abnormal foramen shape and abnormal bone textures. This pattern in osteopathological expression may reflect a part of the complex system of adaptations in the Rhinocerotidae over millions of years, where increased mass, cursoriality, and/or increased life span are selected for, to the detriment of long-term bone health. This work has important implications for the future health of hoofed animals and humans alike.


Introduction
Rhinos diverged from their closest living relative, the tapir, about 50.3 million years ago (Ma) [1,2,3] and quickly increased in abundance and species richness through the mid-Cenozoic. The rhinocerotid lineage is hypothesized to have diversified into four major clades in North America and Eurasia: the Diceratheriinae in the Oligocene, and the Aceratheriinae, Teleoceratinae, and Rhinocerotinae in the Miocene [4,5]. Cursoriality, or the habit of running, has been hypothesized to have been maintained through the majority of these lineages [5,6]. These grazing and browsing lineages were some of the most numerous and widespread large mammals in We initially expected to see a correlation between the severity of pathological expression and an increase in rhino mass and cursorial habits, because of the known correlation between osteopathology and mechanical stress [16,18,20,21,22], as well as previous observations of pathologies in rhinos [4,13,14,15]. We reasoned that an increase in mass would put greater stress on bones and joints, increasing the likelihood that arthritis-like pathologies, such as osteophyte formation and articular surface degradation, would occur. If this were the case, the tendency to develop stress-related osteopathologies would be trackable and predictable. We asked two overarching questions in this study: (1) Do these osteopathologic features exhibit a trend over time? And (2) what is the relationship between mass and osteopathology?

Materials and Methods
To determine the relationship between mass and osteopathology through time, we collected data on osteopathologies from a number of extinct and extant taxa in the family Rhinocerotidae, Table 1, and an outgroup, Hyrachyus eximius, a perissodactyl sister group to the Rhinocerotidae [8,23]. We collected data from localities with a large number of rhino skeletal elements to avoid individual preservation bias as much as possible, forming a series of species-level "snapshots" of the rhino lineage. Fossil species were chosen to span the temporal range of rhinocerotids and for the presence of adequate samples of identified elements, Table 2. Data resolution is not on the order of populations, but species in formations; this lumping of occurrences allows us to achieve a statistically adequate sample. For example, there are 15 different localities that comprise the Hyrachyus eximius sample, but they are all part of the Bridger Formation within Uinta County, Wyoming. No permits were required for the described study, which complied with all relevant regulations.
Phylogenetic Data, R scripts and digital photographs associated with this study are available at Morphobank (project ID: 1238) [27] with permission from The American Museum of Natural History, The University of Washington Burke Museum, The University of Oregon Museum of Natural and Cultural History, and The University of Texas Jackson School of Geosciences Vertebrate Paleontology Laboratory. Digital photographs of fossils from the UCMP are also available through the Calphotos archive (http://calphotos.berkeley.edu/). S1 File of raw pathology scores available online through PLoS One.

Study Species
Hyrachyus eximius, (50. 5-45.4 Ma) (Fig 1) is a sister lineage of both the tapir and rhinos. This species is estimated to have weighed around 36.3 kg (equivalent in mass to a large dog), lacked           horns, and was a cursorial browser [11,23]. The rhinocerotid with the earliest first appearance datum (FAD) included in this study is the basal rhinocerotid Trigonias osborni. T. osborni also lacked horns and was substantially larger than H. eximius at about 677 kg [24,25]. T. osborni is known from the Chadronian (37.2 to 33.9 Ma) and was a cursorial browser. Menoceras arikarense emigrated from Europe in the late Oligocene or early Miocene (24.8-20.43 Ma) and had a mass around 375 kg [24,25]. M. arikarense is notable for two firsts: horns and grazing [5]. Diceratherium niobrarense is larger than M. arikarense (about 1010 kg [5]). Although D. niobrarense also displays laterally paired rostral horns, it is thought be descended from Subhyracodon and is not considered a sister group of M. arikarense [5,18]. This rhinocerotid was present in North America in the early and middle Miocene  and was probably a browser [5]. Both M. arikarense and D. niobrarense show morphologies characteristic of increased graviportality: increased bone robusticity, more vertically-oriented pelvis [26], and widening rib cage [5]. Limb length also decreased relative to mass [5].
Aphelops mutilis and Teleoceras hicksi are similar to modern rhinos in graviportal morphology and robust limbs [5]. A. mutilis was a hornless aceratheriine browser known from the mid-Miocene to the beginning of the Pliocene (10.3-4.9 Ma) and is estimated to have weighed around 1840 kg. T. hicksi (10.3-4.9 Ma) is morphologically similar to aquatic hippos [5], but has highly hypsodont teeth [5] with enamel oxygen isotope ratios similar to terrestrial herbivores [29]. T. hicksi is estimated to have weighed around 1660 kg, is thought to have a small nasal horn and is one of the last rhinocerotids in the North American fossil record [30].
From the five modern taxa we examined in planning this study we chose Diceros bicornis (the black rhino) as the modern exemplar. Diceros (5.3332 Ma to present) weighs 800-1,350 kg [24] and is a browser with a prehensile lip specialized to grab foliage [26].

Data Collection Procedure
Diagnosing specific diseases from osteopathologies (often the only pathologies available for study in fossil taxa) is difficult, but not impossible. Certain recognized diseases and disorders can leave distinctive features; e.g., six fingers in a human skeleton are an indicator of polydactyly. The majority of diseases display a common range of pathologies and it is these unique combinations of pathologies that are most informative. For example, irregular holes in a bone can be caused by abnormal nutrient canals, bone infection, soft-tissue swelling, or preservation damage. Arthritis may cause bones to form these irregular holes as well as bone exostoses or thinning, lipping, and fibrous, candlewax, and lumpy bone textures. Arthritis is often labeled spondylarthropathy in non-human paleopathologic studies [17,31] to acknowledge that arthritis itself is not a specific disease, but can be caused by a range of environmental, genetic, and behavioral factors depending on the system under study [16,18,19]. Each specimen was digitally photographed with a Nikon D90 camera. The camera was hand held approximately perpendicular to the photographic plane. Elongate fossils (e.g. femora or D. bicornis has no blue line because only modern bones were examined. Tree was pruned from Cerdeño's 1998 [10] morphologic phylogeny or Rhinocerotidae and timecalibrated in RStudio using the 'equal' setting in the function timePaleoPhy() in the software package 'Paleotree' [28]. Tree was set to be fully dichotomous and to extend all the way to the LAD. doi:10.1371/journal.pone.0146221.g001 metapodials) were photographed in lateral view and fossils with irregular shapes (i.e. podials) were oriented in medial view. Proximal and distal articular surfaces were photographed as well for limb and foot elements. Vertebral elements were photographed in dorsal, ventral, proximal, and distal views. Extra photos were taken if a unique pathology was observed or for striking examples of specific pathologies.
The specimen number and corresponding photo numbers were recorded digitally and associated with a pathology index scoring and any qualitative observations (see supporting information). The presence or absence of pathology was recorded on site, while pathology severity scorings were determined from the digital photographs.
We quantitatively described the visible surface of each bone using a category, or binning, system. We sorted our initial qualitative descriptions of possible symptoms of disease into seven different categories. The seven categories were divided into ranks from 1 (regular bone) to 4 (severe) (Fig 2). These ranks are artificial, but should allow for consistent scoring. All scoring was completed by the first author. The seven categories are exostoses, lipping, bone texture, cavitation, foramen shape, foramen size, and articular surface modification. All categories except for 'articular surface' refer to the nonarticular surfaces of the bones. The categories were chosen following the methodology of Aufderheide [16], Rothschild [17,31], and Bartosiewicz et al. [21]. Analogs of this procedure have been used for decades in anthropologic [21,22,32] and modern cattle [21] studies. Our goal was to quantitatively describe all irregularities observed in the osteology of the Rhinocerotidae, even if they could not immediately be categorized as a pathology.
Category One: Exostoses. Exostoses are formations of new bone on the surface of a bone, caused by inflammation of the periosteum. Extoses appear as bumps or protuberances on an area of the bone that is expected to be smooth or relatively flat. This category includes ossification of the periosteum, ligaments, or muscle. Bones in rank one do not exhibit any exostoses. Bones in rank two show minor irregular bulging of bone. Bones in rank three show clear protrusions of irregular bone. Bones in rank four show a continuous irregular distortion of the non-articular surface of the bone.
Category Two: Lipping. Lipping occurs when osteophytes (commonly referred to as bone spurs) form as new bone on the margin of articular surfaces. They usually form as a series of merging osteophytes around the joint margin, but can occur singly as well. Bones in rank one do not exhibit any lipping. Bones in rank two show slight bulging of the bone adjacent to the articular surface. Bones in rank three show bulging of the bone surrounding the articular surface to the point where a prominent shelf is beginning to form. Bones in rank four show a prominent shelf adjacent to the articular surface. The shelf may be regular or irregular.
Category Three: Textures. Bone constantly remodels and rebuilds itself in response to localized stress. This can result in characteristic external textures. Care must be taken to not conflate exostoses (which has more to do with shape) with texture. Bones in rank one have a smooth texture. Bones in rank two have an elevated linear texture, termed fibrous. Bones in rank three have an elevated linear texture that is slightly bulging or uneven texture, likened to candle wax. Bones in rank four have an elevated, uneven, nonlinear texture.
Category Four: Cavitation. Cavitation is the first category concerned with loss of bone. A cavitation is a hole in the bone, usually caused by infection and/or decreased blood flow. Unlike the categories of foramen shape and size, these are relatively large areas of the bone that cannot be confused with vascularization. Bones in rank one do not exhibit any cavitation. Bones in rank two show a pockmarked appearance where the bone has lost integrity. Bones in rank three show small cavities. Bones in rank four show large cavities that may be linked together.
Category Five: Articular Surface. The articular surface forms the bony portion of a joint. Bones in rank one do not show any irregularities in the joint surface. Bones in rank two show a pockmarked appearance where the cartilage has been worn away. Bones in rank three show bone loss on the articular surface. Bones in rank four show eburnation of the articular surface and/or osteophyte formation. Category Six: Foramen Shape. We found bone cysts can easily be confounded with vascularization (called 'lucencies' in ), so we decided to describe the degree of foramen deformation instead of labeling all foramina as cysts. Cysting (pockets or holes where localized infections occurred [16,18]) was divided into two categories (foramen shape and size). Rank one consists of circular foramina on the surface of the bone. Rank two consists of elongate or ovoid foramina. Rank three consists of elongate foramina that are irregularly ovoid, but still linear. Rank four consists of irregular, nonlinear (or bent) foramina.
Category Seven: Foramen Size. Rank one consists of foramina of approximately the same size. Rank two consists of foramina which show little variation in size relative to one another. Rank three consists of foramina which show moderate variation in size relative to one another. Rank four consists of foramina which show a high degrees of variation in size relative to one another.
Each bone was also classified as appendicular and axial. To explore whether there were any overt patters of regionalization in pathological expression, all appendicular elements were then divided into the functional categories: hindlimb or forelimb, and also developmental categories: girdle, stylopod, autopod, zeugopod. The overall percent expression of each category was tabulated and then compared relative to the total number of appendicular element.

Data Analysis
Each fossil was given a score of 1-4 for each pathology, and these scores were then averaged for each taxon and pathological category, yielding 49 results. These averaged scores were then added together for each taxon (i.e. the seven pathology categories were added together for each taxon) to create an index of pathology (IPa). The minimum score possible would therefore be seven (all pathological categories in a taxon having a score of one) and the maximum would be twentyeight (all pathological categories in a taxon having a score of four). These average scores for each taxon do not behave as ordinal data, because they are subject to the central limit of means. That is, species averages of non-continuous data behave like continuous data, especially with large sample sizes. The smallest number of specimens we analyzed for one species was 65, more than adequate to produce this effect. Consequently, we decided it was appropriate to analyze these values using continuous-data approaches: linear regression and independent contrasts.
We tested whether mass was associated with increased osteopathologic expression in two ways. First, we ran a series of linear regressions in JMP [33] with estimated mass against each individual categorical score as well as the total index. Mass estimates for extinct taxa were calculated using the total molar length (M1-3) [34] from Radinsky 1967 for Hyrachyus eximius [23] and Prothero 2005 [5] for all other extinct taxa, which we found to be the most reliable of available body mass estimators. Other available proxies (femur length and humerus width) produced unreasonable mass estimates [35] for one or more of the included taxa, likely as a result of the changing degree of graviportality through the history of the rhino lineage.
The second test used Felsenstein's [36] independent contrast (IC) method to examine the influence of shared ancestry on the relationship between mass and pathology. We constructed a fully resolved tree of just the taxa in our study by paring down the results of Cerdeño 1995 [8]. The tree was time-calibrated in RStudio [37] using the packages 'ape' [38] and 'paleotree' [28] with paleotree's function TimePaleoPhy. The r code is available in the S2 File. We used the 'Equal' method within TimePaleoPhy, which prevents zero-length branches, and the setting 'add.term = TRUE', which gave us branch lengths that took LAD into consideration. FAD and LAD for the Equal method were determined by the temporal extent of the formation at the locality where the fossils were excavated. To implement the IC method, we used the package 'ape' [38], to calculate the absolute values of the difference for each pair of nodes for both mass and all seven types of pathologies, as well as for the overall IPa, under a Brownian Motion model. The resulting contrasts for pathologic values were regressed against the contrast for the mass values. The r-squared and p-values for the non-phylogenetic linear regression versus the IC regression analysis were then compared.

Results
Overall, geologically older taxa show the smallest relative abundance of pathologic elements, while the greatest pathologic expression was seen in the more derived taxa, which were also the most massive taxa sampled in North America. The one exception is the extant species, D. bicornis, which appears later, yet is less massive and less pathologic overall than A. mutilis and T. hicksi.
When taxa are considered separately, H. eximius displayed low osteopathologic expression (~28%), most of which was expressed as cysting and exostoses in the podials. T. osborni and D. niobrarense also displayed a greater degree of pathologic expression in the distal elements. The M. arikarense fossil assemblage displayed prominent exostoses. In smaller elements (i.e. podials) the non-articular surfaces would be almost entirely composed of exostoses. A. mutilis and T. hicksi commonly contained large visible cysts and rank three, candlewax, bone texture. Only two fossils displayed eburnation, in A. mutilis on the articular surfaces of a proximal tibia (UCMP F-30266) and in T. hicksi on a distal humerus (UO F-2772). A. mutilis also had the highest percent expression of any one pathology (in this case, foramen shape), while D. bicornis had comparatively more foramen variation adjacent to the articular surfaces and fewer exostoses than the other robust taxa. One specimen (VPL M-8259) had flat 'rice grain' crystals on the proximal articular surface of right radius and ulna, as well as the distal articular surface of the humerus, a possible indication of gout [16,18]. Tendon ossification was only seen in A. mutilis and T. hicksi. Of note, most of the articular surfaces of the synovial joints in both extinct and extant taxa appeared smooth and free of damage.
Overall index of pathology (IPa) scores were between 8 and 18, Table 3. The two oldest lineages, H. eximius and T. osborni, had an overall pathologic score of 8.8 and 11.06 respectively. The next oldest lineage, D. niobrarense, had an overall score of 12.81, while M. arikarense had an overall score of 13.31. T. hicksi had an overall score of 14.26, while A. mutilis had an overall score of 17.57. The modern rhino, D. bicornis, had an overall IPa of 12.23.
When we regressed mass for each taxon against each of the seven osteopathology categories Table 4, four of the pathologies (exostoses, abnormal textures, foramen shape, and foramen size variation) had p-values less than or equal to 0.05. A linear correlation of mass against the overall pathologic scores was found to be significant (p = 0.04) and accounted for about 52% of the variation (r 2 adj.). The IC analysis comparing mass and the seven osteopathology categories was also significant (p 0.05) for the both foramen shape and the overall pathologic index, with mass accounting for 42% of the overall variation.
We were also interested in testing whether certain bones or regions of the appendicular skeleton (which comprises the majority of the data) displayed a greater amount of pathology than other bones or regions of the appendicular skeleton. We divided all appendicular elements into the functional categories: hindlimb or forelimb (Fig 3), and also developmental categories: girdle, stylopod, zeugopod, autopod (Fig 4). For example, if the stresses generating the osteopathology were greater in the distal parts of the limb, one might expect greater pathology in the autopod (manus and pes) than the stylopod (humerus and femur). We found no significant difference in pathological expression between different regions of the appendicular skeleton.

Discussion
In our study we found that mass can explain roughly 50% of the osteopathological expression. A. mutilis, surprisingly, had the highest pathology scores by a wide margin, while T. hicksi, which was close to A. mutilis in estimated mass, had scores similar to the smaller D. niobrarense and M. arikarense. Both the overall expression of pathology and the subcategory of foramen shape were significant when regressed against mass regardless of whether phylogeny was taken into account or not. However, the r 2 value in the vicinity of 0.5 suggests that other factors besides mass, such as bone robusticity, cursoriality or environment, could play a significant role in pathological expression. There might be a tradeoff between a lineage increasing in size or weight and abnormal textures, lipping, and other pathologies that intuitively should be selected against on an evolutionary scale. Lower levels of expression (categories 1 and 2) were more common, but no taxon was entirely pathology-free. The 'maximum operational level' of different pathologies may be similar or drastically different in different vertebrate lineages, which in turn could lead to diverse selection pressures. Longevity could also be a factor in pathological expression. There is a positive correlation in the Mammalia between body mass and lifespan, although there is a great amount of variation [39]. The larger taxa may be living longer, which could increase the likelihood of osteoarthritis, synovitis, traumatic injury, etc. Captive mammals that live longer than their wild counterparts often display these pathologies [4,13,14,15]. Pathology could be a reflection of ontogeny. However, a longer lifespan does not necessarily increase only the geriatric portion of a mammal's life. In an animal with a longer lifespan bone would presumably stay healthy for the same proportion of a mammal's life as the shorter-lived counterpart, but this remains to be tested [40] Our main difficulty in this study was to establish a measurement method for pathology. In anthropology several qualitative and quantitative metrics have been used to study paleopathology [16,18,19]; paleontology also has no universal methodology for identifying and analyzing paleopathologies, but several parallel methods [12,13,14,17,21,27]. Our method uses a scoring system that is focused on parsing out the severity of symptoms, not direct diagnoses of disease. It is possible to apply these separate pathology categories to studies across the vertebrate kingdom. Comparison of pathological expression between these vastly different taxa could lead to new insights into bone repair, species and lineage-level responses to pathology, and the uniformity of bone-related diseases over time.
Synovitis, not arthritis, may be the proximal cause of the pathologies observed. Notably, we found that most of the pathology in the taxa we studied was located immediately adjacent to the articular surface of joints and not in the articular surface itself. That is, the articular surface itself appeared healthy (that is, not scarred or pitted) in all but five individual specimens even in individuals with advanced exostoses and abnormal bone textures. This could indicate that the joints (and therefore the organism) are functional well after pathologies begin to appear and swelling of the synovium caused the observed cortical erosion [16].
The overall picture painted by our results shows a measureable increase in the percentage of elements that display osteopathologies related to arthritis from the older to newer branches in the North American rhinocerotid lineage, consistent with earlier observations [4,13,14,15,16,17]. Our initial hypothesis was that more massive rhino species would display a greater frequency of osteopathologies as increased loading pressure on the limbs caused microfractures that resulted in inflammation and abnormal bone textures. With our current sample, we found a significant correlation of pathology with mass, suggesting that increasing size in the lineage was partially responsible for osteopathologies. This study is not powerful enough to conclude if these pathologies are related primarily to arthritis or is a multifactorial response or a range of diseases. Because size increase is a common adaptation in terrestrial herbivores for both eating low-quality diets [26,41,42] and resisting predation pressure in open environments [26,43], it seems that adaptations for food and predator avoidance may incur a cost in bone stress and osteopathology. The accumulation of pathologies in a lineage may no longer be solely a herald of disaster, but of adaptation as well.