Variation in the shape and mechanical performance of the lower jaws in ceratopsid dinosaurs (Ornithischia, Ceratopsia)

Ceratopsidae represents a group of quadrupedal herbivorous dinosaurs that inhabited western North America and eastern Asia during the Late Cretaceous. Although horns and frills of the cranium are highly variable across species, the lower jaw historically has been considered to be relatively conservative in morphology. Here, the lower jaws from 58 specimens representing 21 ceratopsoid taxa were sampled, using geometric morphometrics and 2D finite element analysis (FEA) to explore differences in morphology and mechanical performance across Ceratopsoidea (the clade including Ceratopsidae, Turanoceratops and Zuniceratops). Principal component analyses and non‐parametric permuted manovas highlight Triceratopsini as a morphologically distinct clade within the sample. A relatively robust and elongate dentary, a larger and more elongated coronoid process, and a small and dorso‐ventrally compressed angular characterize this clade, as well as the absolutely larger size. By contrast, non‐triceratopsin chasmosaurines, Centrosaurini and Pachyrhinosaurini have similar morphologies to each other. Zuniceratops and Avaceratops are distinct from other taxa. No differences in size between Pachyrhinosaurini and Centrosaurini are recovered using non‐parametric permuted anovas. Structural performance, as evaluated using a 2D FEA, is similar across all groups as measured by overall stress, with the exception of Triceratopsini. Shape, size and stress are phylogenetically constrained. A longer dentary as well as a long coronoid process result in a lower jaw that is reconstructed as relatively much more stressed in triceratopsins.


Introduction
Ceratopsidae (Ornithischia, Ceratopsia) represents a group of quadrupedal herbivorous dinosaurs that dominated the terrestrial ecosystems of western North America during the Late Cretaceous (Dodson et al. 2004;, along with a sporadic occurrence in Asia (Xu et al. 2010).
The overall ceratopsid lower jaw is elongate and highly specialized with a robust, toothless and sharp-edged beak, that appears to be suitable for grasping and plucking (Ostrom, 1966;Dodson et al. 2004;Makovicky, 2012). The dentary bears a long tooth battery, with 25-40 doublerooted tooth positions, and multiple replacement teeth at each position (Makovicky, 2012). A ridge across the lateral surface of the dentary terminates at the rostral edge of the coronoid process, and the dentary reaches its maximum thickness at the base of the coronoid process (Dodson et al. 2004). Lower jaw motion is mostly restricted to vertical or near-vertical shear during mastication (Ostrom, 1966;Dodson et al. 2004;Brusatte, 2012).
Ceratopsids were hypothesized to have eaten highly fibrous plants (palms and cycads) (Ostrom, 1966;Krauss, 2001), because the musculoskeletal configuration in these animals is reconstructed to produce high bite forces during the feeding process (Ostrom, 1966;Bell et al. 2009;Tanoue et al. 2009) and also because they possessed a peculiar dentition acting to produce a strict shearing stroke (Mallon & Anderson, 2014a). Ceratopsids are considered as 'concentrate' feeders due to their narrow rostrum shape (Mallon & Anderson, 2014b).
Recent new discoveries (Ryan, 2007;Kirkland & DeBlieux, 2010;Sampson et al. , 2013Sullivan & Lucas, 2010;Farke et al. 2011;Fiorillo & Tykoski, 2012;Ryan et al. 2012;Wick & Lehman, 2013;Longrich, 2014;Brown & Henderson, 2015) have allowed many researchers to investigate the biomechanics of the ceratopsid feeding apparatus (Bell et al. 2009;Tanoue et al. 2009;Maiorino et al. 2014;Mallon & Anderson, 2015). Ostrom (1964Ostrom ( , 1966 modeled the ceratopsid lower jaw as a second-class lever, having a tooth row shifted posterior relative to that in the primitive condition, a depressed glenoid, and a dorsally elongated and laterally displaced coronoid process (Fig. 1). These modifications were accompanied by a great expansion of the frill muscles, where the adductors represents the major jaw muscle group (Holliday, 2009). The combination of the edentulous beak, a jaw articulation well below the level of the dentary tooth row, and a tall coronoid process indicates that high bite forces could be generated along the tooth row (Ostrom, 1966;Dodson et al. 2004;Makovicky, 2012). Recent works highlighted that relative higher bite forces could be generated by ceratopsids at the tooth positions, close to the jaw joint (Tanoue et al. 2009;Mallon & Anderson, 2015). Bell et al. (2009) investigated the structural performance of a lower jaw of Centrosaurus apertus using finite element analysis (FEA). The major stresses were concentrated along the leading edge of the coronoid process and transmitted along the longitudinal ridge. By contrast, the post-dentary region experienced low stress. Thus, Bell et al. (2009) concluded that the morphological adaptations of the lower jaw of Centrosaurus accommodate stress acting on the coronoid process, reducing the risk of breakage of this structure during mastication. Henderson (2010) performed a biomechanical analysis of ceratopsid skulls using beam theory, and hypothesized that centrosaurines differ from chasmosaurines in cranial depth, with correlated differences in the resistance to bending and torsional stresses, possibly facilitating niche partitioning between the two subfamilies and distinct diets. Mallon & Anderson (2013) used linear measurements to investigate cranial and lower jaw parameters of potential biomechanical relevance for niche partitioning. They confirmed some differences in rostral shape and cranial depth, but not in measurements directly taken from (or indirectly related to) the lower jaw.
Traditionally, centrosaurines and chasmosaurines have been considered to possess a similar lower jaw framework with no major morphological or mechanical differences beyond the predentary. In centrosaurines, the predentary's triturating surface is strongly canted laterally, contrasting with a more horizontally oriented surface in chasmosaurines . Despite the recent descriptions of new ceratopsid species (Ryan, 2007;Sampson et al. , 2013Farke et al. 2011;Fiorillo & Tykoski, 2012;Brown & Henderson, 2015) and new specimens (Mallon et al. 2014), a few contributions have highlighted possible lower jaw shape differences between centrosaurines and chasmosaurines and within each clade, particularly for the dentary or post-dentary bones. The sole exception concerns the observations by Mallon et al. (2011Mallon et al. ( , 2014, who noted that some anatomical features in Anchiceratops and Arrhinoceratops lower jaws (a prominent lateral ridge of the dentary) are shared with those of centrosaurines, whereas they are lacking in other chasmosaurines.
In this work, geometric morphometrics (GM) and FEA were combined to study a large dataset of ceratopsoid lower jaws (n = 58 representative of 21 species), in a phylogenetic context, to explore any morphological differences characterizing the several clades under investigation and to infer if mechanical performance differs among the selected ceratopsid clades. Two specific questions are addressed here: (i) what is the range of variation of the lower jaws in lateral view within Ceratopsidae? And (ii) can the lower jaw morphology be used to distinguish ceratopsid taxa? The null hypothesis is that centrosaurines and chasmosaurines have similar lower jaw morphologies, as generally accepted.

Materials
Images of 58 lower jaws in lateral view of ceratopsoid dinosaurs from both original photos and published pictures were collected. Twenty-one ceratopsid species are represented in the sample (one centrosaurine, three centrosaurins, nine non-triceratopsin chasmosaurines, four triceratopsins, four pachyrhinosaurins) and one ceratopsoid, Zuniceratops christopheri. Tables S1 and S2 in Fig. 1 Diagram of the ceratopsid lower jaw lever system (modified after Ostrom, 1964). a, horizontal distance between the jaw joint and the coronoid process; F, applied force of adductor muscle; h, height of the coronoid process above the jaw joint; l, lower jaw length; R, resistance force applied at the mid-point of the tooth row.
'Supporting Information' report the specimen list and the number of individuals for each species as well as the list of institutions where the several images were collected. A Canon 400D camera was used to collect most pictures for the sample; where necessary some images were taken from literature. The protocols of Marcus et al. (2000) and Mullin & Taylor (2002) were followed to minimize parallax and measurement error on the photographs.

Grouping ceratopsid clades
Historically, Centrosaurinae and Chasmosaurinae represent the two most consistently recognized clades within Ceratopsidae, but recent phylogenetic work has recognized additional clades within Centrosaurinae and Chasmosaurinae. Along with the description of the chasmosaurine Titanoceratops ouranos, Longrich (2011) proposed to name Triceratopsini for all species closer to Triceratops horridus than to Anchiceratops ornatus or Arrhinoceratops brachyops. Even if Titanoceratops (Longrich, 2011) was considered as Pentaceratops in this work, as suggested by Wick & Lehman (2013), the other indications of Longrich (2011) in considering derived Maastrichtian chasmosaurines Triceratops, Torosaurus, Nedoceratops, Ojoceratops and Eotriceratops pertaining to the Triceratopsini group were followed.
In 2012, Fiorillo & Tykoski suggested the name Pachyrhinosaurini, originally proposed by Sternberg (1950) as Pachyrhinosauridae, to define all centrosaurine ceratopsids more closely related to Pachyrhinosaurus canadensis than to Centrosaurus apertus. Therefore, Achelousaurus, Einiosaurus and Pachyrhinosaurus as members of Pachyrhinosaurini were considered here. In addition, recent cladistic analyses have highlighted how this group is phylogenetically supported (Sampson et al. 2013). Moreover, Centrosaurus, Coronosaurus and Styracosaurus (along with Rubeosaurus and Spinops) are here considered as Centrosaurini and defined as all centrosaurine ceratopsids more closely related to Centrosaurus apertus than to Pachyrhinosaurus canadensis. Its phylogenetic position is supported by the cladistic analysis performed by Sampson et al. (2013). The remaining chasmosaurines are indicated in the text using the formal designation 'non-triceratopsin chasmosaurines'.

GM
Geometric morphometrics is a suitable method to quantify morphological changes and to analyze phenotypic differences among taxa (Adams et al. 2004;Zelditch et al. 2012). Twelve landmarks and 21 semi-landmarks were digitized on each lower jaw in lateral view ( Fig. 2; Table S3) using TPSDIG2 v2.16 software (Rohlf, 2010). Scale bars were used to scale each digitized specimen.
Semi-landmarks are useful to capture morphological information of outlines where no homologous points can be detected. Curves or contours are assumed to be homologous among specimens (Bookstein et al. 2002;Perez et al. 2006). Semi-landmarks were digitized at equal distance along outlines drawn on the specimens.
Generalized Procrustes analysis (GPA; Bookstein, 1991) implemented using the procSym() function in the 'Morpho' R package (Schlager, 2013) was used to analyze shape differences among specimens in the lower jaw sample. GPA scales, aligns and rotates each landmark configuration to the unit centroid size (CS = the square root of sum of squared differences between landmarks from their centroid; Bookstein, 1986). CS represents the individual size of specimens in all analyses.
Taphonomic distortion presents a major issue in a GM study. It can produce deformations or damage to skulls or lower jaws, and therefore affect the results. Lower jaws are relatively flat bones and are usually less affected by post mortem distortion when compared with more complicated structures such as the horns and frills of ceratopsids. Landmarks were digitized only on the lower jaw profile, and the predentary was removed from the original configuration because it is often lost, damaged or restored with plaster. The articular bone in the landmark configuration was not included because it is often damaged or not easily distinguished in visual inspection. In a few cases, where the postdentary bones are missing, the function fixLMtps() from the 'Morpho' R package was used (Schlager, 2013) to estimate missing landmarks based on the three closest complete specimens, to minimize errors during the digitization process. Specimens lacking postdentary bones, such as Anchiceratops ornatus TMP (Royal Tyrrell Museum of Paleontology) 1983.01.01, Coahuilaceratops magnacuerna CPC (Coleccion Paleontologica de Coahuila) 276, Einiosaurus procurvicornis MOR373-7-15-6-13 and Ojoceratops fowleri SMP (State Museum of Pennsylvania) VP 1875, were also included in the sample to investigate lower jaw shape variation. However, the morphological considerations between taxa have been related to the dentary shape only, and not to the surangular-angular complex, because each species in the dataset is represented by at least one complete dentary.
After GPA, a between-group principal components analysis (bgPCA) was performed on the averaged Procrustes shape variables to identify orthogonal axes of maximal variation in the dataset. The bgPCA provides a projection of the data onto the PCs of the group means, leading to an ordination of the shape variables between the group means. The new axes are orthogonal and can be computed even when the data are not of full rank, as often happens for Procrustes shape coordinates (Mitteroecker & Bookstein, 2011). This method offers a good performance when the number of cases is smaller than the number of variables (Boulesteix, 2005;Sansalone et al. in press).  Table S3 for landmark definitions.

Phylogeny
A synthetic phylogenetic tree was built (Fig. 3), including all valid ceratopsoid species using the software MESQUITE 2.75 (Maddison & Maddison, 2011), based on previous published cladistic analyses. Branch lengths were calibrated in millions of years (Ma) based on published ages or age estimates for the formations yielding the fossils (see Appendix S1 and Tables S4 and S5 for details about topology and branch length calibration).

FEA
To reconstruct the lower jaw geometry, average values of Procrustes coordinates were used to generate a smooth and closed contour per species by spline approximation using the software COM-SOL MULTIPHYSICS 5.0 (http://www.comsol.com). A 2D solid geometry was used to study the mechanical behavior of the lower jaws by means of FEA. FEA is a popular mathematical tool among morphologists and vertebrate bio-mechanists to evaluate the strain and stress state within a solid, given its material properties, and under the actions of appropriate loads and constraints that represent a particular functional or behavioral scenario (McHenry et al. 2007;Rayfield, 2007Rayfield, , 2011Farke, 2008;Snively & Cox, 2008;Slater & Van Valkenburgh, 2009;Farke et al. 2010;Young et al. 2010;Tseng et al. 2011;Piras et al. 2012Piras et al. , 2013Piras et al. , 2015.
Here, linear elasticity in a 2D plane strain approximation was used to model static loading problems and quantitatively assess the biomechanical performance of the 2D geometries of lower jaws for comparative purposes in a specific phylogenetic context, that is which lower jaw shape is more stressed, within Ceratopsoidea, when constrained and loaded in a static condition?
Although 2D approximation of a 3D object is undoubtedly a simplification, it is argued that it is an appropriate (and necessary) approximation for the purposes of this study looking at general trends across taxa. The current arguments follow those enacted for a similar study on the lower jaw across the fish/tetrapod transition (Neenan et al. 2014). In the current work, it is difficult if not impossible to generate a number of 3D models on practical grounds, outlined below.
Specimen size as well as the inaccessibility of many specimens, within exhibit mounts, prohibited photography from additional views (ventral and postero-dorsal) or performing computed tomographic scans of lower jaws pertaining to many of the ceratopsoid species under study. The overall mediolateral shapes of the lower jaw appear quite similar across ceratopsoids, in terms of a laterally placed coronoid process and a ridge down the lateral surface of the jaw. Crushing can obscure these features in some cases, too. Therefore, 2D models offer the best chance to assemble a large sample.
Moreover, it is argued that 2D models are appropriate on the basis of an inferred orthopalinal power stroke for ceratopsid dinosaurs (Varriale, 2011;Mallon & Anderson, 2014a). Although torsional or transverse forces undoubtedly occurred in ceratopsid jaws, it is posited that these were far less important than the dorsoventral or craniocaudal movements implied by biomechanical studies. Essentially the jaw system is being modeled as a simple lever, and thus the 2D plane strain assumption reveals those strains that are likely of primary importance here.
The predentary was excluded from the landmark configuration because it is often missing, damaged or restored with plaster, generating an incomplete geometry. This assumption has been made in light of the available dataset.
To evaluate the stress on the lower jaws, a pressure load (i.e. a force per unit area) was applied, acting in the ventral direction, on the mid-point of the tooth row to simulate the chewing process. The bite force value calculated by Bell et al. (2009Bell et al. ( ) was used, 1131.06 N, and the same value was applied at the same point to all FE models for a comparative purpose. To ensure consistency in position, the rostral and caudal tips of the tooth row were identified from the original material, and the tooth row length calculated from photos using TPSDIG2. Then, two landmarks with a distance of 1 cm between each other on the mid-point of the tooth row were digitized, and the loading value applied on the same segment of each FE model. Isotropic material properties corresponding to bovine Haversian bone (Young's modulus: 10 GPa; Poisson's ratio: 0.41; Rayfield et al. 2001) were assumed for all FE models in this study.
Two distinct elastic constraints were placed on the lower jaw: the first was placed on the dorso-caudal margin of the coronoid process where the musculus pseudotemporalis (MPsT) superficialis and musculus adductor mandibulae externus (MAME) profundus insert (major muscles involved during mastication); and the second one was placed on the dorso-caudal margin of the surangular where musculus pterygoideus, musculus adductor mandibulae posterior, MPsT profundus and MAME superficialis and medialis insert (Ostrom, 1966;Bell et al. 2009;Holliday, 2009;Tanoue et al. 2009). All the remaining portions of the boundary are free. Figure 4 shows the 2D geometry, the portion of the boundary with the applied biting load and the portion with the elastic constraints. Each specimen underwent the same numerical experiment (loading and boundary conditions) using COMSOL MULTIPHSICS 5.0 software.
The lower jaw geometries were rescaled to the same unit size to exclude any scale effect among ceratopsoid species in assessing the mechanical performance variables. Due to the comparative approach applied, scaling has no direct effect on the results. This means that using a different bite force returns just differently scaled results (Piras et al. 2013).
The global von Mises stress was chosen, calculated on the whole lower jaw area, as an effective structural variable to be used in the subsequent comparative analyses.

Linear models and comparative methods
The potential effects of shared ancestry on the interspecific phenotypic variables were accounted for, and significant phylogenetic signal in size (CS), medio-lateral thickness and stress (von Mises) were tested for using the function phylosig() in 'phytools' R package (Revell, 2012). Shape consists of multivariate data, so the function physignal() in 'geomorph' R package (Adams & Otarola-Castillo, 2013) was used to estimate the degree of phylogenetic signal present in shape data for a given phylogenetic tree.
Then, differences in shape, size, thickness and stress between four distinct groups (non-triceratopsin Chasmosaurinae, Centrosaurini, Triceratopsini and Pachyrhinosaurini) were tested for. Pair-wise comparisons on size, thickness and stress were performed using non-parametric permuted ANOVA (PERANOVA), and on shape using nonparametric permuted MANOVA (PERMANOVA). Those tests were performed using the function adonis() of the 'vegan' R package (Oksanen et al. 2011), which manages unequal sample sizes (1000 permutations).
The relationships between shape (as dependent; average values per species), size, thickness and stress (as independent; average values per species) were accounted for by using an ordinary least square (OLS) linear regression in order to explore the influence of allometry and structural performance in the dataset. Being shape multivariate, canonical correlation analysis (CCA) was used, just for the sake of visualization, in order to plot the regression score representing shape against the independent variables, that is size (CS) or stress (von Mises values), respectively. Moreover, to ensure that those associations were not primarily the results of a phylogenetic effect, phylogenetic generalized least squares (PGLS; Rohlf, 2001Rohlf, , 2006 regression was applied. This method is equivalent to the phylogenetic independent contrasts (PICs) proposed by Felsenstein (1985). Data are rarely independent in biology as assumed in standard regressions. Taking phylogeny into account in regressions permits a reduction in the standard error due to a shared ancestry, over a model that neglects phylogeny (Felsenstein, 1985;Garland, 1992;Garland et al. 2005). Moreover, the influence of the maximum thickness of the lower jaw on stress variables was investigated, adopting an OLS model, to take into account if thickness variation affects the structural performance of the lower jaws. Linear measurements (cm) were taken of medio-lateral thickness at the base of the coronoid process from scaled photographs of each specimen. The average thickness values per-species were used in the analyses performed below. Unfortunately, there are no measurements of medio-lateral thickness for a few species in the sample (Achelousaurus horneri, Einiosaurus procurvicornis, Ojoceratops fowleri and Coahuilaceratops magnacuerna).

Lower jaw shape variation within Ceratopsoidea
The first four PCs of the bgPCA, performed on the lower jaws, in lateral view, explain collectively 98% of total shape variance. Figure 5a shows the relationship between PC1 (49.77% of the total shape variance) and PC2 (23.15% of the total shape variance), and Fig. 5b shows the relationship between PC1 and PC3 (corresponding to 19.57% of the total shape variation).
Negative PC1 values are associated with a lower jaw having a long, straight and massive dentary, a narrow rostral portion of the dentary, dorsally elongated coronoid process, large coronoid and short and low angular. This morphological arrangement is typical of triceratopsins. Positive PC1 values are associated with a shorter and dorso-ventrally thicker dentary, having a wider rostral portion of the dentary, a slightly concave profile of the dentary ventral edge, a less elongated and smaller coronoid process, and a larger angular.
At positive PC2 values the lower jaw is long, straight and robust, having a large coronoid, whereas at negative PC2 values the lower jaw possesses a smaller coronoid, a more massive dentary, having a concave ventral edge, a narrow rostral portion of the dentary, and a dorsal portion of surangular shifted dorso-rostrally.
The lower jaw has a slender and long dentary with a small coronoid process at positive PC3 values, whereas it is dorso-ventrally thicker and massive with a large coronoid process at negative PC3 values.
In summary, triceratopsins vary mainly along the negative PC1 values. The non-triceratopsin chasmosaurines Arrhinoceratops and Anchiceratops lie close to triceratopsins at negative PC1 values. Pachyrhinosaurini, Centrosaurini and the other non-triceratopsin chasmosaurines lie close to each other at positive PC1 values. The basal ceratopsoid Zuniceratops christopheri and the basal centrosaurine Avaceratops lammersi appear distinct with respect to other ceratopsid lower jaws. Figure S1 shows the 3D relationship between PC1, PC2 and PC3.

Allometric shape variation
An OLS regression between shape (dependent variable) and size (independent variable, quantified by CS) was performed for the entire sample to explore evolutionary allometry within Ceratopsoidea. Table 1 reports R 2 -and P-values for the pooled dataset, and for the individual clades under investigation. PGLS was performed only on the pooled dataset. The linear regression performed on the pooled dataset is not significant, and neither are the linear regressions performed within each clade, with the exception of Pachyrhinosaurini. PGLS results instead show a significant relationship between shape and size (Table 1). These results suggest the presence of evolutionary allometry within Ceratopsoidea. Figure 6 shows the relationship, via CCA, between the lower jaw shape changes, in lateral view, and CS. At low CS values, the lower jaw bears a relatively large and hookshaped coronoid, with a long, straight and massive dentary. At high CS values, the lower jaw is long and dorso-ventrally thicker, having a slightly concave ventral edge and bearing a slightly smaller and rounded coronoid process. Table 1 reports the P-values of the OLS regression between shape and von Mises stress, highlighting a significant result for the pooled sample under investigation. By contrast, all relationships tested within each clade are non-significant, with the exception of the non-triceratopsin chasmosaurines (P = 0.005). When accounting for phylogeny (PGLS), the relationship between shape and stress is statistically significant. Similar results are found when regressing shape with thickness, with the sole exception of the pooled sample OLS regression, which is not significant at the current cut-off level. Figure 7 shows the relationship, via CCA, between lower jaw shape and stress variable and the associated shape changes. At high von Mises values, the lower jaw has a long and more slender dentary, a moderately small angular, and a dorsally elongated coronoid process. A lower jaw bearing a thicker and massive dentary, a moderately large angular and a shorter hooked coronoid process is associated with low stress values.

Structural behavior of lower jaws
When regressing thickness and von Mises stress, the relationship is not significant in the pooled sample (R 2 = 0.053, P = 0.376), suggesting that thickness variation is not correlated with von Mises stress. The taxa with highly stressed morphologies do not possess a lower jaw with similar medio-lateral thickness (Fig. 8).

Linear models and comparative methods
The pair-wise PERMANOVA performed on shape variables between groups highlights significant differences between groups. Triceratopsins have a different lower jaw morphology relatively to that of Pachyrhinosaurini, Centrosaurini and non-triceratopsin Chasmosaurinae. Other groups have similar lower jaw shapes (Table 2).
Pair-wise PERANOVAs show a clear difference in size between Triceratopsini and other clades. Non-triceratopsin chasmosaurines show significant differences in size relative to Centrosaurini and Pachyrhinosaurini. The clade Table 1 OLS models for pooled dataset and individual clades for shape-size, shape-von Mises stress and shape-thickness relationships, and PGLS models for pooled sample. Pachyrhinosaurini has a size value similar to that seen in Centrosaurini (Table 3). Pair-wise PERANOVAs on thickness show Triceratopsini as significantly different from Centrosaurini. This latter clade results are distinct from those for non-triceratopsin chasmosaurines (Table 3). Concerning stress, triceratopsins show significant differences from Centrosaurini and Pachyrhinosaurini, whereas other clades share similar stress values (Table 3).

Discussion
Despite recent overviews suggesting no lower jaw differences between clades within Ceratopsidae (Dodson et al. 2004;Makovicky, 2012), the current investigation identifies some subtle differences of phylogenetic relevance, while confirming other areas of morphological conservatism. As such, this presents another excellent case study of how a GM approach can highlight features sometimes missed in subjective visual examination.
Centrosaurini (Centrosaurus apertus + Coronosaurus brinkmani + Styracosaurus albertensis) has a similar lower jaw morphology to that of Pachyrhinosaurini (derived centrosaurines) and non-triceratopsin Chasmosaurinae, as  highlighted by the results of the pair-wise PERMANOVAs ( Table 2). All of these taxa lie close to each other at positive PC1 values (Fig. 5a,b). Therefore, all centrosaurines investigated in this study possess similar lower jaw shapes.
The only exception is potentially represented by Avaceratops lammersi, a basal centrosaurine from the Late Campanian Judith River Formation of Montana (Dodson, 1986). It seems to have a different lower jaw shape from other   centrosaurines and partially similar to that of Anchiceratops and Arrhinoceratops. It has a short coronoid process with a sharp forward hook, and short, robust and thick dentary with a gently convex ventral edge (Penkalski & Dodson, 1999). However, this taxon is represented by only one lower jaw in the dataset (Academy of National Science of Philadelphia, ANSP 15981) pertaining to a sub-adult specimen. Thus, new fossil material of adult specimens is needed to confirm this result. Centrosaurine lower jaw shape variation resembles nontriceratopsin chasmosaurine morphology. In addition, all non-triceratopsin chasmosaurines have similar lower jaws.
Zuniceratops christopheri is a Turonian ceratopsoid, from the Moreno Hill Formation of New Mexico (Wolfe & Kirkland, 1998), with a different shape relative to that of ceratopsids. The dentary, bearing single-rooted teeth, is long, slender and possesses an elongated and thin coronoid process lacking the distinct 'hook' of ceratopsids (Wolfe et al. 2010). Unfortunately, Zuniceratops is represented by a single lower jaw (holotype) in the sample, and even if the referred material (Arizona Museum of Natural History, MSN P3201 and MSN P3202, not included in this work due to poor preservation) has similar morphology, new fossil material is needed to provide new insights on the non-ceratopsid condition of its lower jaw.
Triceratopsini (derived chasmosaurines) represents a clade having a lower jaw shape different from that of other ceratopsids. Triceratopsins are distinct in having a relatively thinner and longer dentary, with a straight ventral edge, a larger and more elongated coronoid process, an upper contact with the predentary shifted dorso-rostrally, and a small and dorso-ventrally compressed angular (Fig. 9). Other groups lie in a different position in the morphospace, suggesting a distinct lower jaw shape from triceratopsins. The only possible exceptions are represented by Anchiceratops ornatus and Arrhinoceratops brachyops, which lie close to the triceratopsin morphospace (Fig. 5a,b). Those two non-triceratopsin chasmosaurines represent the sister taxa of Triceratopsini in the current phylogeny (Fig. 3), and thus they share a somewhat similar lower jaw shape with triceratopsins (Fig. S1). By contrast, Mallon et al. (2011) argued that Anchiceratops shares more synapomorphies with Chasmosaurus than Triceratops, whereas Arrhinoceratops is much closer to the more derived chasmosaurines (Triceratopsini). Similarly, Mallon et al. (2014) highlighted potential poor resolution in the chasmosaurine phylogeny, in contrast with other hypotheses Mallon et al. 2011;Wick & Lehman, 2013;Longrich, 2014). Recently, Brown & Henderson (2015) proposed a new cladistic analysis, along with a description of the new chasmosaurine Regaliceratops peterhewsi, that shows Anchiceratops as sister to Arrhinoceratops and both taxa basal to Triceratopsini group. At least for the current sample at hand, the result of the physignal() function highlights a strong phylogenetic influence in the shape, therefore indicating related species have similar morphologies.
Size variation is phylogenetically structured, and it presents a discriminant trait within Ceratopsidae. Triceratopsini is the clade with the largest size values with respect to other groups (Fig. 6). Triceratops and its closest relatives represent the largest ceratopsids in the evolutionary history of Ceratopsia (Dodson et al. 2004;Brusatte, 2012). Non-triceratopsin chasmosaurines are larger than all centrosaurines (Pachyrhinosaurini and Centrosaurini), even if Utahceratops and Pentaceratops possess a similar size to that of triceratopsins. By contrast, Pachyrhinosaurini and Centrosaurini share similar size values (Table 3). Moreover, a strong evolutionary allometry results from this study (shape-size relationship, PGLS P = 0.006). However, when exploring the shape-size relationship within each clade, only Pachyrhinosaurini show a significant result (P = 0.04; Table 1). However, the large pachyrhinosaurin Pachyrhinosaurus sp., from the Late Campanian Dinosaur Park Formation of Alberta, has an unclear specific status. It lacks the portions of the frill with taxonomically diagnostic ornamentation, preventing a precise identification (Ryan et al. 2010). When excluding this specimen from the dataset, no allometric signal appears within Pachyrhinosaurini (P = 0.33).
Investigating the structural performance of the lower jaws by means of 2D FEA reveals that mechanical performance is conservative within Ceratopsidae (Fig. 7); this is unsurprising given the overall conservation of morphology. Triceratopsins (derived chasmosaurines) represent a clade within Ceratopsidae that differs in stress from Centrosaurini and Pachyrhinosaurini. No significant differences on stress values have been highlighted among other clades (Table 3), and the structural performance is phylogenetically constrained. Moreover, thickness variation does not significantly influence the structural performance within Ceratopsidae.
The elongation of the coronoid process associated with a longer and leaner dentary is characterized by high stress values, whereas a massive dentary and a shorter coronoid process is associated with a low stress framework (Fig. 7). However, lower jaw shape variation within Ceratopsidae is not associated with significant stress differences among groups (Table 3). The only exception is represented by Triceratopsini, in being relatively more stressed than Centrosaurini and Pachyrhinosaurini. A longer dentary, and a larger and more elongated coronoid process, typical of triceratopsins, could have produced a lower jaw morphology much more subject to stress. A long coronoid process increases the length of the moment arm (Fig. 1), thus potentially generating a higher bite force. Nevertheless, this statement appears counterintuitive; generating high bite forces, despite the mechanically stressed jaws, would overload the jaw itself. This could have occurred along the caudalmost part of the tooth row only (not investigated here), where the lower jaw acts as a second class lever (advantageous lever; Tanoue et al. 2009;Mallon & Anderson, 2015). Mesial to the coronoid process, the lower jaw acts as a third class lever (disadvantageous lever; Ostrom, 1966;Tanoue et al. 2009), therefore a longer dentary (i.e. longer lever) would have produced a comparative mechanical disadvantage (Mallon & Anderson, 2015). In addition, larger size requires that the coronoid process is larger as well, providing a wider muscular attachment area (for MPsT superficialis and MAME profundus) for triceratopsins when compared with centrosaurin and pachyrhinosaurin lower jaws. Therefore, triceratopsins could probably exert a relatively higher force by the jaw adductors (Ostrom, 1964) to compensate for the mechanical disadvantage in having a longer dentary. This higher magnitude of the applied force does not necessarily imply a greater bite force when chewing (Fig. 10).
Several anatomical characters might be evolved to compensate a mechanical disadvantage when chewing in ceratopsids. A prominent bony ridge that protrudes from the lateral surface of the dentary at the base of the coronoid process in some ceratopsids seems to be evolved to strengthen the lower jaw and avoid breakage (Bell et al. 2009;Mallon et al. 2011Mallon et al. , 2014. Centrosaurins, such as C. apertus and S. albertensis, and pachyrhinosaurins, such as A. horneri and P. lakustai, possess this structure on the lower jaw (Fig. 9). On the other hand, among non-triceratopsin chasmosaurines, only A. brachyops and A. ornatus appear to have this lower jaw feature (Mallon et al. 2011(Mallon et al. , 2014. Triceratopsins seem to not possess this prominent ridge (Fig. 9), therefore potentially possessing a lower jaw much more subjected to stress.
Thickness is not correlated with the structural behavior of lower jaws, suggesting that the increase of thickness within Ceratopsidae may not be a mechanical response to strengthen the structure (Fig. 8).
In the current FE experiments, it was not possible to consider the role of the enamel in ceratopsid teeth. The lingual enamel and dentine provide a tough layer on tooth useful to strengthen the dental structure and allowing to eat a tougher fodder (Hwang, 2005(Hwang, , 2011. Recent assessments highlighted that ortho-and vaso-dentine served to reduce the contact areas between maxillary and dentary teeth, thus reducing the friction created between them and decreasing the amount of force generated when chewing (Kay, 2014). Bell et al. (2009) pointed out that enamel acted to prevent transmission of stress into the dentine. However, modeling the occurrence of enamel layer in the current FE models, representing teeth with different mechanical properties, would return essentially the same results, because the loading effects would decay very fast over a very short distance in the neighborhood of the loaded area.
Mallon & Anderson (2014a) investigated tooth morphology and wear of the Campanian ceratopsids from the Dinosaur Park Formation (Alberta, Canada), and suggested that centrosaurines ate a more abrasive food than chasmosaurines. The current results partially confirm those findings in having Chasmosaurus belli more stressed than other chasmosaurines and centrosaurines from the Dinosaur Park Formation. By contrast, the other Campanian ceratopsids from the Dinosaur Park Formation share similar stress values. However, triceratopsins could reflect a similar pattern to that highlighted by Mallon & Anderson (2014a).
Triceratopsini lived in the North American landmass during the Late Maastrichtian, when the Western Interior Seaway (WIS) had almost entirely retreated. By contrast, the Centrosaurini in the current sample lived in Alberta and Montana, during the Campanian stage, when Laramidia was separated from Appalachia by the WIS , whereas Pachyrhinosaurini were spread in the Northern Laramidia (Alberta and Alaska) during the Late Campanian and early Maastrichtian (Fiorillo & Tykoski, 2012).
Arens & Allen (2014) studied the flora composition from the Hell Creek Formation (where Triceratops is one of the most common dinosaur) and highlighted that angiosperms, such as 'Vitis' stantoni, Dryophyllum subfalcatum and 'Celastrus' taurinensis, were the most abundant plants, whereas gymnosperms were instead rare. Similarly,  provided evidence on the relative abundance of angiosperms rather than gymnosperms at the Cretaceous-Paleogene boundary in the Hell Creek Formation, even if angiosperms were characterized by a species decline during the Late Maastrichtian. Recent works (Lupia et al. 2000;Eaton & Kirkland, 2008) point out the wide diversity of angiosperms during the Maastrichtian stage.
By contrast, the Campanian stage was characterized by the gymnosperms as the most abundant plants rather than the flowering plants (angiosperms). In 2007, Chin studied coprolites from the Two Medicine Formation and reported wood as the main fodder for the hadrosaur Maiasaura with possibly conifers. Recent assessments of the floral assemblage of the Dinosaur Park Formation (Braman & Koppelhus, 2005;Koppelhus, 2005), where centrosaurines and chasmosaurines lived, provided evidence that gymnosperms were more abundant than angiosperms and represented the main dietary item (Currie et al. 1995;. Additionally, ferns and pteridophites appear to be abundant in the Late Campanian (~73 Ma) of Wyoming (Wing et al. 1993(Wing et al. , 2012. Those disparate scenarios, characterized by different floral associations (gymnosperms in the Campanian and angiosperms in the Maastrichtian) could have affected the ecology of those animals and thus be related in part to lower jaw shape. Usually gymnosperms represent a tougher fodder for the diet than angiosperms (Farlow, 1987;Friis et al. 2011). This issue requires deeper investigations in the future. Bell et al. (2009) found in their 3D FEA of C. apertus that its lower jaw experiences stress concentrations mainly along the leading edge of the coronoid process. Muscle force placed a high degree of stress on the coronoid process, but several morphological adaptations such as the longitudinal Fig. 10 Schematic diagram of a third-class lever lower jaw with two different lever lengths. R represents the resistance force, and it is the same in the two cases (when having two distinct lever lengths). m + a is the moment arm of the resistance force as seen in Centrosaurini and Pachyrhinosaurini, whereas m + a + a 0 is the moment arm of the resistance force in Triceratopsini. m 0 is the moment arm of the applied forces, and F represents the magnitude of the applied force (the typeletter size is related to the magnitude of the applied force by the musculature). Black arrows indicate the short-lever condition (non-triceratopsin ceratopsids), whereas the gray arrows indicate the long-lever condition (triceratopsins).
ridge of the dentary buttressed the structure. Similar results are reported by this study: ceratopsids report high stress values in the dentary (see Fig. S2), where the prominent bony ridge is located, suggesting, as previously stated, that this structure evolved as a mechanical response to strengthen the lower jaw during chewing, although it cannot be a priori excluded that high stressed areas in the current FE models could represent an artifact of the 2D approach.

Limitations of the current work
Given its 2D approach, the current results and interpretations should be treated with caution, and considered as hypotheses for future tests on biomechanical performance with more detailed 3D FE models. Three-dimensional models are needed to account for differences in material properties, medio-lateral thickness of the lower jaw and possible bending moments in the transverse plane (Rayfield, 2005(Rayfield, , 2012, although it is rather difficult to apply this methodology to the sample presented here. Recent work has also highlighted the validity of the 2D method when exploring the mechanical performance of the theropod skulls, with results similar to those found in the 3D approach (Rayfield, 2012). In addition, medio-lateral thickness was included as a 3D variable. Therefore, it is concluded that the current results are reasonable within the limits of the assumptions.
Further studies should focus on bolstering the framework for studying the mechanical performance of the ceratopsid lower jaws. Dynamic models (rather than static), as well as 3D models, could prove useful information, particularly with the inclusion of the elements not considered here (e.g. articular, splenial and predentary). In addition, a complete lower jaw is not known for some species; thus, new fossil material will undoubtedly provide new information on the systematics, palaeoecology and evolution of biomechanical performance in this group in response to specific functional demands.

Conclusions
Triceratopsins (the most derived chasmosaurines) possess a lower jaw morphology distinct relative to that of other ceratopsids. Triceratopsins are characterized in having a robust and longer dentary, a larger and more elongated coronoid process, and a small and dorso-ventrally compressed angular. In addition, the members of this clade are also distinct from other ceratopsids in having the largest size within the sample. Other ceratopsid groups, such as Centrosaurini, Pachyrhinosaurini and non-triceratopsin Chasmosaurinae, have similar lower jaw shapes. Moreover, these morphological differences are accompanied by different structural performance of the lower jaws. Triceratopsins possess a more stressed lower jaw than centrosaurins and pachyrhinosaurins, possibly correlated with distinct dietary regimes.
The possible inclusion of Anchiceratops and Arrhinoceratops in the clade Triceratopsini appears to be supported here, even if this statement is just related to the lower jaw morphology. Mallon et al. (2011Mallon et al. ( , 2014 noted that the chasmosaurine phylogeny is still controversial and unresolved. A further investigation and comparison of cranial material of those taxa with that of triceratopsins could provide new insights on this issue. Unfortunately, several triceratopsins (i.e. Nedoceratops hatcheri, Eotriceratops xerinsularis and Torosaurus latus) are poorly represented by lower jaws in the fossil record. New data will be fruitful to provide information on the origin of the unique jaw shape and structural performance of triceratopsins, in particular, and ceratopsids in general.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Appendix S1. Detailed description of phylogenetic relationships and complete literature corpus upon which the ceratopsoid phylogeny was built. Table S1. List of Institutional abbreviations. Table S2. List of ceratopsoid material used for this study and references. Table S3. Landmark definitions for the lower jaw in lateral view (see Fig. 2). Table S4. List of ceratopsoid taxa, stratigraphic and geographical distribution, and references considered in this study. Table S5. Age and references for phylogenetic positions and stratigraphic ranges of nodes and OTUs on the phylogenetic tree depicted in Fig. 3. Fig. S1. Dynamic 3D plot of PCA of lower jaws in lateral view. Hulls represent morphospaces for Triceratopsini (black); Pachyrhinosaurini (red); Centrosaurini (green); and non-triceratopsin Chasmosaurinae (blue). The light blue point represents Avaceratops lammersi and the purple represents Zuniceratops christopheri. Points dimensions are proportional to specimen CS.