RESEARCH ARTICLE Why ruminating ungulates chew sloppily: Biomechanics discern a phylogenetic pattern Zupeng Zhou1, Daniela E. Winkler2, Josep Fortuny ID3,4, Thomas M. Kaiser5, Jordi Marcé- Nogué ID3,5* 1 School of Mechanical and Electrical Engineering, Guilin University of Electronic Technology, Guilin, China, 2 Johannes Gutenberg University Mainz, Institute of Geosciences, Mainz, Germany, 3 Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, c/ Columnes s/n, Barcelona, Spain, 4 Centre de Recherches en Paléontologie de Paris, Muséum National d’Histoire Naturelle, Bâtiment de Paléontologie, Paris, France, 5 Centrum für Naturkunde, University of Hamburg, Hamburg, a1111111111 Germany a1111111111 a1111111111 * Jordi.marce.nogue@uni-hamburg.de a1111111111 a1111111111 Abstract There is considerable debate regarding whether mandibular morphology in ungulates pri- marily reflects phylogenetic affinities or adaptation to specific diet. In an effort to help resolve OPEN ACCESS this debate, we use three-dimensional finite element analysis (FEA) to assess the bio- Citation: Zhou Z, Winkler DE, Fortuny J, Kaiser TM, mechanical performance of mandibles in eleven ungulate taxa with well-established but dis- Marcé-Nogué J (2019) Why ruminating ungulates chew sloppily: Biomechanics discern a tinct dietary preferences. We found notable differences in the magnitude and the distribution phylogenetic pattern. PLoS ONE 14(4): e0214510. of von Mises stress between Artiodactyla and Perissodactyla, with the latter displaying https://doi.org/10.1371/journal.pone.0214510 lower overall stress values. Additionally, within the order Artiodactyla the suborders Rumi- Editor: Carlo Meloro, Liverpool John Moores nantia and Tylopoda showed further distinctive stress patterns. Our data suggest that a University, UNITED KINGDOM strong phylogenetic signal can be detected in biomechanical performance of the ungulate Received: August 28, 2018 mandible. In general, Perissodactyla have stiffer mandibles than Artiodactyla. This differ- Accepted: March 15, 2019 ence is more evident between Perissodactyla and ruminant species. Perissodactyla likely rely more heavily on thoroughly chewing their food upon initial ingestion, which demands Published: April 17, 2019 higher bite forces and greater stress resistance, while ruminants shift comminution to a later Copyright: © 2019 Zhou et al. This is an open state (rumination) where less mechanical effort is required by the jaw to obtain sufficient dis- access article distributed under the terms of the Creative Commons Attribution License, which integration. We therefore suggest that ruminants can afford to chew sloppily regardless of permits unrestricted use, distribution, and ingesta, while hindgut fermenters cannot. Additionally, our data support a secondary degree reproduction in any medium, provided the original of adaptation towards specific diet. We find that mandibular morphologies reflect the masti- author and source are credited. catory demands of specific ingesta within the orders Artiodactyla and Perissodactyla. Of Data Availability Statement: All of the raw and particular note, stress patterns in the white rhinoceros (C. simum) look more like those of a processed CT-scan data of: Aepyceros melampus, general grazer than like other rhinoceros’ taxa. Similarly, the camelids (Tylopoda) appear to Alcelaphus buselaphus, Ceratotherium simum, Dicerorhinus sumatrensis, Diceros bicornis, Equus occupy an intermediate position in the stress patterns, which reflects the more ancestral quagga and Litocranius walleri, has been uploaded ruminating system of the Tylopoda. to MorphoSource.org at the following link: https:// www.morphosource.org/Detail/ProjectDetail/Show/ project_id/730. Funding: J.F. is supported by the postdoc grant “Beatriu de Pinos” 2014 – BP-A 00048 from the Generalitat de Catalunya, and by the Spanish Ministerio de Economı́a, Industria y Competitividad PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 1 / 21 Why ruminating ungulates chew sloppily and the European Regional Development Fund of Introduction the European Union (MINECO/FEDER EU, project CGL2017-82654-P) and by the Catalonian The main function of the mammalian mandible is to transmit the forces generated by mastica- Research Group 2017 SGR 86. D.E.W. was funded tory muscles to the teeth for food processing [1]. The shape of the mandible has been proposed by the European Research Council (ERC) under the to reflect the frequency and magnitude of muscle activation, bite force, and the properties of European Union’s Horizon 2020 research and the food ingested, while ensuring structural integrity under all loads without dissipation or fail- innovation programme (ERC CoG grant agreement ure [2]. In general, the greater the force required to fracture ingesta, and the more repeatedly no. 681450 to T. Tütken). J.M-N. was supported by the Deutsche Forschungsgemeinschaft (DFG, such forces need to be produced, the stiffer the mandible has to be to maintain its structural German Research Foundation, KA 1525/9-2). This integrity [3]. For example, it has been demonstrated that among primates, hard-food consum- research is publication no. 93 of the DFG Research ers have stiffer mandibles compared to soft-food consumers [4]. Unit 771 ‘Function and performance enhancement The relationship between bite force, size and shape of the mandible in relation to ingesta in the mammalian dentition—phylogenetic and has been studied previously in several mammalian clades including Artiodactyla [5–7], Chir- ontogenetic impact on the masticatory apparatus’. This research paper is a contribution to the CERCA optera [8–12], Primates [13–15] and Carnivora [16–19]. Due to the close interaction between programme (Generalitat de Catalunya). the mammalian feeding mechanism and diet, the biomechanical study of extant species can illuminate ecomorphological adaptations to provide accurate predictions of extant species and Competing interests: The authors declare that they have no competing interests. potentially acquire valuable tools for the reconstruction of oral behaviour in extinct taxa [20– 22]. In ungulates, the hypsodonty index (HI), among other morphological trait, and, to a lesser extent, the moment arms of jaw muscles have been used to infer dietary traits in extinct species [7,23]. The hypsodonty index relates to ingesta abrasiveness and clearly increase in grazers due to higher dietary abrasive component [24]. For example, Varela and Fariña [7] demonstrated a relationship between the HI, masseter moment arm, and the amount of grass in the habitat. It is further well established that grazers tend to have more robust mandibles and larger masse- ters compared to browsing species, because they have to comminute the relatively tough, highly fibrous ingesta more thoroughly [25]. Mendoza and Palmqvist [26] demonstrated that the length of the mandible and the HI are higher in ungulates from open (hence dustier) habi- tats, whereas a wider muzzle is common in bulk, non-selective grass feeders. These observa- tions suggest that mandible shape variation between species can be linked to a different mode of feeding especially when comparing specialized browsers with grazers. However, mandible movement types are also constrained by developmental processes so that interspecific varia- tion might reveal strong phylogenetic signal in morphological design and performance [6,17,27–29]. This phylogenetic signal is found to be stronger in Perissodactyla than Artiodac- tyla [30]. The degree to which mandibular morphology is influenced by diet, reflects phyloge- netic affinities or both is thus still debated. Biomechanical data, which better characterizes ungulate mandibular diversity, may shed more light on this relationship and help to disentan- gle contradictory interpretations. Finite element analysis (FEA) is an ideal technique to provide a comparative perspective on mandibular biomechanics. Recently, FEA was used to assess stress during mastication in the mandibles of different vertebrate taxa using 3D models [31–34] and plane models [35–38]. In a comparative context, FEA studies have employed a variety of approaches including qualita- tive and quantitative approaches. Qualitative approaches include stress distribution plots, while quantitative methods involve examining stress at particular points, comparing model means and combining FEA with geometric morphometrics to compare model distortion [39]. Both, quantifying stress data at particular points [38,40–46] as well as taking into account the strength of the whole model by computing von Mises stress averages, have proven useful in ecomorphological analyses [37,47–54]. FEA outputs have been used to explore functional mor- phology, ecomorphology and macroevolution, by applying standard statistical methods [35,55] or geometric morphometrics to analyse deformations after load application [56]. PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 2 / 21 Why ruminating ungulates chew sloppily In this study, we employed “weight-meshed values” and “quasi-ideal meshes” to compute the values of stress recently proposed by Marcé-Nogué et al. [4,34,57]. This method is more robust compared to other comparative approaches because it considers the effects caused by the different size and homogeneity of the mesh elements in the FEA model. We also employed the new intervals’ method to analyse data from FEA models in a comparative multivariate framework [58,59]. This method allows for considerably more effective comparison of FEA models, and consequently more precise distinction between feeding strategies. Both, the quan- tification of average values of stress [57] and the use of the intervals’ method [58], have been previously employed to discriminate dietary preferences in the study of the biomechanical traits of armadillo mandibles by separating the results of stress as a function of diet. In the present study we compared 3D finite element models of mandibles of Artiodactyla (Alcelaphus buselaphus, Aepyceros melampus, Camelus dromedarius Giraffa Camelopardalis, Lama glama and Litocranius walleri) and Perissodactyla (Ceratotherium simum, Diceros bicor- nis, Dicerorhinus sumatrensis, Equus quagga and Tapirus terrestris) selected because they cover different examples of proportions of grass in their diets [60,61]. We analysed lateral and orthal biting action and explore the degree to which morphology reflects adaptation towards diet or underlying phylogenetic relationships. Materials and methods Specimens A total of 11 extant ungulate taxa were analysed. All specimens examined are housed either at the Centrum für Naturkunde in Hamburg (ZMH), Germany, the Museum für Naturkunde Berlin (ZMB), Germany, the Museum Victoria in Melbourne (NMV), Australia or in the Texas Memorial Museum (TMM), USA. All the specimens are adult individuals with no pathology. See Table 1 for further details. Species were classified according to dietary trait as follows [60,62]: GGR: general grazer, GBR: general browser, OMF: open-habitat mixed-feeder and HBR: high-level browser (Table 1). Geometric reconstruction High resolution CT scans of Alcelaphus buselaphus, Aepyceros melampus, Equus quagga and Litocranius walleri were obtained at Steinmann-Institut für Geologie, Mineralogie und Paläon- tologie (Universität Bonn, Germany) using a CT scanner (Philips Brilliance 64). Specimen dimensions ranged between approximately 20–40 cm in length resulting in voxel sizes between 0.456 mm and 0.793 mm. Scans of Diceros bicornis and Ceratotherium simum were performed at Universitätsklinikum Hamburg Eppendorf (UKE), Germany, on a medical Philips Brilliance 64 CT scanner (voxel size 0.977 mm and 0.651 mm respectively). Lama glama (voxel size 0.1768 mm) and Tapirus terrestris (voxel size 0.5918 mm) were scanned on a medical CT at The University of Texas High-Resolution X-ray CT Facility (UTCT) (see www.digimorph.org for further details). Giraffa camelopardalis (voxel size of 0.976 mm) and Camelus dromedarius (voxel size 0.804 mm) were scanned in a medical CT (Siemens Sensation 64 scanner) at St. Vincent’s Hospital, Melbourne, Australia, and Dicerorhinus sumatrensis was scanned in a medical CT at Yxlon facilities in Hamburg, Germany using a Y.CT Modular (YXLON Interna- tional GmbH) equipped with a Y.TU450-D11 x-ray tube (voxel size of 0.23 mm). CT datasets were imported to the software Avizo 7.0 (FEI-VSG Company). The scans were manually segmented and a surface model was reconstructed with multiple materials. Molar teeth (M1, M2 and M3) were reconstructed as a separate material from the mandibular bone, as they bear the highest functional load in comminution [63]. Molars were then included in the models as auxiliary geometry to apply the force in the bite point. During this step, PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 3 / 21 Why ruminating ungulates chew sloppily Table 1. List of ungulate species investigated. SPECIES ID ORDER DIET AGE/SEX STATUS Alcelaphus buselaphus ZMH-S-7487 Artiodactyla GGR AF Captive Aepyceros melampus ZMH-S-10162 Artiodactyla OMF YAM Captive Camelus dromedarius NMV R 5444 Artiodactyla OMF AM Captive Ceratotherium simum ZMH-S-2552 Perissodactyla GGR SAF Wild Diceros bicornis ZMH-S-9379 Perissodactyla GBR Adult-U Wild Dicerorhinus sumatrensis ZMH-S-8273 Perissodactyla GBR Adult-U Wild Equus quagga ZMB-Mam-70335 Perissodactyla GGR U Wild Giraffa camelopardalis NMV C 32820 Artiodactyla HBR Adult-U Captive Lama glama TMM M-2052 Artiodactyla OMF Female U Litocranius walleri ZMB-Mam-39663 Artiodactyla HBR U Wild Tapirus terrestris TMM M-16 Perissodactyla GBR Subadult-U U Museum Acronyms: ZMH = Centrum für Naturkunde, Hamburg, Germany, ZMB = Museum für Naturkunde, Berlin, Germany, TMM = Texas Memorial Museum, Austin, USA and NMV = Museum Victoria, Melbourne, Australia); Diet classification after Mendoza et al. [60] GGR: general grazer, GBR: general browser, OMF: open-habitat mixed-feeder and HBR: high-level browser. YAM: young adult male, SAF: Subadult female, AF: Adult female, AM: Adult male, U: Unknown. https://doi.org/10.1371/journal.pone.0214510.t001 irregularities in the surface resulting from model generation during scanning were repaired using Rhinoceros 5.0 software (McNeel & associates) refinement and smoothing tools. The segmented models were converted to a CAD models following Marcé-Nogué et al. [64]. A CAD model is a vector-graphic standard for three-dimensional design used in mechanical design and engineering which does not involve any tetrahedralization of the geometry such as stereolithographic formats. The CAD-step allowed us to create the mesh of finite elements in ANSYS FEA Package. All models were standardized and oriented in anatomical position with maxilla and mandible articulated to facilitate comparison. Model properties A structural static analysis was performed to evaluate the biomechanical performance of the mandible during biting using the finite element package ANSYS 17.1 (Ansys, Canonsburg, USA) and a Dell Precision Workstation T7910 with 256GB. Our primary interest was the com- parison of von Mises stress distribution throughout the mandible under loading conditions defined by masticatory movement. The von Mises criterion is the most accurate value for pre- dicting fracture location when isotropic material properties are applied to cortical bone [65]. For this study, elastic, linear and homogeneous material properties were assumed for bone tis- sue. Although the specific mechanical properties of bone for the investigated individuals are unknown, Gil et al. [66] demonstrated that the use of general values for bone does not alter results to an extent that a relative comparison between models is inappropriate. Similarly, Strait et al. [67] found that models with heterogeneous mechanical properties of bone closely matched models assuming homogenous properties. Therefore, values from bovine haversian bone: E [Young´s modulus] = 10 GPa and v [Poisson ratio] = 0.4)[68] and for tooth enamel E = 50.14 GPa and v = 0.3 [69] were used and homogenous mechanical properties were assumed for all bones in this study. The high stiffness of the molars compared with the bone makes the molars act as a rigid body that helps the transmission of the forces to the mandible. The mandibles were meshed using the ANSYS mesh module with an adaptive mesh of hexahe- dral elements [70]. The mesh of each model was approximately between 0.8–1.5 million elements. PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 4 / 21 Why ruminating ungulates chew sloppily Boundary and loading conditions Two bite configurations: lateral and orthal, were simulated based on configurations described by Fortuny et al. [71] and Fletcher et al. [36], respectively (Fig 1). Bite positions were placed on the first, the second, and the third molar separately, generating six cases per individual. The lateral bite configuration was defined using extrinsic loading and the application of force in the z-direction in the plane of the molar row, resulting in a lateral chewing motion [72]. It needs to be noted that postcanine orthal biting is a less significant component of masticatory movements in herbivorous ungulates where lateral power stroke component predominate in comminution. Orthal biting action however would apply on ingestion when using incisors. Previous works on ungulates [36] tested digestive physiology versus mandible stiffness in orthal bite but did not investigate lateral movements, although lateral occlusal movements are the basic chewing strokes in ungulates. For the orthal bite configuration, boundary and loading conditions were defined using intrinsic loading and including only the masseter and temporalis muscles. The insertion areas of both muscles were defined in the model in order to apply the muscle forces during biting. The direction of these forces was defined by the line that joins the centroid of the area of origin on the cranium with its attachment point on the mandible. Because the force that a muscle can produce depends on its active cross section, the size of the attachment area for each muscle was used to define the application of muscle tension at the point of attachment. However, as the interest of the study is a comparative one, instead of calculating absolute values for forces, we used these areas to quantify the proportional amount of force developed by each muscle assuming a muscular pressure of 0.3 MPa [73]. To simulate a bite, a fixed displacement bound- ary condition in y-direction was applied in the bite locations to model the moment at which cranium and mandible contact the food. The six bite cases simulated were: 1. Lateral first molar biting, 2. Lateral second molar bit- ing, 3. Lateral third molar biting, 4. Orthal first molar biting 5. Orthal second molar biting and, 6. Orthal third molar biting. In all cases, restrictions in displacement were applied at the occipital condyles to fix the mandible in the x and y-direction in relation to the vertebral column (Fig 1). Muscle force: Scaling the models To compare the performances of structures that differ in shape and size, the values of muscular contraction transmitted to bone were calculated according to the methodology developed by Fig 1. Free-body diagram of the biomechanical problem with boundary conditions, muscular forces, area of insertion, and bite position for each molar (M1, M2, M3) in L. walleri for a) lateral and b) orthal biting configurations. https://doi.org/10.1371/journal.pone.0214510.g001 PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 5 / 21 Why ruminating ungulates chew sloppily Marcé-Nogué et al. [74] and adjusted for 3D models by Fortuny et al. [75]. This methodology uses allometric principles to relate mandibular volume for each specimen with the muscular forces applied by a 2/3 power relationship and agrees with the allometric proportions of the species [41,76,77]. To obtain calibration values for force (lateral bite configuration) and muscular pressure (orthal bite configuration), Litocranius walleri was employed as a reference model, with a value of 0.3 MPa assumed for muscular contraction pressure [73]. Models for each of the remaining species were scaled against L. walleri according to their respective difference in mandibular volume. According to the equation proposed by Fortuny et al. [75] muscular pressure (P) is defined as force divided by the area of muscular insertion (MI) for each muscle. The muscular pressure of both models A and B are related to variation in the volume (V) of the mandible (Eq 1 and Table 2). sffiffiffiffiffi !2 MIB 3 VA PA ¼ PB ðEq 1Þ MIA VB Average values and quasi-ideal mesh The distribution of von Mises stress is useful in biomechanics because von Mises stress pat- terns indicate relative strength. For instance, specimens that are characterized by higher values of stress are considered less resistant to the effect of the forces. In this study, we used von Mises stress maps of the mandibles, the average values of von Mises stress and box-plots of the stress distribution in each taxon to investigate mandibular strength. Average values of von Mises stress were evaluated following the method proposed by Farke [47], who recommends plotting stress distributions as quantitative data. However, the use of box-plots for stress distri- butions and the statistics derived from them (e.g. percentiles or whiskers) requires the use of a quasi-ideal mesh (QIM) which necessitates corrections for the non-uniformity of the mesh. For non-uniform meshes (where different elements have different sizes), new statistics, which take into account non-uniformity such as the mesh-weighted arithmetic mean (MWAM) [57] have been proposed. The phylogenetic signal was estimated for the biomechanical data using a Table 2. Geometric properties and muscular pressure. Volume of the mandible [mm3], masseter and temporalis area [mm2] and muscular pressure [MPa] applied in the FEA model for each species. Model Volume of the mandible Lateral Bite Cases 1,2 Orthal Bite Cases 4, 5 and 6 [mm3] and 3 Lateral force [N] Masseter Area Temporalis Area Masseter Pressure Temporalis Pressure [mm2] [mm2] [MPa] [MPa] A. buselaphus 372240 204.91 5334 1932 0.3226 0.3333 A. melampus 114826 93.55 1551 510 0.5063 0.5768 C. dromedarius 967170 387.28 8366 1542 0.3887 0.7894 C. simum 5275100 1199.98 26769 3137 0.3764 1.2026 D. bicornis 3266000 871.70 12406 1972 0.5900 1.3897 D. sumatrensis 1531700 526.18 10254 1169 0.4309 1.4146 E. quagga 802150 341.87 12671 1059 0.2265 1.0151 G. 1133000 430.37 10824 2763 0.3339 0.4896 camelopardalis L. glama 198380 134.69 3104 959 0.3643 0.4416 L. walleri 20852 30.00 841 314 0.3000 0.3000 T. terrestris 823910 348.02 9414 2883 0.3104 0.3795 https://doi.org/10.1371/journal.pone.0214510.t002 PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 6 / 21 Why ruminating ungulates chew sloppily mathematical generalization of the K-statistic appropriate for multivariate data (i.e. Kmult) [78]. For the mandibles, the MWAM, and the stress quartile values given in the boxplots, the per- centage error of the arithmetic mean (PEofAM) and the percentage error of the median (PEofM) (described in Marcé-Nogué et al. [57]), as well as the values needed to calculate both percentage errors, can be found in Supplementary information (S1 and S2 Tables) All analyses were performed in Matlab R2017a (Mathworks, Massachusetts, USA). Cluster analysis Cluster analyses (Ward’s method for agglomerative-hierarchical analysis) were performed using the Von Mises stress distribution (MWAM and percentiles) data. Two different cluster analyses were performed: one using the average MWAM values and the second using the per- centiles of the six biting cases. Euclidean distances were used as the similarity index [79], cophenetic correlation coefficient to show the stability of the cluster [80] and a dendogram was generated using the results to determine affinities between taxa in the biomechanical data. The data was scaled and centred to improve numerical stability. All analyses were performed in R v. 3.4.4. [81]. Intervals’ method A new variable hereafter referred to as the “vector for stress intervals” representing a different interval of stress values was created following the interval method [58]. The interval method involves a new organization of the stress values of the models subdividing stress values obtained into defined intervals, with each interval representing an subarea (as a percentage) of the original model and reflecting a specific range of stress. The interval method described by Marcé-Nogué et al. [58] requires the definition of a fixed upper threshold FTupper = 15 MPa for the three lateral bite cases and a fixed upper threshold FTupper = 25 MPa for the three orthal bite cases as well as a specific number of intervals (N = 100). The number of intervals for our analyses were determined using the convergence procedure described by Marcé-Nogué et al.[58] (S1, S2, S3, S4, S5 and S6 Figs and S3 Table of the supplementary information). Although fewer intervals could have been be used for the orthal cases, we decided to standardize all the analysis to N = 100, which is suitable for both orthal and lateral cases. The values of each interval when N = 100 for each specimen and for each biting case can be found in the supplementary information (S1 Document). Subsequently, these new vectors for stress intervals were analysed using multivariate meth- ods such as principal component analysis (PCA) and plotted in a biplot. All analyses were per- formed using Matlab (Mathworks, Massachusetts, USA) and R v. 3.4.4. [81]. The phylogenetic signal was also estimated for the intervals’ method data using a mathematical generalization of the K-statistic appropriate for multivariate data. Results Von mises stresses A general pattern distinguishes the orders Artiodactyla and Perissodactyla. Artiodactyla man- dibles are characterized by more areas of high stress in the ramus whereas Perissodactyla exhibit lower stress values, particularly across large areas of the corpus (Fig 2). The three lateral bite cases 1, 2 and 3 showed high stresses in the ramus with the highest stress occurring in the condyle (Fig 2). The orthal bite cases 4, 5 and 6 exhibited an important area of lower stress in the region between the bite position (molars) and the mandibular symphysis. All mandibles PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 7 / 21 Why ruminating ungulates chew sloppily Fig 2. Phylogenetic tree and modelled von mises stress distribution in the mandible and the molars for each taxon for each bite configuration and bite position. Lateral cases 1) in M1 2) in M2 and 3) in M3. Orthal cases 4) in M1 5). In M2 and 6) in M3. Phylogenetic tree from http://10ktrees.nunn-lab.org/. https://doi.org/10.1371/journal.pone.0214510.g002 showed a high level of stress at the mandibular notch as well as from the condyle through the ramus in a descending gradient. Quantitative methods were utilized to better understand the intensity of the stresses in the mandibles. Boxplots of stress distributions (assuming Quasi-Ideal Mesh (QIM)) are depicted by taxon for both bite configurations in Fig 3. Average MWAM values and the quartiles values for the boxplots are listed in S1 Table (lateral cases) and S2 Table (orthal cases) of the supple- mentary information: M25 corresponds to the 25th percentile, M50 to the 50th percentile, M75 to the 75th percentile and M95 to the 95th percentile. An important caveat to the FEA method is that singular and unusually high stress values can appear where the boundary conditions are set. These numerical singularities are an arte- fact of the applied mathematical approach, inflated by the constraints imposed on the model [70]. In these areas, stresses have the tendency to increase in value towards infinity. However, these singularities are not related to any biological process and thus should be excluded from qualitative analyses of stress. The suggestions of Walmsley et al. [82], who recommend using the 95th percentile (M95) as the peak stress value in analyses in order to avoid these artificially high values are followed here. PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 8 / 21 Why ruminating ungulates chew sloppily Fig 3. Box-plots of modelled von mises stress distributions in the mandible by taxon for each bite configuration and bite position. Lateral cases 1) in M1 2) in M2 and 3) in M3. Orthal cases 4) in M1 5). In molar M2 and 6) in M3. Quasi-Ideal Mesh (QIM) was assumed. Upper and lower whiskers indicate 90% thresholds of the von Mises stress distributions. https://doi.org/10.1371/journal.pone.0214510.g003 Comparison of averaged stress values (MWAM) (Fig 4) and peak stress values (Fig 5), reveal differences between Artiodactyla and Perissodactyla in both bite configurations and suggest that Perissodactyla have stiffer mandibles than Artiodactyla. However, both averaged stress values (MWAM) and peak stress values (M95) also appear to capture differences in dietary traits (Figs 4 and 5). For instance, the mandibles of D. bicornis, D. sumatrensis and T. terrestris, all categorized as general browsers in this comparison, are characterized by lower stress values than general grazers and mixed feeders (intermediate stress values) and high-level browsers (highest stress values), suggesting that these taxa have stiffer mandibles despite their browsing activity). Cluster analyses show that species group according to their taxonomic group if averaged MWAM are used (Fig 6A) although the camel and the lama both cluster with perissodactyls. The same result indicating affinities between perissodactyls (Perissodactyla) and camelids (Tylopoda) is obtained if percentiles are used (Fig 6B). Cophenetic correlation coefficient for both clusters was obtained (RMWAM = 0.64 and Rpercentiles = 0.66). Fig 7 shows the results of the ancestral state reconstruction for the averaged MWAM stresses in the six bite scenarios. These values were mapped on the phylogeny using a maxi- mum-likelihood ancestral character estimation method based on a Brownian motion model of PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 9 / 21 Why ruminating ungulates chew sloppily Fig 4. Boxplots of averaged stress MWAM values of all species grouped by order (Artiodactyla and Perissodactyla) and by dietary categories (GGR: General grazer, GBR: General browser, OMF: Open-habitat mixed-feeder and HBR: High-level browser). Lateral cases 1) in M1 2) in M2 and 3) in M3. Orthal cases 4) in M1 5) in M2 and 6) in M3. Central line is the median, whiskers represent the range. https://doi.org/10.1371/journal.pone.0214510.g004 evolution. The figure depicts the same results of cluster analysis where species presented simi- lar results according to their taxonomic group but camelids with the order of magnitude of perissodactyls. Finally, Significant phylogenetic signal was found for MWAM data of the six bite cases (Kmult = 1.06; p-value = 7e-4; 9999 permutations) and percentile data of stress (Kmult = 1.24; p-value = 3e-4; 9999 permutations). Intervals’ method Biplots of stress space for the first two principal components of the PCA indicate that multiple variables contribute importantly to the PCs analysed (Fig 8). Variance of the first two compo- nents explain more than 80% of the variance in lateral and orthal cases (see S4 Table for the percent of the variance explained by the first ten components in all six cases). Perissodactyla are characterized by negative PC1 values which represents overall lower stress. Tapirus terres- tris exhibiting also negative scores in PC2, because it has large areas of very low stress whereas other Perissodactyla with fewer and smaller areas of low stress present positive PC2 values. Artiodactyla exhibit primarily positive scores in PC1 because their mandibles are characterized by larger areas of high and moderately high stress. A notable exception to this pattern are the Tylopoda (C. dromedarius and L. glama) which, for the lateral bite case, present negative PC1 PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 10 / 21 Why ruminating ungulates chew sloppily Fig 5. Boxplots of the peak stress values (assuming the 95th percentile as the peak stress value) for all species grouped by order (Artiodactyla and Perissodactyla) and by dietary categories (GGR: General grazer, GBR: General browser, OMF: Open-habitat mixed-feeder and HBR: High-level browser). Lateral cases 1) in M1 2) in M2 and 3) in M3. Orthal cases 4) in M1 5) in M2 and 6) in M3. Central line is the median, whiskers represent the range. https://doi.org/10.1371/journal.pone.0214510.g005 Fig 6. Hierarchical clustering analysis using Ward’s method for a) averaged stress values (MWAM) values and b) percentiles in all the six biting cases. https://doi.org/10.1371/journal.pone.0214510.g006 PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 11 / 21 Why ruminating ungulates chew sloppily Fig 7. Stress values MWAM on the phylogeny for the six biting cases: Lateral cases 1) in M1 2) in M2 and 3) in M3. Orthal cases 4) in M1 5) in M2 and 6) in M3. The colour ranges from red representing higher average stress values, to blue, representing lower stress values. https://doi.org/10.1371/journal.pone.0214510.g007 scores (more areas of intermediate-low stress) close to some of the Perissodactyla. For the orthal bite case, L. glama plots near the Ruminantia (A. buselaphus, A. melampus, G. Camelo- pardalis and L. walleri) and C. dromedaries, with large areas of intermediate values of stress, plots between the Ruminantia and the Perissodactyla. Finally, Significant phylogenetic signal was found for MWAM data of the six bite cases (Kmult = 0.78; p-value = 0.0012; 9999 permu- tations). The PC plots do not clearly differentiate dietary traits. Discussion: Diet-related stress patterns versus phylogenetic patterns In this study, we generated FEA models to test whether the biomechanical performance of the ungulate mandible reflects more feeding adaptations or historical contingencies. Our data include eleven species of ungulate with distinct, well-established dietary traits but belonging to two different clades, Perissodactyla and Artiodactyla. Our results reveal clear differences between Perissodactyla and Artiodactyla (Figs 2, 4 and 5). The mandibles of these taxa differ in biomechanical performance during orthal and lateral PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 12 / 21 Why ruminating ungulates chew sloppily Fig 8. Biplot of PCAs based on the correlation matrix for the six biting cases: Lateral cases 1) in M1 2) in M2 and 3) in M3. Orthal cases 4) in M1 5) in M2 and 6) in M3. The loadings for each variable are coloured according stress intensity. Red colours indicate high stress, blue indicate low stress. X-axis: PC1. Y- axis: PC2. https://doi.org/10.1371/journal.pone.0214510.g008 biting and in their mandibular movements during chewing and suggesting that Perissodactyla have stiffer mandibles than Artiodactyla. However, differences are also evident between taxa with distinct dietary traits. For instance, von Mises stress distributions (Figs 4 and 5), indicate that general browsers have comparatively stiffer mandibles than general grazers and mixed feeders while highly specialized browsers possess the most fragile mandibles. PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 13 / 21 Why ruminating ungulates chew sloppily These patterns are not necessarily reflective of dietary adaptations only, but may be biased by phylogeny, particularly when closely related taxa have undergone an adaptive radiation into different dietary niches. The latter scenario is suggested by the significant differences observed in this study between members of Artiodactyla (Ruminantia and Tylopoda) and Perissodactyla in both average and peak von Mises stress values which seem independent of dietary trait assignment (S1 and S2 Tables of the Supplementary information). Cluster analysis of stress val- ues (Fig 6A and 6B) and the intervals’ method (Fig 8) shed further light on this issue. Cluster analysis of averaged stress values (MWAM) and the percentiles both resulted in two distinct clusters: one including Tylopoda and Perissodactyla and the other including the Ruminantia (Fig 6). Similarly, intervals method reveal that Perissodactyla and Ruminantia are different in mandibular stress distribution because they are placed in a different area of the PC plots. Apparent phylogenetic differences are more evident between Perissodactyla and ruminant species. The strong phylogenetic signal supported by Kmult values in mandibular stress levels we observe between ruminants and Perissodactyla is likely related to major divergence in the digestive physiology of the two groups. Perissodactyla are hindgut fermenters and have shorter gut passage times compared to Ruminantia and Tylopoda of similar body mass and feeding trait. Hindgut fermenters thus need to compensate for this with higher food intake rates [83– 85]. Fletcher et al. [36] suggest that extant hindgut fermenters have more robust mandibles than ruminants—regardless of their dietary traits—as a result of their higher food intake rates. Conversely, the more gracile mandibles of ruminants likely result from their ability to re-mas- ticate their food several times after initial ingestion by regurgitating and chewing their cud. Rumination results in re-hydration and softening of food as well as internal ‘washing’ of ingesta in the rumen [86], which removes hard, exogenous grid from food surfaces. Moreover, splitting mastication into several cycles may reduce the overall mechanical demand of chewing by permitting less thorough mastication and necessitating lower bite forces upon initial inges- tion. All of these variables may contribute to a release on masticatory demands. Although Fletcher et al. [36] investigated orthal bite configurations only; our FEA results are nonetheless consistent with their data in both orthal and lateral bite configurations. In general, it would appear that ruminants can afford to chew sloppily upon initial digestion, while hindgut fer- menters cannot, necessitating stronger mandibles with higher stress resistance in the latter group. The separation of Ruminantia from Tylopoda, both of which belong to the order Cetartio- dactyla, suggest a non-phylogenetic influence on the evolution of even-toed ungulate mandib- ular morphology. For example, C. dromedarius and L. glama clearly differ from ruminants in their stress distribution. In fact, both cluster statistics and intervals method PC plots indicate affinities between perissodactyls (Perissodactyla) and camelids (Tylopoda). Interestingly, clus- ter analyses using averaged stress values (MWAM) suggest that inside the rhinocerontids these affinities in mandibular morphologies could reflect dietary preferences (Fig 6). In the PC plots of lateral biting configurations, both Tylopoda taxa have intermediate stress values, which place them closer to the Perissodactyla, but also high average stress values, align- ing them with the ruminants. For orthal biting configurations, L. glama plots close to the rumi- nants in all analyses. These results are consistent with the cluster analysis, which plots the Tylopoda together with the Perissodactyla. Considering the dominance of the phylogenetic signal as well as the specific digestive physiology, it is no surprise that both tylopods show spe- cific traits, different from the Ruminantia. The lower stress values observed in the Tylopoda mandibles suggest greater morphological robusticity and may indicate greater reliance on comminution immediately following ingestion than in other ruminating taxa. Ruminatia meanwhile, appear to shift comminution to later rumination cycles, when cuds are soft and water soaked, thus requiring less chewing force. In fact, differences in the digestive-washing PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 14 / 21 Why ruminating ungulates chew sloppily system between Ruminatia and Tylopoda could also explain the different biomechanical response of the mandible [86]. Further, our FEA data suggest that within the order Perissodactyla taxa also cluster accord- ing to dietary trait. This pattern is most notable amongst the Rhinocerotidae; D. bicornis and D. sumatrensis (general browsers) share a subcluster while C. simum and E. quagga (general grazers) also show close affinities. These results suggest that diet has influenced the bio- mechanical signal in Rhinocerotidae to a greater degree than any of the other ungulate families investigated in this study. Interestingly, tapirs (T. terrestris) differ from the other perissodactyla. Cluster analysis sug- gests that T. terrestris has affinities with this group, possibly due to a strong phylogenetic signal [30], however it is almost consistently isolated in the lower left space of all intervals method biplots (but see Fig 8, Case 3), suggesting significant differences with the rest of the Perissodac- tyla. In particular, T. terrestris has the lowest average and peak stress values (see Fig 3 and S1 and S2 Tables of the Supplementary information). Moreover, the distribution of stress (Fig 2) shows more areas of low stress values, especially in the lateral bite configuration. Fig 8 sustains these results showing low values of average stress in tapir in both lateral and orthal bite. But, compared the other species in lateral occlusion, values are not too much lower (specifically compared with D. sumatrensis, D. bicornis or C. dromedarius) while in orthalbite values are the lowest by far. These results are due to a lack of she-cutting functionality in tapirs that need to be compensated with high orthal loads that require large masseters and masseter insertions areas. The high ortal loads result in significant attritional wear in tapirs that does not occur in the perissodactyls with a lateral power stroke like equids and rhinocerontids. These data con- firms that the mandible of T. terrestris can be not particularly adapted to perform the lateral mode of power stroke and thus a more robust mandible should be necessary for comminution its cuts performing the orthal occlusion during chewing [87].Finally, average and peak stress values (see Fig 3 and S1 and S2 Tables of the Supplementary information) indicate that Litocra- nius walleri has the most fragile mandible. This finding is not surprising because L. walleri, a highly selective folivorous browser [88] forages on soft and easy to comminute. In browsing ruminants, compression of the bolus extraction of cell contents is accomplished by lateral occlusal movements which take advantage of the compression basins formed between enamel ridge 2 and 3 (e.g. [89]). The mandible (and brachydont cheek teeth) of L. walleri are hence assumed to not be adapted to withstand high bite forces or loads, and thus show the highest stress patterns when compared to mixed-feeding or grazing taxa among the ruminating group. In summary, we conclude that the biomechanics of the ungulate mandible are influenced by both diet and phylogeny although more taxa might be needed to expand this conclusion. In general, a stronger phylogenetical signal seems to explain why Perissodactyla have stiffer man- dibles than Artiodactyla. This difference is most pronounced between Perissodactyla and ruminant species. Perissodactyla rely heavily on thoroughly chewing their food upon inges- tion, demanding the generation of high bite forces, while ruminants shift comminution to later rumination cycles when cuds are soft and water soaked, requiring lower bite forces to obtain physiologically sufficient mechanical disintegration. We therefore suggest that rumi- nants can afford to chew sloppily regardless of ingesta, while hindgut fermenters cannot. Additionally, our data suggest a second adaptive layer resulting from diet suggesting that, within the orders Artiodactyla and Perissodactyla, mandibular morphologies also reflect the masticatory demands of specific ingesta (e.g. Typlopoda, Rhinocerotidae). A similar observa- tion was made by [27], who found that mandible movement reflects a mammal’s phylogenetic history as well as its current feeding behaviour. Our findings thus agree with classical interpre- tations of ungulate mandibular morphology [6,17,27–29]. PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 15 / 21 Why ruminating ungulates chew sloppily Supporting information S1 Table. FEA results for the lateral biting cases. Number of mandible mesh elements and statistics: Arithmetic Mean (AM), Mesh-Weighted Arithmetic Mean (MWAM), Percentage Error of the Arithmetic Mean (PEofAM), Median (M), Mesh-Weighted Median (MWM), Per- centage Error of the Median (PEofM) and the value quartiles (M25, M50, M75 and M95) according to Marcé-Nogué et al. 2016 (DOCX) S2 Table. FEA results for the orthal biting cases. Number of mandible mesh elements and statistics: Arithmetic Mean (AM), Mesh-Weighted Arithmetic Mean (MWAM), Percentage Error of the Arithmetic Mean (PEofAM), Median (M), Mesh-Weighted Median (MWM), Per- centage Error of the Median (PEofM) and the value quartiles (M25, M50, M75 and M95) according to Marcé-Nogué et al. 2016 (DOCX) S3 Table. Convergence of the R2 values of the PC scores. Each value is the R2 for a different pair of PCAs of the correlation matrix. Each PC was correlated with the equivalent PC of the PCA developed using a larger number of intervals. (DOCX) S4 Table. Percent of the variance of the PC scores. For the six biting cases (DOCX) S1 Document. Intervals’ method data of Von Mises stress when N = 100 M for the six bit- ing cases. (XLSX) S1 Fig. Plots displaying the first two PCs of the different PCAs for N = 10, 25, 50, 75 and 100 for case 1. Lateral biting in the first molar. The species are coloured by order: blue: Peri- ssodactyla and brown Cetartiodactyla. The axes of each pair of PCs are in the same scale. (TIF) S2 Fig. Plots displaying the first two PCs of the different PCAs for N = 10, 25, 50, 75 and 100 for case 2. Lateral biting in the second molar. The species are coloured by order: blue: Perissodactyla and brown Cetartiodactyla. The axes of each pair of PCs are in the same scale. (TIF) S3 Fig. Plots displaying the first two PCs of the different PCAs for N = 10, 25, 50, 75 and 100 for case 3. Lateral biting in the third molar. The species are coloured by order: blue: Peri- ssodactyla and brown Cetartiodactyla. The axes of each pair of PCs are in the same scale. (TIF) S4 Fig. Plots displaying the first two PCs of the different PCAs for N = 10, 25, 50, 75 and 100 for case 4. Orthal biting in the first molar. The species are coloured by order: blue: Peri- ssodactyla and brown Cetartiodactyla. The axes of each pair of PCs are in the same scale. (TIF) S5 Fig. Plots displaying the first two PCs of the different PCAs for N = 10, 25, 50, 75 and 100 for case 5. Orthal biting in the second molar. The species are coloured by order: blue: Peri- ssodactyla and brown Cetartiodactyla. The axes of each pair of PCs are in the same scale. (TIF) S6 Fig. Plots displaying the first two PCs of the different PCAs for N = 10, 25, 50, 75 and 100 for case 6. Orthal biting in the third molar. The species are coloured by order: blue: PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 16 / 21 Why ruminating ungulates chew sloppily Perissodactyla and brown Cetartiodactyla. The axes of each pair of PCs are in the same scale. (TIF) Acknowledgments The authors want to thank Wencheng Ge and Jie Zhou (Guilin University of Electronic Tech- nology) for their help building the three-dimensional models, Dr. Jennifer Leichliter for the English review and comments in the manuscript, Dr. Thomas A. Püschel and Dr. Soledad de Esteban for his comments in the manuscript and the colleagues at Steinmann Institut Bonn as well as Dr. Jan-Hendrik Buhk at Universitätsklinikum Eppendorf for CT scanning of speci- mens. We want to thank YXLON International GmbH for the CT-scan of Dicerorhinus suma- trensis. The authors acknowledge the Museum für Naturkunde Berlin for specimen loan and A.C. Sharp for the original files of the scanned G. Camelopardalis and C. dromedarius. L. glama and T. terrestris are hosted in DigiMorph.org and were scanned by The University of Texas High-Resolution X-ray CT Facility (UTCT), with the NSF grants IIS-0208675 and IIS- 9874781. Last but not least, comments and review from Dr. Alana Sharp and Dr. Gabriele San- salone as well as the editor Dr. Carlo Meloro greatly improved the previous version of this manuscript. Author Contributions Conceptualization: Daniela E. Winkler, Jordi Marcé-Nogué. Methodology: Zupeng Zhou, Daniela E. Winkler, Josep Fortuny, Thomas M. Kaiser, Jordi Marcé-Nogué. Writing – original draft: Zupeng Zhou, Daniela E. Winkler, Josep Fortuny, Thomas M. Kai- ser, Jordi Marcé-Nogué. Writing – review & editing: Zupeng Zhou, Daniela E. Winkler, Josep Fortuny, Thomas M. Kaiser, Jordi Marcé-Nogué. References 1. Maynard-Smith J, Savage R. The mechanics of mammalian jaws. Sch Sci Rev. 1959; 289–301. 2. Hylander WL. Mandibular function in Galago crassicaudatus and Macaca fascicularis: An in vivo approach to Stress Analysis of the mandible. J Morphol. 1979; 159: 253–296. https://doi.org/10.1002/ jmor.1051590208 PMID: 105147 3. Hogue AS. Mandibular corpus form and its functional significance: Evidence from marsupials. In: Vin- yard C, Ravosa M, Wall C, editors. Primate Craniofacial Function and Biology. New York, NY: Springer Academic Publishers; 2008. pp. 329–356. https://doi.org/10.1007/978-0-387-76585-3 4. Marcé-Nogué J, Püschel TA, Kaiser TM. A biomechanical approach to understand the ecomorphologi- cal relationship between primate mandibles and diet. Sci Rep. 2017; 7: 8364. https://doi.org/10.1038/ s41598-017-08161-0 PMID: 28827696 5. Hogue AS, Ravosa MJ. Transverse masticatory movements, occlusal orientation, and symphyseal fusion in selenodont artiodactyls. J Morphol. 2001; 249: 221–241. https://doi.org/10.1002/jmor.1051 PMID: 11517466 6. Spencer LM. Morphological Correlates of Dietary Resource Partitioning in the African Bovidae. J Mam- mal. 1995; 76: 448–471. https://doi.org/10.2307/1382355 7. Varela L, Fariña RA. Masseter moment arm as a dietary proxy in herbivorous ungulates. J Zool. 2015; 296: 295–304. https://doi.org/10.1111/jzo.12246 8. Freeman PW. Specialized insectivory: beetle-eating and moth-eating molossid bats. J Mammal. 1979; 60: 467–479. https://doi.org/10.2307/1380088 9. Freeman PW. Correspondence of food habits and morphology in insectivorous bats. J Mammal. 1981; 62: 166–173. https://doi.org/10.2307/1380489 PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 17 / 21 Why ruminating ungulates chew sloppily 10. Freeman PW. Functional cranial analysis of large animalivorous bats (Microchiroptera). Biol J Linn Soc. 1984; 21: 387–408. https://doi.org/10.1111/j.1095-8312.1984.tb01601.x 11. Freeman PW. Frugivorous and animalivorous bats (Microchiroptera): dental and cranial adaptations. Biol J Linn Soc. 1988; 33: 249–272. https://doi.org/10.1111/j.1095-8312.1988.tb00811.x 12. Freeman PW. Macroevolution in Microchiropteria: Recoupling morphology and ecology with phylogeny. Evol Ecol Res. 2000; 2: 317–335. 13. Bouvier M. Biomechanical scaling of mandibular dimensions in New World Monkeys. Int J Primatol. 1986; 7: 551–567. https://doi.org/10.1007/BF02736661 14. Bouvier M. A biomechanical analysis of mandibular scaling in old world monkeys. Am J Phys Anthropol. 1986; 69: 473–482. https://doi.org/10.1002/ajpa.1330690406 15. Taylor AB. Feeding behavior, diet, and the functional consequences of jaw form in orangutans, with implications for the evolution of Pongo. J Hum Evol. 2006; 50: 377–393. https://doi.org/10.1016/j.jhevol. 2005.10.006 PMID: 16413045 16. Radinsky LB. Evolution of skull shape in carnivores 1. representative modern carnivroes. Biol J Linn Soc. 1981; 15: 369–388. https://doi.org/10.1111/j.1095-8312.1981.tb01657.x 17. Radinsky LB. Evolution of skull shape in carnivores: 2. Additional modern carnivores. Biol J Linn Soc. 1981; 16: 337–355. https://doi.org/10.1111/j.1095-8312.1981.tb01657.x 18. Biknevicius AR, Ruff CB. The structure of the mandibular corpus and its relationship to feeding behav- iours in extant carnivorans. J Zool. 1992; 228: 479–507. https://doi.org/10.1111/j.1469-7998.1992. tb04450.x 19. Christiansen P, Wroe S. Bite forces and evolutionary adaptations to feeding ecology in carnivores. Ecol- ogy. 2007; 88: 347–358. https://doi.org/10.1890/0012-9658(2007)88[347:BFAEAT]2.0.CO;2 PMID: 17479753 20. Tseng ZJ, Su DF, Wang X, White SC, Ji X. Feeding capability in the extinct giant Siamogale melilutra and comparative mandibular biomechanics of living Lutrinae. Sci Rep. Springer US; 2017; 7: 15225. https://doi.org/10.1038/s41598-017-15391-9 PMID: 29123190 21. Tseng ZJ, Grohé C, Flynn JJ. A unique feeding strategy of the extinct marine mammal Kolponomos: convergence on sabretooths and sea otters. Proc R Soc B Biol Sci. 2016; 283: 20160044. https://doi. org/10.1098/rspb.2016.0044 PMID: 26936242 22. Tseng ZJ, Flynn JJ. Are cranial biomechanical simulation data linked to known diets in extant taxa? A method for applying diet-biomechanics linkage models to infer feeding capability of extinct species. PLoS One. 2015; 10: 1–25. https://doi.org/10.1371/journal.pone.0124020 PMID: 25923776 23. Mihlbachler MC, Rivals F, Solounias N, Semprebon GM. Dietary Change and Evolution of Horses in North America. Science (80-). 2011; 331: 1178–1181. https://doi.org/10.1126/science.1196166 PMID: 21385712 24. Clauss M, Kaiser TM, Hummel J. The Morphophysiological Adaptations of Browsing and Grazing Mam- mals. In: Gordon Iain J., Prins HHT, editor. The ecology of browsing and grazing. Berlin: Springer; 2008. pp. 47–88. https://doi.org/10.1007/978-3-540-72422-3_3 25. Janis CM. Correlations between craniodental morphology and feeding behavior in ungulates: reciprovo- cal illumanation between living and fossil taxa. Functional morphology in vertebrate paleontology. 1995. pp. 76–98. 26. Mendoza M, Palmqvist P. Hypsodonty in ungulates: an adaptation for grass consumption or for foraging in open habitat? J Zool. 2008; 274: 134–142. https://doi.org/10.1111/j.1469-7998.2007.00365.x 27. Goldberg LJ, Gerstner GE. Species-Specific Morphology of Masticatory Jaw Movements. Behaviour. 1994; 128: 229–253. https://doi.org/10.1163/156853994X00271 28. Pérez-Barberı́a FJ, Gordon IJ. The functional relationship between feeding type and jaw and cranial morphology in ungulates. Oecologia. 1999; 118: 157–165. https://doi.org/10.1007/s004420050714 PMID: 28307690 29. Pérez-Barberı́a FJ, Gordon IJ. Relationships between oral morphology and feeding style in the Ungu- lata: a phylogenetically controlled evaluation. Proc Biol Sci. 2001; 268: 1023–32. https://doi.org/10. 1098/rspb.2001.1619 PMID: 11375086 30. Raia P, Carotenuto F, Meloro C, Piras P, Pushkina D. The shape of contention: adaptation, history, and contingency in ungulate mandibles. Evolution (N Y). 2010; 64: 1489–1503. https://doi.org/10.1111/j. 1558-5646.2009.00921.x PMID: 20015238 31. Gill PG, Purnell MA, Crumpton N, Brown KR, Gostling NJ, Stampanoni M, et al. Dietary specializations and diversity in feeding ecology of the earliest stem mammals. Nature. Nature Publishing Group; 2014; 512: 303–305. https://doi.org/10.1038/nature13622 PMID: 25143112 PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 18 / 21 Why ruminating ungulates chew sloppily 32. Gröning F, Liu J, Fagan MJ, O’Higgins PM. Why do humans have chins? Testing the mechanical signifi- cance of modern human symphyseal morphology with finite element analysis. Am J Phys Anthropol. 2011; 144: 593–606. https://doi.org/10.1002/ajpa.21447 PMID: 21404235 33. Tseng ZJ, Binder WJ. Mandibular biomechanics of Crocuta crocuta, Canis lupus, and the late Miocene Dinocrocuta gigantea (Carnivora, Mammalia). Zool J Linn Soc. 2010; 158: 683–696. https://doi.org/10. 1111/j.1096-3642.2009.00555.x 34. Püschel TA, Marcé-Nogué J, Kaiser TM, Brocklehurst RJ, Sellers WI. Analyzing the sclerocarpy adap- tations of the Pitheciidae mandible. Am J Primatol. 2018; 80: e22759. https://doi.org/10.1002/ajp.22759 PMID: 29664191 35. Piras P, Maiorino L, Teresi L, Meloro C, Lucci F, Kotsakis T, et al. Bite of the cats: Relationships between functional integration and mechanical performance as revealed by mandible geometry. Syst Biol. 2013; 62: 878–900. https://doi.org/10.1093/sysbio/syt053 PMID: 23925509 36. Fletcher TM, Janis CM, Rayfield EJ. Finite Element Analysis of Ungulate Jaws: Can mode of digestive physiology be determined? Palaeontol Electron. 2010; 13. Available: http://www.uv.es/pardomv/pe/ 2010_3/234/234.pdf 37. Neenan JM, Ruta M, Clack J a., Rayfield EJ. Feeding biomechanics in Acanthostega and across the fish-tetrapod transition. Proc R Soc B Biol Sci. 2014; 281: 20132689–20132689. https://doi.org/10. 1098/rspb.2013.2689 PMID: 24573844 38. Serrano-Fochs S, De Esteban-Trivigno S, Marcé-Nogué J, Fortuny J, Fariña RA. Finite Element Analy- sis of the Cingulata Jaw: An Ecomorphological Approach to Armadillo’s Diets. PLoS One. 2015; 10: e0120653. https://doi.org/10.1371/journal.pone.0120653 PMID: 25919313 39. O’Higgins PM, Cobb SN, Fitton LC, Gröning F, Phillips R, Liu J, et al. Combining geometric morpho- metrics and functional simulation: An emerging toolkit for virtual functional analyses. J Anat. 2011; 218: 3–15. https://doi.org/10.1111/j.1469-7580.2010.01301.x PMID: 20880075 40. Fortuny J, Marcé-Nogué J, DE Esteban-Trivigno S, Gil L, Galobart Á. Temnospondyli bite club: ecomor- phological patterns of the most diverse group of early tetrapods. J Evol Biol. 2011; 24: 2040–54. https:// doi.org/10.1111/j.1420-9101.2011.02338.x PMID: 21707813 41. Attard MRG, Parr WCH, Wilson LAB, Archer M, Hand SJ, Rogers TL, et al. Virtual reconstruction and prey size preference in the mid cenozoic thylacinid, Nimbacinus dicksoni (Thylacinidae, Marsupialia). PLoS One. 2014; 9. https://doi.org/10.1371/journal.pone.0093088 PMID: 24718109 42. Fortuny J, Marcé-Nogué J, Steyer J-S, de Esteban-Trivigno S, Mujal E, Gil L. Comparative 3D analyses and palaeoecology of giant early amphibians (Temnospondyli: Stereospondyli). Sci Rep. 2016; 6: 30387. https://doi.org/10.1038/srep30387 PMID: 27457883 43. Piras P, Sansalone G, Teresi L, Moscato M, Profico A, Eng R, et al. Digging adaptation in insectivorous subterranean eutherians. The enigma of Mesoscalops montanensis unveiled by geometric morpho- metrics and finite element analysis. J Morphol. 2015; 276: 1157–1171. https://doi.org/10.1002/jmor. 20405 PMID: 26134243 44. Sharp AC. Comparative finite element analysis of the cranial performance of four herbivorous marsupi- als. J Morphol. 2015; 276: 1230–1243. https://doi.org/10.1002/jmor.20414 PMID: 26193997 45. Jones MEH, Gröning F, Dutel H, Sharp AC, Fagan MJ, Evans SE. The biomechanical role of the chon- drocranium and sutures in a lizard cranium. J R Soc Interface. 2017; 14: 20170637. https://doi.org/10. 1098/rsif.2017.0637 PMID: 29263126 46. Panagiotopoulou O, Iriarte-Diaz J, Wilshin S, Dechow PC, Taylor AB, Mehari Abraha H, et al. In vivo bone strain and finite element modeling of a rhesus macaque mandible during mastication. Zoology. Elsevier; 2017; 0–1. https://doi.org/10.1016/j.zool.2017.08.010 PMID: 29037463 47. Farke AA. Frontal sinuses and head-butting in goats: a finite element analysis. J Exp Biol. 2008; 211: 3085–3094. https://doi.org/10.1242/jeb.019042 PMID: 18805807 48. Tseng ZJ. Cranial function in a late Miocene Dinocrocuta gigantea (Mammalia: Carnivora) revealed by comparative finite element analysis. Biol J Linn Soc. Wiley Online Library; 2009; 96: 51–67. https://doi. org/10.1111/j.1095-8312.2008.01095.x 49. Aquilina P, Chamoli U, Parr WCH, Clausen PD, Wroe S. Finite element analysis of three patterns of internal fixation of fractures of the mandibular condyle. Br J Oral Maxillofac Surg. British Association of Oral and Maxillofacial Surgeons; 2013; 51: 326–331. https://doi.org/10.1016/j.bjoms.2012.08.007 PMID: 22981343 50. Figueirido B, Tseng ZJ, Serrano-Alarcon FJ, Martin-Serra A, Pastor JFF, Serrano-Alarcón FJ, et al. Three-dimensional computer simulations of feeding behaviour in red and giant pandas relate skull bio- mechanics with dietary niche partitioning. Biol Lett. 2014; 10: 20140196–20140196. https://doi.org/10. 1098/rsbl.2014.0196 PMID: 24718096 PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 19 / 21 Why ruminating ungulates chew sloppily 51. Fish JF, Stayton CT. Morphological and mechanical changes in juvenile red-eared slider turtle (Tra- chemys scripta elegans) shells during ontogeny. J Morphol. 2014; 275: 391–397. https://doi.org/10. 1002/jmor.20222 PMID: 24301373 52. Parr WCH, Wroe S, Chamoli U, Richards HS, McCurry MR, Clausen PD, et al. Toward integration of geometric morphometrics and computational biomechanics: New methods for 3D virtual reconstruction and quantitative analysis of Finite Element Models. J Theor Biol. Elsevier; 2012; 301: 1–14. https://doi. org/10.1016/j.jtbi.2012.01.030 PMID: 22342680 53. Lautenschlager S, Gill PG, Luo Z-X, Fagan MJ, Rayfield EJ. The role of miniaturization in the evolution of the mammalian jaw and middle ear. Nature. Springer US; 2018; https://doi.org/10.1038/s41586-018- 0521-4 PMID: 30224748 54. Brassey CA, Gardiner JD, Kitchener AC. Testing hypotheses for the function of the carnivoran baculum using finite-element analysis. Proc R Soc B Biol Sci. 2018; 285: 20181473. https://doi.org/10.1098/rspb. 2018.1473 PMID: 30232157 55. Maiorino L, Farke AA, Kotsakis T, Teresi L, Piras P. Variation in the shape and mechanical performance of the lower jaws in ceratopsid dinosaurs (Ornithischia, Ceratopsia). J Anat. 2015; 227: 631–646. https://doi.org/10.1111/joa.12374 PMID: 26467240 56. Gröning F, Fagan MJJMJMJJMJMJJ, O’Higgins PM, O’Higgins P, O’Higgins PM , O’Higgins P, et al. The effects of the periodontal ligament on mandibular stiffness: A study combining finite element analy- sis and geometric morphometrics. J Biomech. Elsevier; 2011; 44: 1304–1312. https://doi.org/10.1016/j. jbiomech.2011.01.008 PMID: 21292267 57. Marcé-Nogué J, De Esteban-Trivigno S, Escrig C, Gil L. Accounting for differences in element size and homogeneity when comparing Finite Element models: Armadillos as a case study. Palaeontol Electron. 2016; 19: 1–22. 58. Marcé-Nogué J, De Esteban-Trivigno S, Püschel TA, Fortuny J. The intervals method: a new approach to analyse finite element outputs using multivariate statistics. PeerJ. 2017; 5: e3793. https://doi.org/10. 7717/peerj.3793 PMID: 29043107 59. Püschel TA, Marcé-Nogué J, Gladman JT, Bobe R, Sellers WI. Inferring locomotor behaviours in Mio- cene New World monkeys using finite element analysis, geometric morphometrics and machine-learn- ing classification techniques applied to talar morphology. J R Soc Interface. 2018; 15: 20180520. https://doi.org/10.1098/rsif.2018.0520 PMID: 30257926 60. Mendoza M, Janis CM, Palmqvist P. Characterizing complex craniodental patterns related to feeding behaviour in ungulates: a multivariate approach. J Zool. 2002; 258: S0952836902001346. https://doi. org/10.1017/S0952836902001346 61. Fraser D, Rybczynski N. Complexity of ruminant masticatory evolution. J Morphol. 2014; 275: 1093– 1102. https://doi.org/10.1002/jmor.20284 PMID: 24753298 62. Janis CM. An Estimation of Tooth Volume and Hypsodonty Indices in Ungulate Mammals, and the Cor- relation of these Factors with Dietary Preference [Internet]. Teeth Revisited: Proceedings of the VIIth International Symposium on Dental Morphology. Paris 20–24 Mai 1986. 1988. pp. 367–387. Available: http://www.worldcat.org/title/teeth-revisited-proceedings-of-the-viith-international-symposium-on- dental-morphology-paris-20-24-mai-1986/oclc/636547678 63. Gailer JP, Kaiser TM. Common solutions to resolve different dietary challenges in the ruminant denti- tion: The functionality of bovid postcanine teeth as a masticatory unit. J Morphol. 2014; 275: 328–341. https://doi.org/10.1002/jmor.20217 PMID: 24136643 64. Marcé-Nogué J, Fortuny J, Gil L, Galobart A. Using reverse engineering to reconstruct tetrapod skulls and analyse its feeding behaviour. Proceedings of the 13th International Conference on Civil, Structural and Environmental Engineering Computing. 2011. pp. 1–12. doi:10.4203/ccp.96.237 65. Doblaré M, Garcı́a JM, Gómez MJ. Modelling bone tissue fracture and healing: a review. Eng Fract Mech. 2004; 71: 1809–1840. https://doi.org/10.1016/j.engfracmech.2003.08.003 66. Gil L, Marcé-Nogué J, Sanchez M, Sánchez M. Insights into the controversy over materials data for the comparison of biomechanical performance in vertebrates. Palaeontol Electron. 2015;18.1.12A: 1–24. Available: http://palaeo-electronica.org/content/2015/1095-controversy-in-materials-data 67. Strait DS, Wang Q, Dechow PC, Ross CF, Richmond BG, Spencer MA, et al. Modeling elastic proper- ties in finite-element analysis: How much precision is needed to produce an accurate model? Anat Rec —Part A Discov Mol Cell Evol Biol. 2005; 283: 275–287. https://doi.org/10.1002/ar.a.20172 PMID: 15747346 68. Reilly DT, Burstein AH. The elastic and ultimate properties of compact bone tissue. J Biomech. 1975; 8: 393–405. https://doi.org/10.1016/0021-9290(75)90075-5 PMID: 1206042 69. Lin CP, Douglas WH. Structure-property relations and crack resistance at the bovine dentin-enamel junction. J Dent Res. 1994; 73: 1072–8. https://doi.org/10.1177/00220345940730050901 PMID: 8006234 PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 20 / 21 Why ruminating ungulates chew sloppily 70. Marcé-Nogué J, Fortuny J, Gil L, Sánchez M. Improving mesh generation in Finite Element Analysis for functional morphology approaches. Spanish J Palaeontol. 2015; 31: 117–132. 71. Fortuny J, Marcé-Nogué J, Konietzko-Meier D. Feeding biomechanics of Late Triassic metoposaurids (Amphibia: Temnospondyli): a 3D finite element analysis approach. J Anat. 2017; 230: 752–765. https://doi.org/10.1111/joa.12605 PMID: 28369819 72. Preuschoft H, Witzel U. Biomechanical investigations on the skulls of reptiles and mammals. Sencken- bergiana Lethaea. 2002; 82: 207–222. https://doi.org/10.1007/BF03043785 73. Alexander RM. Exploring biomechanics. Animals in motion. New York Scientific American Library Dis- tributed by W.H. Freeman; 1992. 74. Marcé-Nogué J, DeMiguel D, Fortuny J, de Esteban-Trivigno S, Gil L. Quasi-homothetic transformation for comparing the mechanical performance of planar models in biological research. Palaeontol Electron. 2013; 16. Available: http://palaeo-electronica.org/content/2013/468-quasihomothetic-transformation 75. Fortuny J, Marcé-Nogué J, Heiss E, Sanchez M, Gil L, Galobart À. 3D Bite Modeling and Feeding Mechanics of the Largest Living Amphibian, the Chinese Giant Salamander Andrias davidianus (Amphibia:Urodela). PLoS One. 2015; 10: e0121885. https://doi.org/10.1371/journal.pone.0121885 PMID: 25853557 76. Meers MB. Maximum Bite Force and Prey Size of Tyrannosaurus rex and Their Relationships to the Inference of Feeding Behavior. Hist Biol A J Paleobiol. 2003; 16: 1–12. https://doi.org/10.1080/ 0891296021000050755 77. Wroe S, McHenry CR, Thomason JJ. Bite club: comparative bite force in big biting mammals and the prediction of predatory behaviour in fossil taxa. Proc R Soc B Biol Sci. The Royal Society; 2005; 272: 619–625. https://doi.org/10.1098/rspb.2004.2986 PMID: 15817436 78. Adams DC. A generalized K statistic for estimating phylogenetic signal from shape and other high- dimensional multivariate data. Syst Biol. 2014; 63: 685–697. https://doi.org/10.1093/sysbio/syu030 PMID: 24789073 79. Toomey D. R for Data Science: Learn and explore the fundamentals of data science with R. 2014. 80. Corliss JO, Sneath PHA, Sokal RR. Numerical Taxonomy: The Principles and Practice of Numerical Classification. Trans Am Microsc Soc. 2006; https://doi.org/10.2307/3225339 81. R Development Core Team. R: A Language and Environment for Statistical Computing. R Found Stat Comput Vienna Austria. 2018; https://doi.org/10.1038/sj.hdy.6800737 PMID: 16106260 82. Walmsley CW, Smits PD, Quayle MR, McCurry MR, Richards HS, Oldfield CC, et al. Why the Long Face? The Mechanics of Mandibular Symphysis Proportions in Crocodiles. Aegerter CM, editor. PLoS One. 2013; 8: e53873. https://doi.org/10.1371/journal.pone.0053873 PMID: 23342027 83. Clauss M, Frey R, Kiefer B, Lechner-Doll M, Loehlein W, Polster C, et al. The maximum attainable body size of herbivorous mammals: Morphophysiological constraints on foregut, and adaptations of hindgut fermenters. Oecologia. 2003; 136: 14–27. https://doi.org/10.1007/s00442-003-1254-z PMID: 12712314 84. Clauss M, Fritz J, Bayer D, Nygren K, Hammer S, Hatt JM, et al. Physical characteristics of rumen con- tents in four large ruminants of different feeding type, the addax (Addax nasomaculatus), bison (Bison bison), red deer (Cervus elaphus) and moose (Alces alces). Comp Biochem Physiol—A Mol Integr Phy- siol. 2009; 152: 398–406. https://doi.org/10.1016/j.cbpa.2008.11.009 PMID: 19049897 85. Gomes VSDM, Correia MCR, De Lima HA, Alves MAS. Potential role of frugivorous birds (Passeri- formes) on seed dispersal of six plant species in a restinga habitat, southeastern Brazil. Rev Biol Trop. 2008; 56: 205–216. https://doi.org/10.1007/978-3-540-72422-3 PMID: 18624238 86. Clauss M, Hummel J. Physiological adaptations of ruminants and their potential relevance for produc- tion systems. Rev Bras Zootec. 2017; 46: 606–613. https://doi.org/10.1590/S1806- 92902017000700008 87. Harris JM. Evolution of feeding mechanisms in the family Deinotheriidae (Mammalia: Proboscidea). Zool J Linn Soc. 1975; 56: 331–362. https://doi.org/10.1111/j.1096-3642.1975.tb00275.x 88. Leuthold W. On the ecology of the gerenuk Litocranius walleri. J Anim Ecol. 1978; 47: 561–580. Avail- able: http://www.jstor.org/stable/3801 89. Kaiser TM, Fickel J, Streich WJ, Hummel J, Clauss M. Enamel ridge alignment in upper molars of rumi- nants in relation to their natural diet. J Zool. 2010; 281: 12–25. https://doi.org/10.1111/j.1469-7998. 2009.00674.x PLOS ONE | https://doi.org/10.1371/journal.pone.0214510 April 17, 2019 21 / 21