Received: 17 March 2020 | Revised: 18 April 2020 | Accepted: 20 April 2020 DOI: 10.1002/ece3.6361 ORIGINAL RESEARCH Modularity and heterochrony in the evolution of the ceratopsian dinosaur frill Albert Prieto-Márquez1,2 | Joan Garcia-Porta3,4 | Shantanu H. Joshi5 | Mark A. Norell6 | Peter J. Makovicky2,7 1 Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Abstract Barcelona, Barcelona, Spain The fossil record provides compelling examples of heterochrony at macroevolution- 2 Integrative Research Center, Field Museum ary scales such as the peramorphic giant antlers of the Irish elk. Heterochrony has of Natural History, Chicago, IL, USA 3 CREAF, Barcelona, Spain also been invoked in the evolution of the distinctive cranial frill of ceratopsian dino- 4 Department of Biology, Washington saurs such as Triceratops. Although ceratopsian frills vary in size, shape, and ornamen- University, St. Louis, MO, USA tation, quantitative analyses that would allow for testing hypotheses of heterochrony 5 Department of Neurology and Ahmanson are lacking. Here, we use geometric morphometrics to examine frill shape variation Lovelace Brain Mapping Center, University of California, Los Angeles, Los Angeles, CA, across ceratopsian diversity and within four species preserving growth series. We USA 6 then test whether the frill constitutes an evolvable module both across and within Division of Paleontology, American Museum of Natural History, New York, NY, species, and compare growth trajectories of taxa with ontogenetic growth series to USA identify heterochronic processes. Evolution of the ceratopsian frill consisted primarily 7 Department of Earth and Environmental of progressive expansion of its caudal and caudolateral margins, with morphospace Sciences, University of Minnesota, Minneapolis, MN, USA occupation following taxonomic groups. Although taphonomic distortion represents a complicating factor, our data support modularity both across and within species. Correspondence Albert Prieto-Márquez, Institut Català de Peramorphosis played an important role in frill evolution, with acceleration operating Paleontologia Miquel Crusafont, Universitat early in neoceratopsian evolution followed by progenesis in later diverging cornos- Autònoma de Barcelona, Carrer de l’Escola Industrial 23, 08201 Sabadell, Barcelona, aurian ceratopsians. Peramorphic evolution of the ceratopsian frill may have been Spain. facilitated by the decoupling of this structure from the jaw musculature, an inference Email: redshore@gmail.com that predicts an expansion of morphospace occupation and higher evolutionary rates Funding information among ceratopsids as indeed borne out by our data. However, denser sampling of the Generalitat de Catalunya; Field Museum of Natural History; National Science meager record of early-diverging taxa is required to test this further. Foundation, Grant/Award Number: FRES 1925884 and PLR 1341645 ; Ministry of KEYWORDS Economy, Industry and Competitivity of Spain, Grant/Award Number: FJCI-2014- dinosaur, evolution, heterochrony, modularity, morphometrics 20380 and RyC-2015-17388 1 | I NTRO D U C TI O N Megaloceros gigantaeus, with its giant antlers, is frequently cited as an iconic example of peramorphosis (Gould, 1974). Nonavian dino- The fossil record provides ample opportunities for the study of the saurs (hereafter dinosaurs), some of which reached extreme body role of heterochrony at macroevolutionary scales, due to the pres- sizes, and many of which exhibited bizarre anatomical traits, have ence of taxa with larger body sizes and more extreme morpholo- also been investigated with respect to heterochronic evolution gies than those known in extant biota. For example, the Irish elk (Long & McNamara, 1995; McNamara & Long, 2012). Many previous This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd 6288 |  www.ecolevol.org Ecology and Evolution. 2020;10:6288–6309. PRIETO-MÁRQUEZ et al. | 6289 studies on heterochrony in dinosaur evolution dealt with life-history Long & McNamara, 1995, 1997; McNamara, 2012; Tumarkin & parameters related to size and growth (Erickson & Druckenmiller, Dodson, 1998). These studies have only dealt with scant quantifi- 2011; Erickson et al., 2004; Lee & Werning, 2008; Weishampel & cation of either frill size or shape, however, and are based on simple Horner, 1994), but a growing body of work has examined the role linear measurements or qualitative observations. of heterochrony in the evolution of individual anatomical regions Here, we take a geometric morphometric approach to investi- and character systems of dinosaurs (e.g., Bhullar et al., 2012; Farke, gate patterns of morphological variation of the ceratopsian frill, Chok, Herrero, Scolieri, & Werning, 2013; Forth, Hedrick, & Ezcurra, bearing on three interrelated questions regarding the ontogeny and 2016; Guenther, 2009; Horner & Goodwin, 2006). Nevertheless, evolution of this structure: quantitative studies on heterochrony are still rare for this clade, de- spite a growing sample of species with increasing ontogenetic sam- 1. Does shape variation vary randomly or are patterns (e.g., phylo- pling and an improving understanding of their development. genetic; functional) evident in frill morphospaces? Which aspects Here, we examine the evolution of the hallmark cranial frill of frill structure exhibit the most variation across species and (Figure 1) of one of the most iconic dinosaur lineages, the Ceratopsia. within species with good ontogenetic sampling? Ceratopsians were a diverse (>70 spp.) ornithischian clade and are one 2. Does the frill constitute a module in the ceratopsian skull across of the best-sampled major clades of nonavian dinosaurs, with some both evolutionary and developmental scales? Heterochronic evo- species representing the most abundant taxa in their fossil depos- lution is thought to occur within semi-independent collections its (Dodson, Forster, & Sampson, 2004; Dodson, 2013). Originating of traits, or modules (Gerber & Hopkins, 2011; Goswami, 2007; in the Middle to Late Jurassic (Xu, Forster, Clark, & Mo, 2006; You Gould, 1977; Klingenberg & Marúgan-Lobon, 2013; Olson & & Dodson, 2004; Zhao, Zhengwu, & Xu, 1999), these animals were Rosell, 2006), so establishing whether the frill behaves as a mod- among the more dominant herbivores of Cretaceous terrestrial eco- ule across evolutionary and ontogenetic scales is prerequisite to systems in North America (Dodson et al., 2004; Farke, Maxwell, Cifelli, testing for heterochrony as a process for changes in frill shape. & Wedel, 2014) and Asia (Makovicky & Norell, 2006; Maryanska & 3. Is frill evolution driven by heterochronic processes such Osmolska, 1975; Xu et al., 2006; You, Li, Ji, Lamanna, & Dodson, 2005), as peramorphosis as proposed by previous authors (Long and are also reported from the Late Cretaceous of Europe (Lindgren & McNamara, 1995, 1997; McNamara, 2012; Tumarkin & et al., 2007; Ösi, Butler, & Weishampel, 2010). Members of this clade Dodson, 1998), and is it possible to characterize these processes had remarkably large heads (Sereno, Xijin, Brown, & Lin, 2007), laterally as well as their potential drivers? compressed turtle-like beaks, and shearing dentitions (Ostrom, 1964), which are thought to have contributed to their evolutionary success (Ostrom, 1966). Most notably, ceratopsians are instantly recogniz- able for their exaggerated cranial structures (Knell & Sampson, 2011; 1.1 | Institutional abbreviations Padian & Horner, 2011), particularly the prominent parieto-squamosal frill at the back of the skull, characteristic of the clade as a whole, AMNH, American Museum of Natural History, New York, USA; BMRP, and also facial and frill ornaments in the more derived Campano- Burpee Museum of Natural History, Rockford, Illinois, USA; CMN, Maastrichtian Ceratopsidae (Brown & Henderson, 2015; Dodson & Canadian Museum of Nature, Ottawa, Canada; FMNH, Field Museum Currie, 1990; Dodson et al., 2004; Sampson et al., 2010; Sampson, of Natural History, Chicago, USA; IVPP, Institute of Vertebrate Lund, Loewen, Farke, & Clayton, 2013). Paleontology and Paleoanthropology, Beijing, China; MOR, Museum of The ceratopsian frill is a neomorphic structure that overhangs the the Rockies, Bozeman, Montana, USA; MPC, Mongolian Paleontological occiput to varying degrees and is formed by the parietal along the Collection, Ulaanbaatar, Mongolia; MDP, Mongolian Dinosaur Project; midline with contributions of the squamosals along its lateral edges MPM, Milwaukee Public Museum, Milwaukee, Wisconsin, USA; (Dodson et al., 2004; Hatcher, Marsh, & Lull, 1907; Figure 1). In small ROM, Royal Ontario Museum, Toronto, Ontario, Canada; TMP, Royal bodied, early-diverging ceratopsians it is relatively short and narrow, Tyrrell Museum of Paleontology, Drumheller, Alberta, Canada; UCMP, but in large bodied (>1,000 kg) ceratopsids, the frill alone can be over a University of California Museum of Paleontology, Berkeley, California, meter in length and width and constitutes more than half the length of USA; UMNH, Natural History Museum of Utah, Salt Lake City, Utah, the skull (Dodson et al., 2004). It is perforated by a pair of fenestrae in USA; YPM, Yale Peabody Museum of Natural History, New Haven, most neoceratopsians. Although the frill margin is smooth or relatively Connecticut, USA. unadorned in most early-diverging species, the large and speciose cer- atopsids bear epiparietal and episquamosal ossifications that form a spectacular diversity of structures projecting from the frill, and that 2 | M ATE R I A L S A N D M E TH O DS distinguish individual species from one another (Figure 1). Unsurprisingly, heterochronic processes have been invoked 2.1 | Sampled taxa and specimens to explain this phylogenetic variation and its general relation- ship with the evolution of body size (~5 kg–10 tons), with various Our study sample consisted of 25 species representing all taxa authors interpreting frill evolution as a peramorphic process (e.g., for which reliable data on frill shape are available and nearly half 6290 | PRIETO-MÁRQUEZ et al. the known taxonomic diversity of Ceratopsia. The sampled spe- Table 1). We conducted two-dimensional geometric morphometric cies comprise three basal Ceratopsia, four basal Neoceratopsia, 10 analyses to quantify frill shape variation throughout ceratopsian chasmosaurine, and eight centrosaurine ceratopsids (Figure 1d and evolution. Although skulls are three-dimensional structures, the F I G U R E 1 Landmark configuration used in the Procrustes superimposition shape analysis. (a) Line drawing of the dorsocaudal view of the skull of the ceratopsid Centrosaurus apertus (CMN 8795). (b) Dorsocaudal view of CMN 8795. (c) Dorsocaudal view of the skull of the basal neoceratopsian Protoceratops andrewsi (FMNH PR 14045). (d) Simplified phylogeny of Ceratopsia with the species included in this study (Table 1). The tree is based on those of Xu et al. (2006), Farke et al. (2014), Longrich (2014), Brown and Henderson (2015), and Evans and Ryan (2015). Yellow stars and white circles correspond to landmarks and semilandmarks, respectively. Landmark/semilandmark numbering is in the order in which they were digitized and indicate: 1, tip of the rostral bone; 2, rostral end of the parietal at the sagittal plane of the skull; 3 and 10, rostrolateral margin of the supratemporal fenestra; 4 and 11, rostral-most margin of the infratemporal fenestra; 5 and 9, caudolateral end of the supratemporal fenestra; 6 and 8, squamosal–parietal contact at the caudolateral margin of the frill or skull; 7, mid-point of the caudal margin of the frill at the interparietal contact; 36 and 40, rostromedial corner of the parietal fenestra; 37 and 43, rostrolateral corner of the parietal fenestra; 38 and 42, caudolateral corner of the parietal fenestra; 39 and 41, caudomedial corner of the parietal fenestra; and 44 and 45, lateroventral tip of the jugal boss PRIETO-MÁRQUEZ et al. | 6291 TA B L E 1 Sample of ceratopsian specimens used as representatives of the main ceratopsian clades in the morphometric analyses Species (Abbreviation) Specimen Achelousaurus horneri (Ach) Centrosaurinae; Sampson, 1995, figure 5b, reconstruction based on MOR 845 and 591 Anchiceratops ornatus (Anc) Chasmosaurinae; Mallon, Holmes, Eberth, Ryan, & Anderson, 2011, figure 2, AMNH 5251 Arrhinoceratops brachyops (Arr) Chasmosaurinae; Mallon, Holmes, Anderson, Farke, & Evans, 2014, figure 4, ROM 796 Centrosaurus apertus (Cen) Centrosaurinae; CMN 8795 Chasmosaurus belli (Cen) Chasmosaurinae; CMN 2245 Diabloceratops eatoni (Dbl) Centrosaurinae; UMNH cast of UMNH VP 16699 Einiosaurus procurvicornis (Ein) Centrosaurinae; Sampson, 1995, figure 4b reconstruction based on MOR 456 Kosmoceratops richardsoni (Kos) Chasmosaurinae; UMMH cast UMNH VP 17000 Leptoceratops gracilis (Lep) Basal Neoceratopsia; Sternberg, 1951, Pl. XLVIII, CMN 8889 Liaoceratops yanzigouensis (Lia) Basal Neoceratopsia; IVPP V12738 Mojoceratops perifania (Moj) Chasmosaurinae; Longrich, 2010, figure 9, including TMP 1983.25.1 Nasutoceratops titusi (Nas) Centrosaurinae; UMNH cast of UMNH VP 16800 Pachyrhinosaurus lakustai (Pac) Centrosaurinae; UMNH cast of TMP 1986.55.258 Pentaceratops sternbergii (Pen) Chasmosaurinae; Forster, Sereno, Evans, & Rowe, 1993, figure 6, AMNH 6325 & MNA Pl.1747 Protoceratops andrewsi (Prt) Basal Neoceratopsia; FMNH PR 14045 Ukhaa Tolgod protoceratopsid gen. et sp. nov. (Ukt) Basal Neoceratopsia; IGM 100/1246 (Ukhaa Tolgod, Mongolia) Psittacosaurus lujiatunensis (Psm) Basal Ceratopsia; FMNH PR 2588 (cast) Psittacosaurus mongoliensis (Psl) Basal Ceratopsia; AMNH 6254 Regaliceratops peterhewsi (Reg) Chasmosaurinae; TMP 2005.055.0001 Styracosaurus albertensis (Sty) Centrosaurinae; CMN 344 Triceratops horridus (Tri) Chasmosaurinae; YPM 1823 (Hatcher et al., 1907, Pl. XXVIII) Utahceratops gettyi (Uta) Chasmosaurinae; Sampson et al., 2010, figure 4, including UMNH VP 16784 Vagaceratops irvinensis (Vag) Chasmosaurinae; Holmes, Forster, Ryan, & Shepherd, 2001, figure 3, CMN 41357 Wendiceratops pinhornensis (Wen) Centrosaurinae; Evans & Ryan, 2015, figure 15, including TMP 2011.051.0009 Yinlong downsi (Yin) Basal Ceratopsia; IVPP V14530 Note: The asterisk denotes those specimens used only in the SRVF analysis. The abbreviations noted between brackets correspond to those used in Figures 2–7. two-dimensional analysis is an appropriate simplification in that Goodwin (2006: p. 2757). With respect to Protoceratops, our rep- it suffices to capture most of the shape of the ceratopsian frill. resentative skull in the all-species sample (FMNH PR 14045) is This is because the frill is a dorsoventrally compressed structure, larger than two skulls (AMNH 6418; AMNH 6466) that have been and as such, its morphology is mostly contained and exposed on aged histologically to be 17 and 19 years and interpreted as adults a dorsal plane of the skull that is facing anterodorsally (Figure 1). in an unpublished study (Makovicky, Sadleir, Dodson, Erickson, & Furthermore, and from a practical point of view, 3D scanning was Norell, 2007). Our representative specimen for P. lujitunensis is of a not an available option and it was not possible to photograph many size (skull length = 19 cm) that overlaps with the oldest sampled in- of the specimens of our sample in all the required views necessary dividuals in a histological study of that species (Erickson, Makovicky, to construct a 3D photogrammetric model of the skull. In particu- Inouye, Keqin, & Zhou, 2009). Thus, for taxa with available skele- lar, it was not uncommon to encounter specimens in exhibit mounts tochronological or histological age data, our choice to use the larg- obscuring some views (e.g., ventral or posterior) of the skull. For the est available specimen is selecting adult or large subadult specimens analyses at interspecific level, each species was represented by one for comparisons. We also selected those specimens that preserve a individual. As a proxy for adulthood, we selected a specimen reach- complete or nearly complete and undistorted frill on at least one side ing the largest (or at least within the range of the largest) known of the skull; taxa lacking sufficiently complete frills were excluded. size for the species. For Triceratops, the select specimen (YPM 1823) The skulls of the sampled species were photographed in dorso- exhibits anatomical criteria such as anterior postorbital horn orienta- caudal view, orthogonal to the dorsal surface of the frill. Particular tion and flattened episquamosal and epiparietal ossifications identi- attention was paid to maintaining that orientation as strictly as fied as representative of adulthood by Horner and Goodwin (2006); possible, given that deviations from a dorsocaudal view can dis- these authors also regard all YPM specimens as adults (Horner and tort the contour of the frills and, potentially, the results of the 6292 | PRIETO-MÁRQUEZ et al. morphometric analyses. In the case of specimens preserving only the epiparietal and episquamosal ornaments are still debated (e.g., one side of the frill and/or displaying substantial taphonomic distor- Brown & Henderson, 2015; Farke et al., 2011; Mallon, Ott, Larson, tion, we symmetrically mirrored the better-preserved side around Iuliano, & Evans, 2016). Exclusion of the ceratopsid ornamentation the sagittal plane of the skull. This correction was applied to one was accomplished by linking the bases of all ornamental epiparietal Liaoceratops yanzigouensis specimen (IVPP V12738), four P. andrewsi and episquamosal structures with straight lines orthogonal to their specimens (AMNH 6425, 6433 and 6439; MPC D100/502), and the long axes, resulting in a polygon following the boundary of the frills ceratopsids Regaliceratops peterhewsi (TMP 2005.055.0001) and (Figure 1b). Landmarks and semilandmarks were digitized on photo- Pachyrhinosaurus lakustai (cast of TMP 1986.55.258). For the anal- graphs of the ceratopsian skulls using the program tpsUtil version yses across species, we sampled only one specimen per species in 1.68 (Rohlf, 2016a). Both sides of the skull were digitized to bet- order to prevent biasing the mean shape of the sample by including ter allow for tests of modularity. The files with the landmark coor- unequal numbers of specimens per taxon (a common situation in the dinates (Appendices S1–S6) were created in tpsDig2 version 2.26 dinosaurian fossil record). (Rohlf, 2016b). For ease of interpretation of morphospace plots, we color coded The landmark configurations of the skull sample were super- the sampled taxa in four groupings: basal ceratopsians comprising imposed by means of the generalized Procrustes procedure, which non-neoceratopsian taxa with short, incipiently developed and un- eliminates all nongeometric information after scaling, rotation, and fenestrated frills (gray); basal noeceratopsians, representing the translation relative to the grand mean (Bookstein, 1991; Klingenberg, nonceratopsid neoceratopsian part of our sample with larger frills Barluenga, & Meyer, 2002). Optimal superimposition of configu- with unornamented margins (orange); and the two main ceratopsid rations of landmark coordinates used least-squares estimation of lineages Chasmosaurinae (blue) and Centrosaurinae (red), both of translation, rotation, and other scaling parameters (Slice, 2007). The which have heavily ornamented frills. Although the primary function generalized Procrustes procedure was conducted using the function of these designations is to facilitate the interpretation of morpho- “gpagen” in the R package Geomorph 3.0.4 (Adams and Otarola- space plots, they correspond approximately to previously proposed Castillo, 2013). categories of frill evolution (e.g., Ostrom, 1966). The resulting shape variables were used to perform a phyloge- netic principal components analysis (pPCA). This is a widely used dimension reduction technique that represents sample variation 2.2 | Shape analysis using procrustes on a reduced number of axes, taking into account phylogenetic in- superimposition formation (Revell, 2009). We performed the pPCA using a correla- tion structure derived from a Brownian motion model by means of The ceratopsian frill and skull table were digitized using 21 homol- the function “phyl.pca” in the R package Phytools 0.6-0.0 (Revell, ogous landmarks (Figure 1b,c). Preservational issues and fusion of 2012). Because there is no published phylogenetic analysis that cranial bones prevented the identification of additional landmarks fully resolves the relationships of all known species of basal cera- of the skull roof and lateral portions of the facial skeleton exposed topsians, centrosaurines, and chasmosaurines, we constructed an in dorsocaudal view. We added 38 evenly spaced semilandmarks informal supertree (Figure 1d) that combined the topologies of Farke (Figure 1b,c) to more completely capture the shape of the frill and et al. (2014, resolving basal neoceratopsian relationships), Longrich skull. In most neoceratopsians, the frill is perforated by a pair of (2014, focusing on chasmosaurines), and Evans and Ryan (2015, fo- parietal fenestrae, which were digitized using four landmarks and cusing on centrosaurines); the phylogeny of Brown and Henderson included in our analysis. However, the parietal fenestra is absent in (2015) was used to position Regaliceratops. The tree was trimmed to non-neoceratopsian Ceratopsia, Leptoceratops, hatchling specimens only include the taxa in our GM sample and calibrated using the first of Protoceratops, and in all but the oldest individuals of Triceratops and last appearance datum (Appendix S7) for each species via imple- (Scannella & Horner, 2010; Torosaurus of some authors). In order mentation of the R package “strap” (Bell & Lloyd, 2015). We used the to include these taxa in our study, we also used an iteration of the function DatePhylo in the package strap to estimate branch lengths Procrustes superimposition analysis in which the parietal fenestra using the method of Brusatte, Benton, Ruta, and Lloyd (2008) for use was excluded from the landmark configuration. in the pPCA and other comparative analyses. Within Neoceratopsia, the more derived Ceratopsidae exhibit Morphospace occupation by the different ceratopsian clades and hypertrophied frill ornamentation, often in the form of spikes, grades was visualized by plotting the shapes under consideration on hooks, and other types of processes of greatly varying morphologies combinations of the first three pPCA axes. Shape variation contained (e.g., Figure 1a). Notable examples include the long parietal spikes in each plot was represented by estimating the mean landmark coor- of Styracosaurus albertensis (Ryan, Holmes, & Russell, 2007) and dinates of all species and calculating the shape differences existing Diabloceratops eatoni (Kirkland & DeBlieux, 2010). We excluded such between the mean and each of the species lying at the extremes ornamental structures from the landmark configuration (Figure 1b) of each pPC axis (functions “mshape” and “plotRefToTarget” in the in order to compare the morphology of the landmark configura- R package Geomorph 3.0.4 (Adams and Otaróla-Castillo, 2013). tion of the margins of the ceratopsid frills to that of unornamented Testing for significance in the morphological difference among the ceratopsid outgroups, and also because the homologies between two grades and five clades of ceratopsians under consideration (i.e., PRIETO-MÁRQUEZ et al. | 6293 basal ceratopsians, Neocertopsia, nonceratopsid neoceratopsians, the parietal fenestra, were conducted for the clades Neoceratopsia, Carnosauria, Ceratopsia, Centrosaurinae, and Chasmosaurinae) was Ceratopsidae, and Chasmosaurinae to allow for inclusion of the accomplished by means of pairwise permutational ANOVAs along chasmosaurine Triceratops horridus. each of the first three pPC axes, which relied on 1,000 randomiza- Intraspecifically, we assessed frill modularity in growth series tions to compute p-values (this was calculated with R custom scripts of Liaoceratops yanzigouensis, Protoceratops andrewsi, a new unde- relying on the function “lm”; Tables 4 and 5). scribed protoceratopsid from Ukhaa Tolgod, and Triceratops horridus. Only five exemplars of L. yanzigouensis were available, but these spanned a substantial range of sizes (Table 2). P. andrewsi is known 2.3 | Evaluation of the frill modularity hypothesis from numerous specimens, many of which are exceptionally pre- served, including a nearly complete developmental series ranging Certain traits in organisms may become compartmentalized into from hatchlings through full-grown adults (Brown & Schlaikjer, 1940; modules, clusters of characters that evolve in a highly correlated Dodson, 1976; Fastovsky et al., 2011; Hone, Wood, & Knell, 2016; manner within a cluster, but relatively independently from charac- Lull, 1908; Russell, 1935).The growth series of 21 P. andrewsi and the ters in other clusters (Klingenberg, 2008). In contrast to their strong nine specimens of the new Ukhaa Tolgod protoceratopsid spanned internal integration, intermodular correlations are hypothesized much of the development of these species, with specimens rang- to be weaker (Adams, 2016; Klingenberg, 2009; Klingenberg & ing from a few dozens of mm to more than 300 mm in basal skull Marúgan-Lobon, 2013). Because the ceratopsian frill is a neomor- length (Table 2). T. horridus is the only available ceratopsid preserv- phic structure among archosaurs and exhibits considerable varia- ing an ontogenetic series of skulls with frills spanning much of its tion in shape, size, and ornamentation among species that otherwise development (Goodwin, Clemens, Horner, & Padian, 2006; Horner exhibit reltively minor differences in postcranial and dental charac- & Goodwin, 2006, 2008), although the extensive MOR sample of ters, we hypothesize that it may have evolved in modular fashion, Triceratops was not available to us for study. We followed Scannella that is, relatively independent of the rest of the skull. Accordingly, and Horner (2010) in considering Torosaurus a junior synonym of we evaluated the a priori hypothesis of frill modularity by asking Triceratops in our primary analysis. The T. horridus ontogenetic series whether there is stronger covariation among the landmarks defining was less complete, ranging from 450 mm to nearly 1,150 mm in basal the frill than between any subset of them and subsets of landmarks skull length (n = 8, Table 3). The T. horridus specimens UCMP 154452 from other skull partitions. We employed the covariance ratio (CR; and MPM VP6841 were excluded from the modularity analyses, be- Adams, 2016) as a measure of the covariation between subsets of cause they do not preserve the rostrum. landmarks. The CR measures the overall covariation between land- Positive modularity results could result from including a greater mark partitions (or modules) relative to the overall covariation within number of evenly spaced semilandmarks than primary landmarks partitions (Adams, 2016), resulting in a coefficient that ranges from (Klingenberg, pers. Communication). To ensure that positive and sta- zero to positive values. CR values ranging from zero to one are indic- tistically significant results were not just an artifact of such a bias, we ative of modularity, whereas values above one characterize greater devised a sensitivity analysis in which all of the above CR analyses covariation between partitions than within them (Adams, 2016). were repeated with varying numbers (n = 10, 20, 28) of semiland- Although higher within-partition covariation does not necessarily marks removed. imply modularity, it is a testable prediction of the modularity hy- pothesis (Klingenberg, 2009). This reasoning also provides a basis for accepting the null hypothesis: If covariation between the frill parti- 2.4 | Heterochrony testing tion and other partitions of the skull is not weaker than among most or all (e.g., >95%) of the other partitions considered, the hypothesis We considered the definition of heterochrony as changes in the of frill modularity can be rejected. Notably, the CR coefficient is not relative timing and rate of developmental events between ances- affected by sample size nor the number of variables (Adams, 2016). tor and descendant species, as well as the onset and rate at which The CR analysis was implemented in the R package Geomorph the processes responsible for such changes take place (Alberch, (Adams and Otárola-Castillo, 2013). The analysis compared the CR Gould, Oster, & Wake, 1979; Gould, 1977, 1992; Klingenberg, 1998; coefficient of the landmark partition corresponding to the frill with McKinney & McNamara, 1991). Typically, heterochronic processes that of 10,000 random alternative subsets of landmarks across the are identified by analyzing allometric variation between taxa rela- whole skull. The output of this analysis was summarized as a his- tive to developmental age (Meyer, 1987; Ryan and Semlitsch, 1998; togram that displays the distribution of CR coefficients for all the Webster and Zelditch, 2005). In the absence of absolute specimen randomized partitions in the skull, thus providing a means for sta- ages, we follow more recent studies (Drake and Klingenberg, 2007; tistically evaluating the robustness of the observed CR coefficient. Forth et al., 2016; Klingenberg & Marugán-Lobón, 2013) in using We used the landmark configuration from the Procrustes analyses centroid size as the independent variable for investigating ontoge- described above to test for the modularity of the frill at various netic allometry in the four species. We performed an initial regres- evolutionary and ontogenetic levels. Two variations of the analyses, sion of shape scores on log-transformed centroid size for all terminals that is, including and excluding those landmarks corresponding to to determine whether there was evidence of allometry (Drake and 6294 | PRIETO-MÁRQUEZ et al. TA B L E 2 Sample of basal Basal skull Frill Relative neoceratopsian and protoceratopsid Species Specimen length (mm) inclination frill length specimens used in this study L. yanzigouensis Uncatalogued #1 66 — — L. yanzigouensis IVPP V12633 93 — — L. yanzigouensis IVPP V12738 99 — — L yanzigouensis Uncatalogued #2 137 — — L. yanzigouensis Uncatalogued #3 165 — — P. andrewsi MPC-D100/530 28 — — P. andrewsi AMNH 6419 76 22° 0.61 P. andrewsi AMNH 6431 150 22° 1.04 P. andrewsi AMNH 6409 191 43° 1.02 P. andrewsi AMNH 6444 210 27° ? P. andrewsi AMNH 6408 235 27° 0.65 P. andrewsi MPC-D100/502 250 31° 1.22 P. andrewsi AMNH 6433 261 34° 0.69 P. andrewsi AMNH 6429 269 33° 0.91 P. andrewsi AMNH 6439 271 49° 0.73 P. andrewsi MDP 530 297 ? 1.06 P. andrewsi AMNH 6425 313 43° 0.86 P. andrewsi AMNH 6413 314 40° 0.73 P. andrewsi FMNH PR 14045 330 37° 1.05 P. andrewsi AMNH 6414 341 57° 0.67 P. andrewsi AMNH 6466 357 36° 0.69 Ukhaa Tolgod IGM 100/3655 24 — — protoceratopsid Ukhaa Tolgod IGM 100/1008 28 — — protoceratopsid Ukhaa Tolgod IGM 100/10020 37 — — protoceratopsid Ukhaa Tolgod IGM 100/1013 47 — — protoceratopsid Ukhaa Tolgod IGM 100/1019 69 5° — protoceratopsid Ukhaa Tolgod IGM 100/3658 146 18° 0.90 protoceratopsid Ukhaa Tolgod IGM 100/3596 220 16° 0.96 protoceratopsid Ukhaa Tolgod IGM 100/3501 323 12° 1.01 protoceratopsid Ukhaa Tolgod IGM 100/1246 420 32° 1.02 protoceratopsid Note: Measurements for Protoceratops andrewsi were taken from Dodson (1976, table 1), except for FMNH PR 14045, MPC D100/502, and MDP 530, which were measured for this study. Basal skull length was measured from the tip of the rostral to the caudal end of the quadrate. Frill inclination was measured as the angle between the sagittal crest and the maxillary tooth row. Relative frill length is the maximum length of the frill divided by the length of the remainder of the skull. The hyphen indicates that no well-developed frill is present in the specimens. The asterisk denotes estimated length. All AMNH and FMNH specimens come from the Bayn Dzak locality, whereas the MDP and MPC specimens were collected at the Tugrikin Shire locality. Klingenberg, 2007) in the form of a significant slope. Following rec- (PIC) for both shape scores on log-transformed centroid sizes using ommendations in Klingenberg and Marugán-Lobón (2013), we also the tree shown in Figure 1. Slope and F test calculations were con- performed a regression of the phylogenetic independent contrasts ducted in the software PAST (Hammer, Harper,& Ryan, 2008). PRIETO-MÁRQUEZ et al. | 6295 We followed analyses outlined in Forth et al. (2016) in order to a punctuated model of evolution, to assess how sensitive resulting determine whether allometric shifts and possibly heterochrony oc- ontogenetic vectors are to branch lengths in the tree. curred in the evolutionary history of the ceratopsian frill. A pooled Allometric shifts were determined by subtracting adult regres- multivariate regression of Procrustes coordinates against log-trans- sion score species values from reconstructed ancestral ones with formed centroid size was conducted on the sample of adult skulls positive differences interpreted as peramorphosis, and negative (Tables 2 and 3), along with the smallest skull for each of the four species that were ontogenetically sampled. Ontogenetic vectors TA B L E 3 Sample of Triceratops horridus specimens used in this were reconstructed for those four species by connecting the small- study est and largest sampled skulls. Ancestral ontogenetic vectors were Basal skull Frill Relative subsequently reconstructed for Neoceratopsia, Coronosauria, and Specimen length (mm) inclination frill length Protoceratopsidae by optimizing adult and juvenile values for re- UCMP 154452 280 ? 0.27** gression scores and log centroid size on the time-scaled phylogeny BMRP cast of 496 31° 0.65 used for the pPCA analyses in the software Mesquite version 3.51 MOR 2951 (Maddison & Maddison, 2018) using square-changes parsimony. BMRP 854 32°** 0.78 Because our sample of ceratopsian skull shapes is heavily skewed 2006.4.1.2 toward Late Cretaceous Ceratopsidae with long ghost lineages YPM 1821 905 21° 0.91 among more basal neoceratopsians, implementation of the Brusatte YPM 1823 955 27° 0.72 et al. (2008) branch smoothing algorithm for adjusting zero-length FMNH PR 1,054 20° 0.79 branches leads to very deep branches subtending each species of 12300 protoceratopsid (~20 million years) and also Liaoceratops (>10 mil- MPM VP6841 1,110** ? ? lion years). Ongoing studies on basal neoceratopsian diversity (He AMNH 5116 1,147 38° 0.89 et al., 2015; Morshhauser, You, Li, & Dodson, 2019) and ages and Note: Measurements for MOR 1120, and UCMP 154452 from Horner correlations of Asian Late Cretaceous red beds (Dashzeveg et al., and Goodwin (2006), YPM 1823 from Hatcher et al. (1907), and YPM 2005; Fanti, Currie, & Badamgarav, 2012) suggest that turnover 1821 from Scannella and Horner (2010). The basal skull length is the among protoceratopsids may have occurred more rapidly within the distance between the tip of the rostrum and the caudal margin of the Campano-Maastrichtian (Makovicky, 2008). We therefore also gen- quadrate. Frill inclination was measured as the angle between the sagittal crest and the maxillary tooth row. Frill length is the maximum erated ancestral ontogenetic vectors from optimizations on a tree length of the frill relative to the length of the remainder of the skull. A in which branches without implied ghost lineage durations were single asterisk denotes those specimens used only in the SRVF analysis; corrected in Mesquite by addition of 1 million years, representing double asterisks indicate estimated values. F I G U R E 2 Results the morphometric analyses of the ceratopsian sample, excluding the parietal fenestra to allow inclusion of basal ceratopsians and Triceratops. (a) PCA bivariate plot summarizing the results of the Procrustes superimposition analysis for PC 1 and PC 2. (b) PCA bivariate plot for PC 1 and PC 3 of same 6296 | PRIETO-MÁRQUEZ et al. ones as paedomorphosis. Following Forth et al. (2016), we consid- the distribution of one hundred randomly generated angles to see ered allometric shifts between individual ancestor-descendant pairs whether they fall within a probability tail, or not. as significant if they were greater than 1.5 times the confidence Because the above method only samples the largest and smallest interval of the mean difference between all ancestor-descendant individual of each of the four ontogenetically sampled species with- pairs in our taxonomic sample. Differences between parallel vectors out accounting for the large individual variance in shape and size in represent either predisplacement if the descendant vector is shifted the full samples, we performed a second pooled multivariate regres- toward the origin of the biplot, or postdisplacement if the reverse sion based on the full range of specimens for just the four species holds. with growth series. We compared slopes for the ontogenetic vectors We estimated the rate of allometric change by dividing regres- using F tests, but did not consider regression score differences as sion score differences between adjacent nodes and/ or tips by the scores and log centroid size values were heavily biased by the differ- length of the branch separating them, and we also employed com- ences in size range and sample sizes between the individual species parative methods to examine whether the rate of frill evolution co- (which is why sampling was limited to two individuals in the previous varies with phylogeny, and whether it remained constant or varied analysis). Slope and F test calculations were conducted in the soft- over ceratopsian evolution. Comparative analyses were conducted ware PAST (Hammer et al., 2008). in the software BayesTraits Ver. 3 (http://www.evolu​tion.rdg.ac.uk) using the same tree topology and branch lengths as in the analyses for PIC and for the pPCA. We first calculated maximum likelihood 3 | R E S U LT S scores for the dataset under several different model settings to de- termine which model best applies to it in its entirety, and used AICc 3.1 | Interspecific frill variation in Ceratopsia and values to estimate the best-fitting overall model. Models evaluated Neoceratopsia were random walk, directional trend, and three models each favor- ing one of the three branch and/ or tree scaling parameters κ, δ, and Overall, variation in the ceratopsian frill (excluding frill ornamenta- λ. λ is a measure of whether trait evolution covaries with phylogeny tion and parietal fenestrae) is concentrated along its caudal and cau- with λ = 0 signifying evolution is independent of phylogeny and λ = 1 dolateral margins (Figure 2). For both Ceratopsia and Neoceratopsia, represents trait evolution covarying with phylogeny according to a the first three PC axes accounted for about 60% of the variation (with random walk model. κ is a measure of how the length of individual PC through PC10 encompassing over 95% of the variation), and PC1 branches scales to the rate of evolution. δ represents the scaling of encapsulated only 25% of the overall variation (Table 4). Basal cera- trait evolution to path length (patristic distance from root to tips), topsians showed a nonoverlapping (albeit nonsignificant, Table 5) with δ < 1 representing an “early burst” evolutionary scenario, δ > 1 distribution with neoceratopsians, appearing more segregated in often interpreted as accelerating evolution or “evolution via species morphospace (Figure 2a,b). Significant (Table 5) nonoverlapping adaptation,” and δ = 1 representing gradual evolution (Pagel, 1999). morphospace distributions were found between Centrosaurines and Individual analyses investigating the effects of branch scaling pa- both basal Ceratopsia and basal Neoceratopsia. Nonceratopsid ne- rameters allow for a finer scaled interpretation of the dynamics of frill oceratopsians and chasmosaurine ceratopsids showed overlapping shape evolution. Following Hernandez et al. (2013), we performed distributions when considering variation of along PC2 (Figure 2a), separate analyses for each of the scaling parameters to estimate the due to the relatively greater lateral expansion of the frill of these optimal value. We then ran constrained analyses for each parameter two groups. set to either 0 or 1 representing the various evolutionary scenarios de- fined above and compared the fit of these to the values to the uncon- TA B L E 4 Variance components for the first three PC axes for strained analyses using Bayes factors. In order to determine whether the Procrustes superimposition analyses for the two ceratopsian evolutionary rates vary across the phylogeny, the Variable Rates proto- groups considered col in BayesTraits ver. 3 was performed and marginal likelihoods were PC axis Ceratopsia Neoceratopsia estimated for both models with rate variation and without. Allometric differences related to shifts in ontogenetic vector an- PC1 25.4 24.4 gles between ancestor and descendant represent allometric accelera- PC2 19.1 18.3 tion or deceleration. Shifts in slope between ancestral and descendant PC3 17.2 16.5 ontogenetic vectors were measured as differences in their angles PC4 9.8 9.4 (Collyer & Adams, 2007) calculated as the arccosine of the dot product PC5 6.7 8.3 of the two vectors divided by the product of their lengths. We devised PC6 5.5 5.6 a simple randomization test to evaluate the significance of the angles. PC7 3.9 4.4 Using the largest and smallest components of the observed vectors as PC8 3.5 3.4 bounds, we generated 200 random vectors and measured the angles PC9 2.9 2.8 between successive pairs of random vectors using the same formula as PC10 1.7 2.5 for the observed vectors. The observed angles were compared against PRIETO-MÁRQUEZ et al. | 6297 TA B L E 5 Results of the pairwise p- permutational ANOVAs computed pPC axis Group 1 Group 2 F-value df1 df2 Value between ceratopsian groups on each of the PC axes PC1 Centrosaurinae Chasmosaurinae 29.15 1 14 0 PC1 Centrosaurinae Basal Neoceratopsia 10.70 1 8 .007 PC1 Centrosaurinae Basal Ceratopsia 149.64 1 8 .008 PC1 Chasmosaurinae Basal Neoceratopsia 0.02 1 10 .926 PC1 Chasmosaurinae Basal Ceratopsia 14.68 1 10 .004 PC1 Basal Neoceratopsia Basal Ceratopsia 5.87 1 4 .142 PC2 Centrosaurinae Chasmosaurinae 2.01 1 14 .186 PC2 Centrosaurinae Basal Neoceratopsia 1.85 1 8 .219 PC2 Centrosaurinae Basal Ceratopsia 1.01 1 8 .351 PC2 Chasmosaurinae Basal Neoceratopsia 0.09 1 10 .752 PC2 Chasmosaurinae Basal Ceratopsia 0.01 1 10 .926 PC2 Basal Neoceratopsia Basal Ceratopsia 0.18 1 4 .617 PC3 Centrosaurinae Chasmosaurinae 4.99 1 14 .053 PC3 Centrosaurinae Basal Neoceratopsia 9.66 1 8 .021 PC3 Centrosaurinae Basal Ceratopsia 72.78 1 8 .008 PC3 Chasmosaurinae Basal Neoceratopsia 9.68 1 10 .012 PC3 Chasmosaurinae Basal Ceratopsia 27.38 1 10 .001 PC3 Basal Neoceratopsia Basal Ceratopsia 5.34 1 4 .104 Note: p-Values are estimated by means of 1,000 randomizations. In this case, the parietal fenestra of the frill was excluded from the morphometric analysis. Abbreviation: df, degrees of freedom. F I G U R E 3 Results of the Procrustes superimposition analysis of neoceratopsians possessing a parietal fenestra. (a) PCA bivariate plot summarizing the results of the analysis for PC 1 and PC 2. (b) PCA bivariate plot for PC1 and PC3 There is a general trend of progressive caudal elongation and caudomedial extension of the caudal border (PC2, Figure 2a) and rostrolateral expansion of the caudal and lateral margins of the caudolateral projection of the caudolateral corner (PC3, Figure 2b) frill, respectively, from basal ceratopsians to centrosaurine cera- of the frill. Basal ceratopsians occupied the region of the morpho- topsids (Figure 2a,b). Other components of frill variation include space defined by a rostrocaudally abbreviated caudal region of the 6298 | PRIETO-MÁRQUEZ et al. TA B L E 6 Results of the pairwise p- permutational ANOVAs computed pPC axis Group 1 Group 2 F-value df1 df2 Value between ceratopsian groups on each of PC1 Centrosaurinae Chasmosaurinae 21.79 1 13 0 the PC axes PC1 Centrosaurinae Basal Neoceratopsia 8.83 1 8 .005 PC1 Chasmosaurinae Basal Neoceratopsia 0.05 1 9 .816 PC2 Centrosaurinae Chasmosaurinae 9.32 1 13 .004 PC2 Centrosaurinae Basal Neoceratopsia 2.06 1 8 .206 PC2 Chasmosaurinae Basal Neoceratopsia 0.98 1 9 .359 PC3 Centrosaurinae Chasmosaurinae 4.66 1 13 .05 PC3 Centrosaurinae Basal Neoceratopsia 6.69 1 8 .047 PC3 Chasmosaurinae Basal Neoceratopsia 7.35 1 9 .015 Note: p-Values are estimated by means of 1,000 randomizations. In this case, the parietal fenestra of the frill was included in the morphometric analysis. Abbreviation: df, degrees of freedom. skull, whereas centrosaurine ceratopsids are at the other end of the Evidence for frill modularity was weaker, though still statistically spectrum featuring a rostrocaudally elongate frill. significant, in the ontogenetic series of Liaoceratops yanzigouen- Morphospace occupation patterns in fenestrated neoceratop- sis, the Ukhaa Tolgod protoceratopsid, Protoceratops andrewsi and sians broadly mirror those found for the dataset that excluded the Triceratops horridus (Figure 6d) as CR coefficients for all taxa ap- fenestra, despite the exclusion of several taxa with solid frills such proached unity (Figure 6). Furthermore, those CR values were lo- as basal ceratopsians and the chasmosaurine ceratopsid Triceratops. cated further from the tails of the distributions (Figure 6) than in Again, substantial overlap in morphospace occupation between the clade analyses reported above (Figures 4 and 5), indicating the chasmosaurines and basal coronosaurs is observed in the PC1 ver- existence of a number of possible landmark partitions in these four sus PC2 plot (Figure 3a), due to the greater lateral expansion of the species with a greater potential for within-module covariation than frill in these two groups. However, when considering PC1 versus in the frill partition. PC3, basal neoceratopsians, chasmosaurines, and centrosaurines The iterative subsampling analyses with removal of 10, 20, and occupy significantly distinct (Table 6) nonoverlapping regions of the 28 semilandmarks resulted in very similar CR values for these cera- morphospace (Figure 3b). Both the frill and parietal fenestra show a topsian groups as were recovered for the full dataset, although the strong trend toward rostrocaudal lengthening from basal neocera- spread of their distributions slightly increased progressively and CR topsians to ceratopsids, particularly in centrosaurines, as captured values migrated slightly away from the tail end of the such distribu- by both PC1 and PC2 (Figure 3a). tions (Appendices S8–S10). The latter effect was more noticeable after removal of 20 and 28 semilandmarks (Appendices S9 and S10). However, results indicate that our overall results were not biased by 3.2 | Tests of frill modularity inclusion of a large number of semilandmarks. The recovered CR coefficients provided varying degrees of support for modularity of the frill in the four species and within different 3.3 | Heterochrony of the frill groupings of ceratopsians, with actual CR values consistently located at the tail end of the distribution of randomized values (Figures 4–6). We found a significant (p-value = .0001) positive (slope = 0.196) al- Weak but significant support for modularity was found with CR val- lometric relationship between regression scores and log centroid ues approaching one in the grade of basal ceratopsians (Figure 4a), size for adults across Ceratopsia, and this relationship was robust Ceratopsia as a whole (Figure 4b) and in basal neoceratopsians to phylogenetic correction using PIC values (p-value = .0053; (Figure 4c). Support for the modularity of the frill was stronger in the slope = 0.112; Figure 7). Reconstructed ontogenetic vectors for clades Neoceratopsia (Figure 4d,e) and Ceratopsidae (Figures 4f and Liaoceratops and the last common ancestor (LCA) of Neoceratopsia 5a), as evidenced by the substantially lower CR and significance val- are flat to slightly negative, but slopes are positive in the LCAs of ues. Within Ceratopsidae, consideration of each subclade separately Coronosauria and Protoceratopsidae, and in the individual species resulted in weaker support for frill modularity (Figure 5b–d) likely they subtend (Figure 8). Positive differences in regression scores due to sample size. However, inclusion of the parietal fenestra in the greater than the significance threshold (Appendix S11) are observed analysis resulted in higher support for modularity in Chasmosaurinae between the LCAs of Ceratopsia and Neoceratopsia, between the (lower CR value located nearer the end of the distribution) relative to LCAs of Neoceratopsia and Coronosauria, and between the latter Centrosaurinae. and the LCA of Protoceratopsidae. Significant positive differences PRIETO-MÁRQUEZ et al. | 6299 6300 | PRIETO-MÁRQUEZ et al. F I G U R E 4 Distribution of covariance ratios (CR) for all the evaluated 10,000 random subsets of landmarks for various paraphyletic and monophyletic ceratopsian groups. The red arrow indicates the CR for the hypothesis of modularity of the frill. The frame diagrams display the landmarks digitized on the dorsal view of the skull; red landmarks correspond to those of the hypothesized frill module are also observed in the branches immediately subtending the and Chasmosaurinae with significant regression score differences lie LCA of Ceratopsidae and in the terminal branches for Liaoceratops, within the two lowest rate estimate quartiles, suggesting that het- Einiosaurus, and Chasmosaurus. Significant negative differences erochronic shifts were achieved by incremental changes over long are calculated for the LCA of Chasmosaurinae and the LCA of periods. By contrast, rapid changes in frill shape are principally ob- (Regaliceratops, Triceratops), in each of those two individual species, served in individual ceratopsid species and shallow subclades, espe- and also in Protoceratops, Vagaceratops, and Pachyrhinosaurus. cially among Chasmosaurinae, which dominate the top quartile of Dividing the Regression Score differences by branch lengths estimated rate values (Appendix S11). provides an estimated rate of heterochronic change (Appendix S11). Likelihood estimation for the five different evolutionary mod- When comparing these rates, it is notable that rates for LCAs of sev- els overwhelmingly favors random walk (= Brownian motion) as eral clades including Coronosauria, Protoceratopsidae, Ceratopsidae, the best-fitting model for frill shape evolution across ceratopsian F I G U R E 5 Distribution of covariance ratios (CR) for all the evaluated 10,000 random subsets of landmarks for various paraphyletic and monophyletic ceratopsian groups. The red arrow indicates the CR for the hypothesis of modularity of the frill. The frame diagrams display the landmarks digitized on the dorsal view of the skull; red landmarks correspond to those of the hypothesized frill module PRIETO-MÁRQUEZ et al. | 6301 F I G U R E 6 Distribution of covariance ratios (CR) for all the evaluated 10,000 random subsets of landmarks for growth series of basal neoceratopsians Liaoceratops, a new protoceratopsid taxon, Protoceratops and the chasmosaurine ceratopsid Triceratops. The red arrow indicates the CR for the hypothesis of modularity of the frill. The frame diagrams display the landmarks digitized on the dorsal view of the skull; red landmarks correspond to those of the hypothesized frill module diversity. All single parameter analyses therefore employed the supported, is not significantly different from either K = 0 or K = 1 random walk model. Results of our scaling parameter analyses in in our dataset, so a gradual (i.e., scaled to branch length) model of BayesTraits are given in Appendix S12. In our unconstrained anal- trait evolution cannot be excluded. The best-supported model for ysis, the estimate for L = 0.88 but is not significantly different from the path length scaling parameter is the one that estimates D = 1.9, either L = 0 (trait evolution independent of phylogeny) or L = 1 (trait but it is not significantly different from 1 as determined by log Bayes evolution tracks phylogeny). Parsing this result at finer taxonomic factor differences, again not allowing us to exclude gradual evolu- levels, L = 0 in centrosaurine and chasmosaurine ceratopsids, but tion. However, testing for variable rates found positive evidence for is closer to one in nonceratopsid taxa indicating an abrupt shift in a rate shift in our data (log Bayes factor difference > 3). Over 28% this parameter during ceratopsian frill evolution. The branch scal- percent of Bayesian Monte Carlo Markov Chain iterations posite a ing parameter was estimated to be 0.465 suggesting that shorter rate shift within Chasmosaurinae, with rates above 65% in the chas- branches (i.e., those among ceratopsids) capture more evolution- mosaurine subclade that includes Chasmosaurus, Vagaceratops, and ary change than long branches. However, this model, while better Mojoceratops. Percentages of iterations within centrosaurines are 6302 | PRIETO-MÁRQUEZ et al. F I G U R E 7 Regression analysis of ontogenetic trajectories from four neoceratopsian species against log-transformed centroid size using the regression score as shape variable. (a) Simplified phylogeny of Ceratopsia with the species included in the regression analyses, based on the trees in Xu et al. (2006), Farke et al. (2014), Longrich (2014), Brown and Henderson (2015), and Evans and Ryan (2015); labeled nodes correspond to those appearing in b. (b) Regression analysis using phylogenetic contrasts. (c) Nonphylogenetically corrected regression analysis. Species abbreviations are as in Table 1 mostly lower (18% < 32%) than among chasmosaurines, but are still Triceratops, even though their angles do not differ substantially. This higher than among any nodes and branches that subtend noncera- shift toward larger values in later diverging (i.e., descendant) taxa is topsid taxa, mirroring the distribution for rate estimates reported in consistent with hypermorphosis according to Alberch et al. (1979) Appendix S11 and supporting one or more shifts in the rate of frill (Figure 8). However, because of the great difference in sampling of evolution in ceratopsids. the youngest specimens between species (neonates in Tengr and Angles between the ontogenetic vectors of Neoceratopsia Protoceratops versus older juveniles in Liaoceratops and Triceratops), and Coronosauria are more than three times greater than the an- we cannot assess whether or postdisplacement may also have con- gles between any other ancestor-descendant pair of vectors in the tributed to these results. analysis regardless of branch length adjustment method indicating heterochrony through an allometric shift in growth rate at this node (Figure 8). The observed angles are greater than 85% of angles gen- 4 | D I S CU S S I O N erated in our randomization tests, suggesting that these differences are substantial, though not necessarily significant at the 5% level. All 4.1 | Shape variation across Ceratopsia other angles between ontogenetic vectors do not differ noticeably from random. Frill shape distributions in morphospace exhibit a phylogenetic Within Neoceratopsia, a progression toward higher centroid pattern even after correction for phylogenetic covariance, with size and regression score is observed in the adult endpoint in the early-diverging forms occupying their own region of morphospace vectors for the LCAs of Coronosauria, Protoceratopsidae, and in distinct from neoceratopsians in all analyses they were included in PRIETO-MÁRQUEZ et al. | 6303 F I G U R E 8 Regression analysis of ontogenetic trajectories from four neoceratopsian species against log-transformed centroid size using the regression score as shape variable. Solid lines represent observed taxa. Ontogenetic vectors for reconstructed ancestors are dashed to distinguish them from species. Short dashed vectors with circular points indicate reconstructions based on trees for which branch lengths were smoothed using the Brusatte et al. (2008) method. Long dashed vectors with triangular data points show reconstructions based on trees calibrated against record with no branch smoothing, but in which 1 Ma is added to taxa with zero-branch length ghost lineages (i.e., oldest sister taxon in a clade). Hollow stars and circles indicate chasmosaurine and centrosaurine ceratopsid species, respectively, and pink ellipses indicate basal ceratopsian species; the abbreviations for all those species are as in Table 1 (Figures 2 and 3). Maiorino, Farke, Kotsakis, and Piras (2017) also & Horner, 2011; Ryan et al., 2007; Sampson & Loewen, 2010). recovered a phylogenetic structuring in the cranial morphospace of These qualitative observations are largely supported by our quan- most ceratopsian clades, although their morphometric analysis con- titative results with ceratopsids occupying unique and markedly sidered the ceratopsian skull in lateral view. We observed further larger part of morphospace than nonceratopsids. Notably, the phylogenetic structuring in morphospace within Neoceratopsia as expansion of frill morphospace in ceratopsids coincides with in- well, with chasmosaurines overlapping in morphospace occupation creased rates of shape change (Appendices S11 and S12), which (PC1/PC2) with early-diverging neoceratopsians (Figures 2a and 3a), can be interpreted as a conservative proxy for frill evolutionary while centrosaurines occupy their own sector to the left of all other rates since the elaboration of episquamosal and epiparietal orna- species (Figures 2a,b and 3). A greater separation between groups is mentation is not accounted for using our methodology. These re- observed in the PC1/PC3 morphospace, with early-diverging taxa sults are in agreement with the accelerated evolutionary rates in occupying the upper right quadrant and progressively later diverging cranial morphology found by Maiorino et al. (2017) in ceratopsids taxa moving to the left along PC3 (caudally wider frills; Figures 2b relative to early-diverging clades. While in our study this pattern and 3b). is undoubtedly related to the uneven sampling across ceratopsian Ceratopsian frills, and in particular those of ceratopsids, have clades and grades with a bias toward ceratopsids, qualitative com- been interpreted as highly variable by multiple authors (e.g., Brown parisons to frill shape and anatomy among unsampled basal cer- & Henderson, 2015; Dodson & Currie, 1990; Dodson et al., 2004; atopsians (Makovicky & Norell, 2006; Morshhauser et al., 2019; Evans & Ryan, 2015; Knell & Sampson, 2011; Ostrom, 1966; Padian Xu, Makovicky, Wang, Norell, & You, 2002) indicates that they are 6304 | PRIETO-MÁRQUEZ et al. generally similar to sampled noncoronosaurian taxa in being rel- acting on the structure under study do so in a modular fashion atively short, narrow, and having a round profile with a gradual (Klingenberg, 2014). rather than stepped parietal-squamosal contact in dorsal aspect. We therefore see no indication that they would markedly increase morphospace occupation of grades stemward of Coronosauria. 4.3 | Heterochrony As one of the first quantitative studies of ceratopsian cranial shape, 4.2 | Modularity we found substantial support for heterochronic processes in the evolution of the hallmark parieto-squamosal frill that diagnoses the We found a strong interspecific signal for modularity (Figures 4 and clade. Acceleration is the main mechanism that accounts for the 5), but a weak (though still significant) signal at intraspecific levels macroevolutionary shift in frill shape and size between the LCA of within ontogeny (Figure 6). Recently, Maiorino et al. (2017) found Neoceratopsia with a short, narrow frill throughout ontogeny, and in ceratopsians a greater degree of integration between the frill and coronosaurs, which have significantly longer and wider frills in adults the mandible than between the frill and rest of the skull, which is (Hone et al., 2016). Furthermore, shifts in Regression Scores corre- consistent with the interspecific modularity of the frill hypothe- lated with shifts in centroid size between early-diverging corono- sized in our study. The weak modularity signal at intraspecific levels saurs and the Triceratops ontogenetic vector are broadly consistent within ontogeny may relate to smaller samples sizes for individual with a peramorphocline between early-diverging coronosaurs and species compared to clades, and also the greater effects of tapho- ceratopsids as suggested by Long and McNamara (1995), although nomic distortion on the more limited morphospace of intraspecific the exact details differ. samples. Other morphometric analyses of ceratopsian species in- Previous discussions of heterochrony in ceratopsians have fo- cluding Psittacosaurus lujiatunesis (Hedrick & Dodson, 2013) and cused heavily on ceratopsids. Because of the lack of growth series Protoceratops andrewsi (Maiorino, Farke, Kotsakis, & Piras, 2015) also for almost all taxa, our analysis cannot quantitatively evaluate these recovered a large amount of taphonomically induced shape variation previous qualitative generalizations (e.g., Long & McNamara, 1995; that overruled earlier studies suggesting greater taxonomic diversity McNamara & Long, 2012; Tumarkin & Dodson, 1998), but we note (You & Dodson, 2004) or sexual dimorphism (Dodson, 1976) for the that high rates of frill shape change accompany both negative and two taxa respectively. Despite the taphonomic effects at species positive allometric shifts between adults of various species and level, our results are congruent with the well-established hypoth- their respective LCAs, indicating that a combination of both pae- esis positing that modularity is requisite for heterochronic evolution domorphic and peramorphic processes underpin the great disparity (Gerber & Hopkins, 2011; Olson & Rosell, 2006; Raff & Raff, 2000). of ceratopsian frill shapes. Such mosaic heterochrony may facilitate Studies on other vertebrates have shown that ontogenetic the evolution of disparate morphologies (Guerrant, 1982), which is patterns of variation in modularity and integration are variable, consistent with the morphospace occupation patterns we observe in depending on the skull components and characters (Goswami, ceratopsids (Figures 2 and 7a). Polly, Mock, & Sánchez-Villagra, 2012; Klenovsek & Jojic, 2016; Although we provide the first strong quantitative evidence for Zelditch & Carmichael, 1989) and taxa (Goswami & Polly, 2010; heterochrony in ceratopsian skull evolution, heterochrony alone Willmore, Leamy, & Hallgrimsson, 2006; Zelditch, Bookstein, & is most likely insufficient as a mechanism to explain frill disparity. Lundrigan, 1992) under consideration. The decrease in frill modu- Innovation was probably an important process in ceratopsian diver- larity within ontogenetic series of select species relative to that we sification, especially through the evolution of the episquamosal and observe at macroevolutionary levels in ceratopsians doubtless re- epiparietal structures that are characteristic of Ceratopsidae. These flects taphonomic distortion as noted above, with the larger shape structures are, however, difficult to homologize for the purposes of differences between exemplars of different species less affected identifying landmarks suitable for geometric morphometric analy- by such artifacts, than are the smaller differences between speci- sis, and the quantitative analysis of their morphology will require mens of a single species. However, we cannot exclude that some of other methods such as outline-based techniques (e.g., Joshi, Prieto- this difference in the robustness of our CR results above and below Márquez, & Parker, 2011) to capture the true disparity of frill forms. the species level could be explained by the “palimpsest” metaphor Overall patterns of morphospace occupation (Figure 2) support (Hallgrímsson et al., 2009). According to this concept, observed pat- the results of our analyses of evolutionary rates (Appendices S11 terns of modularity or integration in a particular specimen of a devel- and S12), which reveal a marked shift in the tempo and mode of opmental series are the cumulative effects of all previous processes frill evolution. Morphospace occupation and evolutionary rates are that occurred prior to that stage (Klingenberg, 2014). Processes tak- constrained in early-diverging nonceratopsid taxa, with frill evo- ing place in later ontogenetic stages may only partially “overwrite” lution happening gradually and largely predicted by phylogeny. the patterns that occurred in previous stages. In this way, the overall Ceratopsid frill evolution, on the other hand, is characterized by result of overlaying modular patterns on each other throughout de- significantly greater disparity driven by evolution at elevated rates velopment may lead to a higher degree of integration where there that is independent of phylogeny. Elucidating the biological drivers is less statistical support for modularity, even though the processes behind this change is challenging, as they are likely multiple and PRIETO-MÁRQUEZ et al. | 6305 complex, but the ceratopsian literature (e.g., Makovicky, 2012; results (Figures 2 and 3) where ceratopsids occupy discrete morpho- Makovicky & Norell, 2006) suggests that the observed accelera- space from basal ceratopsians characterized by more expanded and tion between the neoceratopsian LCA and coronosaurs, as well as variable frills regardless of whether or not fenestrae are excluded. subsequent high rates of evolution in ceratopsids, may relate to a Finally, if the shift in frill function entails a functional release of the shift in frill function. caudal frill border from serving as anchorage for jaw adductor mus- Although myriad hypotheses have been put forth for the function of cles (Dodson, 1996; Makovicky, 2012) and frees it up to evolve in re- the ceratopsian frill, recent studies favor either feeding (e.g., Maiorino sponse to other selective pressures, we may predict that frills exhibit et al., 2017) or intraspecific signaling (e.g., Horner & Goodwin, 2006) greater disparity and an increased rate of evolution. Both of these as the primary function for this structure. The frills of basal neocer- predictions are also met by our results, with our comparative anal- atopsians exhibit a short, deep morphology with muscle scars along yses revealing that the increased disparity of ceratopsid frills was most of the rim providing evidence of the insertion of enlarged jaw achieved by an increase in evolutionary rates and a decoupling of adductor muscles (Makovicky & Norell, 2006; Morshhauser et al., trait evolution from phylogeny. It is worth stressing that these pat- 2019; Xu et al., 2002). Muscle scarring is much less developed in the terns are robust even without considering the highly disparate epi- expanded frills of coronosaurs. Current reconstructions of jaw muscles parietal and episquamosal ossifications that characterize ceratopsid in these animals place the jaw adductor origins closer to the rostral frill ornamentation. edge of the supratemporal fenestra, along the midline keel of the frill, and along the inner surface of the squamosal, with the majority of the frill, including its caudal margin, free of musculature (Dodson, 1996; 5 | CO N C LU S I O N S Nabavizadeh, 2018). These differences in frill anatomy between early and late diverging ceratopsians prompted Makovicky (2012) The parieto-squamosal frill of ceratopsian dinosaurs is a neomorphic and Makovicky and Norell (2006) to propose that the ceratopsian frill structure that has been proposed to have evolved via heterochrony. underwent a shift in its primary function within Neoceratopsia from Through geometric morphometric methods, we demonstrate that providing an increased area for attachment of jaw adductor muscles frill shapes vary by clade and also that the highly variable frills of as (Haas, 1955; Lull, 1908; Ostrom, 1964; Xu et al., 2002), to serving as later diverging ceratopsids occupy a greater area of morphospace a species-specific display structure with possible social and/or sexual relative to the more constrained occupation of earlier diverging functions (Dodson, 1976, 1993, 1996; Farlow & Dodson, 1975; Hone, taxa. Our morphometric analyses demonstrate that the principal Naish, & Cuthill, 2012; Hone et al., 2016; Horner & Goodwin, 2006, mechanism of macroevolutionary frill shape is expansion along the 2008; O'Brien et al., 2018; Padian & Horner, 2011; Sampson, Ryan, & caudolateral “corners” and caudal margin of the frill. In the more de- Tanke, 1997). Insertion of jaw musculature along the majority of the rived ornamented ceratopsians that radiated during the Campanian frill margin in early-diverging taxa such as Liaoceratops, Yamaceratops and Maastrichtian (latest Cretaceous), morphological variation is (Makovicky & Norell, 2006), and Auroraceratops (Morshhauser et al., concentrated along the caudolateral (Chasmosaurinae) and lateral 2019) would likely have constrained allometric growth patterns by (Centrosaurinae) margins of the frill. optimizing muscle orientations and mechanical leverage, a conclusion The frill constituted a module at both macroevolutionary, and to supported by the observed absence of ontogenetic shape change in a lesser extent, ontogenetic levels, and evolved via peramorphosis Liaoceratops and the LCA of Neoceratopsia (Figure 7a). However, within Ceratopsia. This included acceleration of frill shape change both anatomical data (Dodson, 1993; Farke, Chapman, & Andersen, early in ceratopsian history between early-diverging neoceratop- ; Nabavizadeh, 2019; Rieppel, 1981) and biomechanical models (Bell, sians and coronosaurs, followed by peramorphosis, possibly by Snively, & Shychoski, 2009) support the conclusion that most of the hypermorphosis, between early-diverging coronosaurs and the cer- frill margin was free of muscular attachment in ceratopsids, and likely atopsid radiation. Within the latter radiation, fast rates of frill shape all coronosaurs. Thus released from ancestral constraints, the frill was evolution decoupled from phylogeny produced both peramorphic able to respond more rapidly to other selective pressures such as so- and paedomorphic shape changes, indicating that mosaic heteroch- ciosexual signaling (Hone et al., 2012; Horner & Goodwin, 2006, 2008; rony accounts for the great disparity in frill shapes. The constrained Sampson et al., 1997). morphospace occupation and low rates of frill shape evolution While our morphometric and comparative analyses cannot di- among early-diverging ceratopsians when compared to ceratopsids rectly test such functional hypotheses, they can be used to evaluate are consistent with a shift in frill function from primarily serving as predictions arising from them. For the frill to shift function in cera- a platform for enlarged jaw adductor muscles, to being a highly vari- topsian evolution, it has to first represent a discrete functional unit, able display structure. This hypothesis along with better insights on an inference that is consistent with our results for modularity across the heterochronic processes within ceratopsids await a denser sam- multiple clades and grades of ceratopsian evolution and also within pling within ceratopsian growth series and across ceratopsian diver- individual species. Secondly, we predict that a shift in the primary sity. Nevertheless, our study provides the first quantitative evidence function of a modular trait should be mirrored in shifts in morpho- supporting previous assertions that peramorphosis played a key role space occupation. This prediction is also met by our morphospace in the early evolution of the ceratopsian frill. 6306 | PRIETO-MÁRQUEZ et al. AC K N OW L E D G M E N T S REFERENCES We are grateful to Leonardo Maiorino, Andrew Farke and Matt Adams, D. C. (2016). Evaluating modularity in morphometric data: Lamana for providing digital images of Protoceratops and to Joseph Challenges with the RV coefficient and a new test measure. Methods in Ecology and Evolution, 7(5), 565–572. A. Frederickson and Michael K. Brett-Surman for sharing pictures Adams, D. C., & Otárola-Castillo, E. (2013). geomorph: An R package for of Centrosaurus specimens. Carl Mehling, John Scannela, Patricia C. the collection and analysis of geometric morphometric data. Methods Burke, Peter Sheehan, Joshua Mathews, and Dallas Evans provided in Ecology and Evolution, 4, 393–399. access to ceratopsian specimens under their care. Thanks also to Alberch, P., Gould, S. J., Oster, G. F., & Wake, D. B. (1979). Size and shape in ontogeny and phylogeny. Paleobiology, 5, 296–317. Khishigjaw Tsogtbaatar for sharing locality data on MPC specimens Bell, M. A., & Lloyd, G. T. (2015). Strap: An R package for plotting phylog- of P. andrewsi. Christian P. Klingenberg and Dean Adams provided enies against stratigraphy and assessing their stratigraphic congru- advice on the use of MorphoJ and computation of the CR coeffi- ence. Palaeontology, 58, 379–389. cient, respectively, while Andrew Meade kindly provided advice on Bell, P. H., Snively, E., & Shychoski, L. (2009). A comparison of the jaw me- chanics in hadrosaurid and ceratopsid dinosaurs using Finite Element BayesTraits. Ali Navabizadeh provided useful insights on ceratopsian Analysis. Anatomical Record, 292(9), 1338–1351. temporal and jaw musculature. This project was supported by a Bass Bhullar, B.-A., Marugán-Lobón, J., Racimo, F., Bever, G. S., Rowe, T. B., Postdoctoral Fellowship from the Field Museum of Natural History Norell, M. A., & Abzhanov, A. (2012). Birds have paedomorphic di- and the Ramón y Cajal program from the Ministry of Economy, nosaur skulls. Nature, 487(7496), 223–226. https://doi.org/10.1038/ Industry and Competitivity of Spain (RyC-2015-17388) presented to natur​e11146 Bookstein, F. L. (1991). Morphometric tools for landmark data: Geometry APM. Data collection was also supported by NSF EAR 0418648 to and biology. New York, NY: Cambridge University Press. PJM. PJM is currently supported by NSF awards PLR 1341645 and Brown, B., & Schlaikjer, E. M. (1940). The structure and relationships of FRES 1925884. JGP is supported by a Juan de la Cierva Fellowship Protoceratops. Annals of the New York Academy of Science, 40(3), 133– from the Ministry of Economy, Industry and Competitivity of Spain 266. https://doi.org/10.1111/j.1749-6632.1940.tb570​47.x Brown, C. M., & Henderson, D. M. (2015). A new horned dinosaur reveals (FJCI-2014-20380). Additional support was also provided by the convergent evolution in cranial ornamentation in Ceratopsidae. Current CERCA Programme of the Generalitat de Catalunya. Biology, 25(12), 1–8. https://doi.org/10.1016/j.cub.2015.04.041 Brusatte, S. L., Benton, M. K., Ruta, M., & Lloyd, G. T. (2008). Superiority, C O N FL I C T O F I N T E R E S T S competition, and opportunism in the evolutionary radiation of dino- saurs. Science, 321(5895), 1485–1488. https://doi.org/10.1126/scien​ None declared. ce.1161833 Collyer, M. L., & Adams, D. C. (2007). Analysis of two-state multivariate AU T H O R C O N T R I B U T I O N phenotypic change in ecological studies. Ecology, 88, 683–692. Albert Prieto-Márquez: Conceptualization (equal); Data cura- Dashzeveg, D., Dingus, L., Loope, D. B., Swisher, C. C. III, Dulam, T., & Sweeney, M. R. (2005). New stratigraphic subdivision, depositional tion (equal); Formal analysis (equal); Funding acquisition (equal); environment, and age estimate for the Upper Cretaceous Djadokhta Investigation (equal); Methodology (equal); Project administration Formation, southern Ulan Nur Basin, Mongolia. American Museum (equal); Resources (equal); Software (equal); Supervision (equal); Novitates, 3498, 1–31. https://doi.org/10.1206/0003-0082(2005)4 Validation (equal); Visualization (equal); Writing-original draft (equal); 98[0001:NSSDE​A]2.0.CO;2 Dodson, P. (1976). Quantitative aspects of relative growth and sexual di- Writing-review & editing (equal). Joan Garcia-Porta: Formal analysis morphism in Protoceratops. Journal of Paleontology, 50(5), 929–940. (equal); Investigation (equal); Methodology (equal); Software (equal); Drake, A. G., & Klingenberg, C. P. (2007). The pace of morphological Supervision (equal); Validation (equal); Visualization (equal); Writing- change: Historical transformation of skull shape in St Bernard dogs. review & editing (equal). Shantanu H. Joshi: Formal analysis (equal); Proceedings of the Royal Society B, 275, 71–76. Investigation (equal); Methodology (equal); Software (equal). Mark Dodson, P. (1993). Comparative craniology of the Ceratopsia. American Journal for Science, 293-A, 200–234. https://doi.org/10.2475/ajs.293.A.200 A. Norell: Data curation (equal); Resources (equal); Supervision Dodson, P. (1996). The Horned Dinosaurs: A natural history. Princeton, NJ: (equal); Validation (equal); Visualization (equal); Writing-review Princeton University Press. & editing (equal). Peter J. Makovicky: Conceptualization (equal); Dodson, P. (2013). Ceratopsia increase: History and trends. Canadian Data curation (equal); Formal analysis (equal); Funding acquisi- Journal of Earth Sciences, 50, 294–305. Dodson, P., & Currie, P. J. (1990). Neoceratopsia. In D. B. Weishampel, P. tion (equal); Investigation (equal); Methodology (equal); Resources Dodson, & H. Osmólska (Eds.), The Dinosauria (2nd ed., pp. 593–618). (equal); Software (equal); Supervision (equal); Validation (equal); Berkeley, CA: University of California Press. Visualization (equal); Writing-original draft (equal); Writing-review & Dodson, P., Forster, C. A., & Sampson, S. D. (2004). Ceratopsidae. In D. editing (equal). B. Weishampel, P. Dodson, & H. Osmólska (Eds.), The Dinosauria (2nd ed., pp. 494–513). Berkeley, CA: University of California Press. Erickson, G. M., & Druckenmiller, P. S. (2011). Longevity and growth DATA AVA I L A B I L I T Y S TAT E M E N T rate estimates for a polar dinosaur: A Pachyrhinosaurus (Dinosauria: The data used in this research that is not included in the tables are Neoceratopsia) specimen from the North Slope of Alaska showing a provided in the appendices listed above. complete developmental record. Historical Biology, 23(4), 327–334. Erickson, G. M., Makovicky, P. J., Currie, P. J., Norell, M. A., Yerby, S. A., & Brochu, C. A. (2004). Gigantism and comparative life-history pa- ORCID rameters of tyrannosaurid dinosaurs. Nature, 430(7001), 772–775. Albert Prieto-Márquez https://orcid.org/0000-0002-4836-4052 https://doi.org/10.1038/natur​e 02699 PRIETO-MÁRQUEZ et al. | 6307 Erickson, G. M., Makovicky, P. J., Inouye, P. D., Keqin, G., & Zhou, C.- Gould, S. J. (1974). The origin and function of “bizarre” structures: Antler F. (2009). A Life Table for Psittacosaurus lujiatunensis: Initial insights size and skull size in the “Irish elk”, Megaloceros giganteus. Evolution, into ornithischian dinosaur population biology. Anatomical Record- 28(2), 191–220. Advances in Integrative Anatomy and Evolutionary Biology, 292(9), Gould, S. J. (1977). Ontogeny and phylogeny. Cambridge, UK: Harvard 1514–1521. University Press. Evans, D. C., & Ryan, M. J. (2015). Cranial anatomy of Wendiceratops pin- Gould, S. J. (1992). Heterochrony. In E. F. Keller, & E. A. Lloyd (Eds.), hornensis gen. et sp. Nov., a centrosaurine ceratopsid (Dinosauria: Keywords in evolutionary biology (pp. 158–165). Cambridge, UK: Ornithischia) from the Oldman Formation (Campanian), Alberta, Harvard University Press. Canada, and the evolution of ceratopsid nasal ornamentation. . PLoS Guenther, M. F. (2009). Influence of sequence heterochrony on hadro- ONE, 10(7), e0130007. https://doi.org/10.1371/journ​al.pone.0130007 saurid dinosaur postcranial development. Anatomical Record, 292(9), Fanti, F., Currie, P. J., & Badamgarav, D. (2012). New specimens of 1427–1441. Nemegtomaia from the Barrungoyot and Nemegt Formations (Late Guerrant, E. O., Jr. (1982). Neotenic evolution of Dephinium nudicaule Cretaceous) of Mongolia. PLoS ONE, 7(2), e31330. https://doi. (Ranunculaceae): A hummingbird-pollinated larkspur. Evolution, org/10.1371/journ​al.pone.0031330 36(4), 699–712. Farke, A. A., Chok, D. J., Herrero, A., Scolieri, B., & Werning, S. Haas, G. (1955). The jaw musculature in Protoceratops and in other cera- (2013). Ontogeny in the tube-crested dinosaur Parasaurolophus topsians. American Museum Novitates, 1792, 1–24. (Hadrosauridae) and heterochrony in hadrosaurids. PeerJ, 1, e182. Hallgrímsson, B., Jamniczky, H., Young, N. M., Rolian, C., Parsons, T. E., https://doi.org/10.7717/peerj.182 Boughner, J. C., & Marcucio, R. S. (2009). Deciphering the palimp- Farke, A. A., Maxwell, W. D., Cifelli, R. L., & Wedel, M. J. (2014). A cer- sest: Studying the relationship between morphological integration atopsian dinosaur from the Lower Cretaceous of Western North and phenotypic covariation. Evolutionary Biology, 36(4), 355–376. America, and the biogeography of Neoceratopsia. PLoS ONE, 9(12), https://doi.org/10.1007/s1169​2-009-9076-5 e112055. https://doi.org/10.1371/journ​al.pone.0112055 Hammer, Ø., Harper, D. A. T., & Ryan, P. D. (2001). Past: Paleontological Farke, A. A., Ryan, M. J., Barrett, P. M., Tanke, D. H., Braman, D. R., statistics software package for education and data analysis. Loewen, M. A., & Graham, M. R. (2011). A new centrosaurine from Palaeontologia Electronica, 4(1), art. 4, 9 pp. http://palaeo-elect​ronica. the Late Cretaceous of Alberta, Canada, and the evolution of parietal org/2001_1/past/issue1_01.htm ornamentation in horned dinosaurs. Acta Palaeontologica Polonica, Hatcher, J. B., Marsh, O. C., & Lull, R. S. (1907). The Ceratopsia. U. S. 56(4), 691–703. Geological Survey Monograph, 49, 1–300. Farke, A. A., Wolf, E. D. S., & Tanke, D. H. (2009). Evidence of combat He, Y., Makovicky, P. J., Wang, K., Chen, S., Sullivan, C., Han, F., & Xu, in Triceratops. PLoS ONE, 4(1), e4252. https://doi.org/10.1371/journ​ X. (2015). A new leptoceratopsid (Ornithischia, Ceratopsia) with a al.pone.0004252 unique ischium from the Upper Cretaceous of Shandong Province, Farlow, J. O., & Dodson, P. (1975). The behavioral significance of frill and China. PLoS ONE, 10(12), e0144148. https://doi.org/10.1371/journ​ horn morphology in ceratopsian dinosaurs. Evolution, 29(2), 353–361. al.pone.0144148 https://doi.org/10.1111/j.1558-5646.1975.tb002​14.x Hedrick, B. P., & Dodson, P. (2013). Lujiatun psittacosaurids: Fastovsky, D. E., Weishampel, D. B., Watanabe, M., Barsbold, R., Understanding individual and taphonomic variation using 3D Tsogtbaatar, K. H., & Narmandakh, P. (2011). A nest of Protoceratops geometric morphometrics. PLoS ONE, 8(8), e69265. https://doi. andrewsi (Dinosauria, Ornithischia). Journal of Paleontology, 85(6), org/10.1371/journ​al.pone.0069265 1035–1041. Hernandez, C. E., Rodríguez-Serrano, E., Avaria-Llautureo, J., Inostroza- Forster, C. A., Sereno, P. C., Evans, T. W., & Rowe, T. (1993). A com- Michael, O., Morales-Pallero, B., Boric-Bargetto, D., … Meade, A. plete skull of Chasmosaurus mariscalensis (Dinosauria: Ceratopsidae) (2013). Using phylogenetic information and the comparative method from the Aguja Formation (Late Campanian) of West Texas. to evaluate hypotheses in macroecology. Methods in Ecology and Journal of Vertebrate Paleontology, 13(2), 161–170. https://doi. Evolution, 4(5), 401–415. org/10.1080/02724​634.1993.10011498 Holmes, R. B., Forster, C., Ryan, M., & Shepherd, K. M. (2001). A new Forth, C., Hedrick, B. P., & Ezcurra, M. D. (2016). Cranial ontogenetic vari- species of Chasmosaurus (Dinosauria: Ceratopsia) from the Dinosaur ation in the early saurischians and the role of heterochrony in the Park Formation of southern Alberta. Canadian Journal of Earth diversification of predatory dinosaurs. PeerJ, 4, e1589. https://doi. Science, 38(10), 1423–1438. org/10.7717/peerj.1589 Hone, D. W., Naish, D., & Cuthill, I. C. (2012). Does mutual sexual se- Gerber, S., & Hopkins, M. J. (2011). Mosaic heterochrony and lection explain the evolution of head crests and dinosaurs? Lethaia, evolutionary modularity: The trilobite genus Zacanthopsis 45(2), 139–156. as a case study. Evolution, 65(11), 3241–3252. https://doi. Hone, D. W. E., Wood, D., & Knell, R. J. (2016). Positive allometry for org/10.1111/j.1558-5646.2011.01363.x exaggerated structures in the ceratopsian dinosaur Protoceratops Goodwin, M. B., Clemens, W. A., Horner, J. R., & Padian, K. (2006). The andrewsi supports socio-sexual signaling. Palaeontologia Electronica, smallest known Triceratops skull: New observations on ceratopsid 19, 1–13. cranial anatomy and ontogeny. Journal of Vertebrate Paleontology, 26, Horner, J. R., & Goodwin, M. B. (2006). Major cranial changes during 103–112. Triceratops ontogeny. Proceedings of the Royal Society B: Biological Goswami, A. (2007). Cranial modularity and sequence heterochrony in Sciences, 273(1602), 2757–2761. mammals. Evolution and Development, 9(3), 290–298. https://doi. Horner, J. R., & Goodwin, M. B. (2008). Ontogeny of cranial epi-ossifica- org/10.1111/j.1525-142X.2007.00161.x tions in Triceratops. Journal of Vertebrate Paleontology, 28(1), 133–144. Goswami, A., & Polly, P. D. (2010). Methods for studying morphological Joshi, S. H., Prieto-Márquez, A., & Parker, W. C. (2011). A land- integration and modularity. In J. Alroy, & G. Hunt (Eds.), Quantitative mark-free method for quantifying biological shape variation. methods in paleobiology (pp. 213–243). Cambridge, UK: The Biological Journal of the Linnean Society, 104(1), 217–233. https://doi. Paleontological Society. org/10.1111/j.1095-8312.2011.01704.x Goswami, A., Polly, P. D., Mock, O. B., & Sánchez-Villagra, M. R. (2012). Kirkland, J. I., & DeBlieux, D. D. (2010). New basal centrosaurine ceratop- Shape, variance and integration during craniogenesis: Contrasting sian skulls from the Wahweap Formation (Middle Campanian), Grand marsupial and placental mammals. Journal of Evolutionary Biology, Staircase-Escalante National Monument, southern Utah. In M. J. 25(5), 862–872. https://doi.org/10.1111/j.1420-9101.2012.02477.x Ryan, B. J. Chinnery-Allgeier, & D. A. Eberth (Eds.), New perspectives 6308 | PRIETO-MÁRQUEZ et al. on horned dinosaurs: The Royal Tyrrell Museum Ceratopsian symposium Makovicky, P. J. (2008). Telling time from fossils: A phylogeny-based (pp. 117–140). Bloomington, IN: Indiana University Press. approach to chronological ordering of paleobiotas. Cladistics, 24(3), Klenovsek, T., & Jojic, V. (2016). Modularity and cranial integra- 350–371. https://doi.org/10.1111/j.1096-0031.2007.00184.x tion across ontogenetic stages in Martino’s vole, Dinaromys bog- Makovicky, P. J. (2012). Marginocephalia. In M. K. Brett-Surman, T. R. danovi. Contributions to Zoology, 85(3), 275–289. https://doi. Holtz, & J. O. Farlow (Eds.), The complete Dinosaur (2nd ed., pp. 527– org/10.1163/18759​866-08503002 550). Bloomington, IN: Indiana University Press. Klingenberg, C. P. (1998). Heterochrony and allometry: The analysis of Makovicky, P. J., & Norell, M. A. (2006). Yamaceratops dorngobiensis, a new evolutionary change in ontogeny. Biological Review, 73, 79–123. primitive ceratopsian (Dinosauria: Ornithischia) from the Cretaceous Klingenberg, C. P. (2008). Morphological integration and developmental of Mongolia. American Museum Novitates, 3530, 1–42. https://doi. modularity. Annual Review of Ecology, Evolution and Systematics, 39, org/10.1206/0003-0082(2006)3530[1:YDANP​C]2.0.CO;2 115–132. Makovicky, P. J., Sadleir, R., Dodson, P., Erickson, G. M., & Norell, M. Klingenberg, C. P. (2009). Morphometric integration and modularity in A. (2007). Life history of Protoceratops andrewsi from Bayn Zag, configurations of landmarks: Tools for evaluating a priori hypotheses. Mongolia. Journal of Vertebrate Paleontology, 27(Suppl. 3), 109A. Evolution and Development, 11, 405–421. Mallon, J. C., Holmes, R., Anderson, J. S., Farke, A. A., & Evans, D. C. Klingenberg, C. P. (2014). Studying morphological integration and (2014). New information on the rare horned dinosaur Arrhinoceratops modularity at multiple levels: Concepts and analysis. Philosophical brachyops (Ornithischia: Ceratopsidae) from the Upper Cretaceous Transactions of the Royal Society of London B: Biological Sciences, of Alberta, Canada. Canadian Journal of Earth Sciences, 51, 618–634. 369(1649), 20130249. https://doi.org/10.1098/rstb.2013.0249 Mallon, J. C., Holmes, R., Eberth, D. A., Ryan, M. J., & Anderson, J. S. (2011). Klingenberg, C. P., Barluenga, M., & Meyer, A. (2002). Shape analysis of Variation in the skull of Anchiceratops (Dinosauria, Ceratopsidae) symmetric structures: Quantifying variation among individuals and from the Horseshoe Canyon Formation (Upper Cretaceous) of asymmetry. Evolution, 56, 1909–1920. Alberta. Journal of Vertebrate Paleontology, 31(5), 1047–1071. Klingenberg, C. P., & Marugán-Lobón, J. (2013). Evolutionary covariation Mallon, J. C., Ott, C. J., Larson, P. L., Iuliano, E. M., & Evans, D. C. (2016). in geometric morphometric data: Analyzing integration, modular- Spiclypeus shipporum gen. et sp. nov., a boldly audacious new ity, and allometry in a phylogenetic context. Systematic Biology, 62, Chasmosaurine Ceratopsid (Dinosauria: Ornithischia) from the Judith 591–610. River Formation (Upper Cretaceous: Campanian) of Montana, USA. PLoS Knell, R. J., & Sampson, S. D. (2011). Bizarre structures in dino- ONE, 11(5), e0154218. https://doi.org/10.1371/journ​al.pone.0154218 saurs: Species recognition or sexual selection? A response to Maryanska, T., & Osmolska, H. (1975). Protoceratopsidae (Dinosauria) of Padian and Horner. Journal of Zoology, 283(1), 18–22. https://doi. Asia. Palaeontologia Polonica, 33, 133–181. org/10.1111/j.1469-7998.2010.00758.x McKinney, M. L., & McNamara, K. J. (1991). Heterochrony: The evolution of Lee, A. H., & Werning, S. (2008). Sexual maturity in growing dinosaurs Ontogeny. New York, NY: Plenum Press. does not fit reptilian growth models. Proceedings of the National McNamara, K. J. (2012). Heterochrony: The evolution of development. Academy of Sciences of the United States of America, 105(2), 582–587. Evolution: Education and Outreach, 5(2), 203–218. Lindgren, J., Currie, P. J., Siverson, M., Rees, J., Cederström, P., McNamara, K. J., & Long, J. A. (2012). Heterochrony in dinosaur evo- & Lindgren, F. (2007). The first neoceratopsian dinosaur re- lution. In M. K. Brett-Surman, T. R. Holtz, & J. O. Farlow (Eds.), The mains from Europe. Palaeontology, 50(4), 929–937. https://doi. complete Dinosaur (2nd ed., pp. 761–784). Bloomington, IN: Indiana org/10.1111/j.1475-4983.2007.00690.x University Press. Long, J. A., & McNamara, K. J. (1995). Heterochrony in dinosaur evolu- Meyer, A. (1987). Phenotypic plasticity and heterochrony in Cichlasoma tion. In K. J. McNamara (Ed.), Evolutionary change and Heterochrony managuense (Pisces, Chichlidae) and their implications for speciation (pp. 151–168). Chichester, UK: John Wiley and Sons Ltd. in cichlid fishes. Evolution, 41, 1357–1369. Long, J. A., & McNamara, K. J. (1997). Heterochrony: The key to dino- Morshhauser, E. M., You, H., Li, D., & Dodson, P. (2019). Phylogenetic saur evolution. In D. L. Wolberg, E. Stump, & G. D. Rosenberg (Eds.), history of Auroraceratops rugosus (Ornithischia: Ceratopsia) from the Dinofest International (pp. 113–123). Philadelphia, PA: Academy of Lower Cretaceous of Gansu Province, China. Journal of Vertebrate Natural Sciences. Paleontology, 38(Issue sup. 1, memoir 18), 111–147. Longrich, N. R. (2010). Mojoceratops perifania, a new chasmosaurine Nabavizadeh, A. (2018). New reconstruction of cranial musculature in ceratopsids from the Late Campanian of western Canada. Journal of ornithischian dinosaurs: Implications for feeding mechanisms and Paleontology, 84(4), 681–694. buccal anatomy. The Anatomical Record, 303(2), 347–362. https://doi. Longrich, N. R. (2014). The horned dinosaurs Pentaceratops and org/10.1002/ar.23988 Kosmoceratops from the upper Campanian of Alberta and implica- Nabavizadeh, A. (2019). Cranial musculature in herbivorous dinosaurs: tions for dinosaur biogeography. Cretaceous Research, 51, 292–308. A survey of reconstructed anatomical diversity and feeding mech- https://doi.org/10.1016/j.cretr​es.2014.06.011 anisms. The Anatomical Record, 303(4), 1104–1145. https://doi. Lull, R. S. (1908). The cranial musculature and the origin of the frill in org/10.1002/ar.24283 the Ceratopsian dinosaurs. American Journal of Science, 25(149), 399– O'Brien, D. M., Allen, C. E., Van Kleeck, M. J., Hone, D., Knell, R., 417. https://doi.org/10.2475/ajs.s4-25.149.387 Knapp, A., … Emlen, D. J. (2018). On the evolution of extreme struc- Maddison, W. P., & Maddison, D. R. (2018). Mesquite: A modular system tures: Static scaling and the function of sexually selected signals. for evolutionary analysis. Version 3.51. Retrieved from http://www. Animal Behaviour, 144, 95–108. https://doi.org/10.1016/j.anbeh​ mesqu​itepr​oject.org av.2018.08.005 Maiorino, L., Farke, A. A., Kotsakis, T., & Piras, P. (2015). Males resem- Olson, M. E., & Rosell, J. A. (2006). Using heterochrony to de- ble females: Re-evaluating sexual dimorphism in Protoceratops an- tect modularity in the evolution of stem diversity in the plant drewsi (Neoceratopsia, Protoceratopsidae). PLoS ONE, 10, e0126464. family Moringaceae. Evolution, 60(4), 724–734. https://doi. https://doi.org/10.1371/journ​al.pone.0126464 org/10.1111/j.0014-3820.2006.tb011​51.x Maiorino, L., Farke, A. A., Kotsakis, T., & Piras, P. (2017). Ösi, A., Butler, R. J., & Weishampel, D. B. (2010). A Late Cretaceous Macroevolutionary patterns in cranial and lower jaw shape of ceratopsian dinosaur from Europe with Asian affinities. Nature, ceratopsian dinosaurs (Dinosauria, Ornithischia): Phylogeny, mor- 465(7297), 466–468. https://doi.org/10.1038/natur​e 09019 phological integration, and evolutionary rates. Evolutionary Ecology Ostrom, J. H. (1964). A functional analysis of jaw mechanics in the dino- Research, 18, 123–167. saur Triceratops. Postilla, 88(3), 1–35. PRIETO-MÁRQUEZ et al. | 6309 Ostrom, J. H. (1966). Functional morphology of the ceratopsian dino- Slice, D. E. (2007). Geometric morphometrics. Annual Review of saurs. Evolution, 20(3), 290–308. Anthropology, 36(1), 261–281. https://doi.org/10.1146/annur​ Padian, K., & Horner, J. R. (2011). The evolution of “bizarre structures” in ev.anthro.34.081804.120613 dinosaurs: Biomechanics, sexual selection, social selection or species Sternberg, C. M. (1951). Complete skeleton of Leptoceratops grac- recognition? Journal of Zoology, 283(1), 3–17. ilis Brown from the Upper Edmonton Member on Red Deer River, Pagel, M. (1999). Inferring the historical patterns of biological evolution. Alberta. National Museum of Canada Bulletin, 123, 225–255. Nature, 401(6756), 877–884. https://doi.org/10.1038/44766 Tumarkin, A. R., & Dodson, P. (1998). A heterochronic analysis of enig- Raff, E. C., & Raff, R. A. (2000). Dissociability, modularity, evolv- matic ceratopsids. Journal of Vertebrate Paleontology, 18(Suppl. 3), ability. Evolution and Development, 2(5), 235–237. https://doi. 83A. org/10.1046/j.1525-142x.2000.00069.x Webster, M., & Zelditch, M. L. (2005). Evolutionary modifications of on- Revell, L. J. (2012). phytools: an R package for phylogenetic compara- togeny: Heterochrony and beyond. Paleobiology, 31, 354–372. tive biology (and other things). Methods in Ecology and Evolution, 3, Weishampel, D. B., & Horner, J. R. (1994). Life history syndromes, het- 217–223. erochrony, and the evolution of Dinosauria. In K. Carpenter, K. F. Revell, L. J. (2009). Size-correction and principal components for inter- Hirsch, & J. R. Horner (Eds.), Dinosaur eggs and babies (pp. 220–243). specific comparative studies. Evolution, 63, 3258–3268. New York, NY: Cambridge University Press. Rieppel, O. (1981). Die funktion des kragens der Ceratopsia. Willmore, K. E., Leamy, L., & Hallgrimsson, B. (2006). Effects of devel- Paläontologische Kursbücher, 1, 205–216. opmental and functional interactions on mouse canial variability Rohlf, F. J. (2016a). tpsUtil version 1.68. Retrieved from https://life.bio. through late ontogeny. Evolution and Development, 8(6), 550–567. sunysb.edu/morph​/ Xu, X., Forster, C. A., Clark, J. M., & Mo, J. (2006). A basal ceratopsian with Rohlf, F. J. (2016b). tpsDig2 version 2.26. Retrieved from https://life.bio. transitional features from the Late Jurassic of northwestern China. sunysb.edu/morph​/ Proceedings of the Royal Society B: Biological Sciences, 273(1598), Russell, L. (1935). Musculature and functions in the Ceratopsia. National 2135–2140. https://doi.org/10.1098/rspb.2006.3566 Museum of Canada Bulletin, 77(52), 39–48. Xu, X., Makovicky, P. J., Wang, X.-L., Norell, M. A., & You, H. (2002). A cer- Ryan, M. J., Holmes, R., & Russell, A. P. (2007). A revision of the Late atopsian dinosaur from China and the early evolution of Ceratopsia. Campanian centrosaurine ceratopsid genus Styracosaurus from the Nature, 416(6878), 314–317. https://doi.org/10.1038/416314a Western Interior of North America. Journal of Vertebrate Paleontology, You, H., & Dodson, P. (2004). Basal Ceratopsia. In D. B. Weishampel, P. 27(4), 944–962. Dodson, & H. Osmólska (Eds.), The Dinosauria (2nd ed., pp. 478–493). Ryan, T. J., & Semlitsch, R. D. (1998). Intraspecific heterochrony and life Berkeley, CA: University of California Press. history evolution: Decoupling somatic and sexual development in a You, H., Li, D., Ji, Q., Lamanna, M. C., & Dodson, P. (2005). On a new facultatively paedomorphic salamander. Proceedings of the National genus of basal neoceratopsian dinosaur from the Early Cretaceous of Academy of Sciences of the United States of America, 95, 5643–5648. Gansu Province, China. Acta Geologica Sinica, 79(5), 593–597. Sampson, S. D. (1995). Two new horned dinosaurs from the Upper Zelditch, M. L., Bookstein, F. L., & Lundrigan, B. L. (1992). Ontogeney Cretaceous Two Medicine Formation of Montana; with a phyloge- of integrating skull growth in the cotton rat Sigmodon fulviventer. netic analysis of the Centrosaurinae (Ornithischia: Ceratopsidae). Evolution, 46(4), 1164–1180. Journal of Vertebrate Paleontology, 15(4), 743–760. https://doi. Zelditch, M. L., & Carmichael, A. C. (1989). Growth and intensity of org/10.1080/02724​634.1995.10011259 integration through postnatal growth in the skull of Sigmodon Sampson, S. D., & Loewen, M. A. (2010). Unraveling a radiation: a re- fulviventer. Journal of Mammalogy, 70(3), 477–484. https://doi. view of the diversity, stratigraphic distribution, biogeography, and org/10.2307/1381419 evolution of horned dinosaurs. (Ornithischia: Ceratopsidae). In M. J. Zhao, X., Zhengwu, C., & Xu, X. (1999). The earliest ceratopsian from Ryan, B. J. Chinnery-Allgeier, & D. A. Eberth (Eds.), New Perspectives the Tuchengzi Formation of Liaoning, China. Journal of Vertebrate on Horned Dinosaurs: the Royal Tyrrell Museum Ceratopsian Symposium Paleontology, 19(4), 681–691. https://doi.org/10.1080/02724​ (pp. 405–427). Bloomington, IN: Indiana University. 634.1999.10011181 Sampson, S. D., Loewen, M. A., Farke, A. A., Roberts, E. M., Forster, C. A., Smith, J. A., & Titus, A. L. (2010). New horned dinosaurs from Utah S U P P O R T I N G I N FO R M AT I O N provide evidence for intracontinental dinosaur endemism. PLoS ONE, 5(9), e12292. https://doi.org/10.1371/journ​al.pone.0012292 Additional supporting information may be found online in the Sampson, S. D., Lund, E. K., Loewen, M. A., Farke, A. A., & Clayton, K. E. Supporting Information section. (2013). A remarkable short-snouted horned dinosaur from the Late Cretaceous (late Campanian) of southern Laramidia. Proceedings of the Royal Society B: Biological Sciences, 280(1766), 20131186. https:// How to cite this article: Prieto-Márquez A, Garcia-Porta J, doi.org/10.1098/rspb.2013.1186 Joshi SH, Norell MA, Makovicky PJ. Modularity and Sampson, S. D., Ryan, M. J., & Tanke, D. H. (1997). Craniofacial ontogeny in centrosaurine dinosaurs (Ornithischia: Ceratopsidae): Taxonomic heterochrony in the evolution of the ceratopsian dinosaur frill. and behavioral implications. Zoological Journal of the Linnean Society, Ecol Evol. 2020;10:6288–6309. https://doi.org/10.1002/ 121, 293–337. https://doi.org/10.1111/j.1096-3642.1997.tb003​4 0.x ece3.6361 Scannella, J. B., & Horner, J. R. (2010). Torosaurus Marsh, 1891, is Triceratops Marsh, 1889 (Ceratopsidae: Chasmosaurinae): Synonomy through on- togeny. Journal of Vertebrate Paleontology, 30(4), 1157–1168. Sereno, P. C., Xijin, Z., Brown, L., & Lin, T. (2007). New psittacosaurid high- lights skull enlargement in horned dinosaurs. Acta Palaeontologica Polonica, 52(2), 275–284.