The importance of the hadrosaurid dentition and feeding mode on their high diversity in the Late Cretaceous is emphasized by the fact that this clade and neoceratopsians appear to buck the trend of overall dinosaur decline through the last 40 million years of the Cretaceous, recently discovered from Bayesian analysis of overall dinosaurian speciation dynamics37.

Supertree

A species-level phylogenetic tree of Ornithopoda was constructed using source phylogenetic trees from the literature that resolve the relationships of major ornithopod subclades. Five recent trees were selected representing the scope of ornithopod diversity: Butler et al. (2011: Fig. 9)38, which resolves the position of basal ornithopods and provides a phylogenetic context for ornithopods within Ornithischia; Makovicky et al. (2011: Fig. 7B)9, which focuses on some of the more basal ornithopod subclades; Norman (2015: Fig. 51)39 and McDonald (2012: Fig. 1)40, which deal with the relationships of basal iguanodontians; and the recent hadrosaurid phylogeny by Prieto-Márquez (2016: Fig. 2)41. In addition, the rhabdodontid section of the phylogeny by Ősi et al.11 was added to the tree by Butler et al.38 to ensure full representation of this clade. Higher ornithischian outgroups to Ornithopoda taxa were replaced by one or two species and ceratopsians were removed from the tree. Euparkeria was chosen as the most distant outgroup to Ornithopoda42.

Each source tree was manually drawn in Mesquite (version 3.02)43 and subsequently input into the Supertree Toolkit44 for calculation of a supertree using matrix representation with parsimony (MRP). The MRP was subjected to maximum parsimony analysis using TNT 1.1 (ref. 45). The tree was produced using 10,000 replicates with the tree bisection and reconnection algorithm, due to the size of the dataset and efficiency of the search method46. This analysis resulted in 190 MPTs of 177 steps each, an optimum score found 676 times out of the 10,000 replicates. A strict consensus tree was computed from these MPTs (Fig. 1).

Dental characters

We scored 44 discrete dental morphological characters for 112 ornithopod and outgroup taxa (Supplementary Appendix 1). The traits represent dentary, maxillary, and, to a lesser extent, premaxillary teeth (absent in more derived ornithopods47). We chose these attributes as they document ornithopod feeding modes and the scope of morphological variation in these animals25. Specifically, 23 characters were culled from a character-taxon matrix that focused on hadrosaurids6 and eight and 13 were taken respectively from two studies38,39 that focused on non-hadrosaurid ornithopods. Taxon scoring was based on high-resolution photographs of fossil specimens, or using data from the scientific literature when photographs were not available, with a preference towards holotypes. The ontogeny of specimens was also recorded where available to ensure all specimens used were mature and therefore comparable.

Disparity analysis

To ensure that our data were compatible with principal coordinate analysis (PCO) and other disparity analyses, taxa with excessive missing data (e.g. solely known from postcranial remains or lacking teeth) were removed from consideration: Callovosaurus, Elrhazosaurus, Barilium, Tanius, Levnesovia, Nanyangosaurus, Jintasaurus, Jaxartosaurus, the Big Bend UTEP 37.7 hadrosaurid, Shantungosaurus, and Theiophytailia. Furthermore, all non-ornithopods were removed from the analysis.

The remaining taxa were placed in taxon-sorted and time-sorted bins. Taxon-sorted bins were based on the position of the taxa in the supertree, keeping monophyletic groups together where possible and avoiding bins with too few taxa to be representative. Five monophyletic clades were considered: jeholosaurids, rhabdodontids, hadrosaurids, saurolophine hadrosaurids and lambeosaurine hadrosaurids. However, keeping jeholosaurids as a taxon bin would leave the most basal ornithopod, Orodromeus, in a bin by itself. Orodromeus was therefore added to the ‘jeholosaurid’ bin. Taxon-sorted bins included: ‘jeholosaurids’, ‘pre-iguanodontids and post-jeholosaurids’, ‘rhabdodontids’, ‘post-rhabdodontids and dryosaurids’, ‘basal ankylopollexians and basal styracosterns’, ‘basal styracosterns’, ‘basal hadrosauroids’, ‘hadrosaurids’, ‘lambeosaurines’ and ‘saurolophines’. Some taxa that are traditionally placed within hadrosauroids are a part of ‘basal styracosterns’ rather than ‘basal hadrosauroids’, based on their position within the supertree.

The time-sorted bins represent ~10–15 Ma each and were based on first appearances of taxa according to the Paleobiology Database (https://paleobiodb.org). Five bins contain two stratigraphic stages, in order to equalise the numbers of taxa in each bin as much as possible and avoid unrepresentative bins with few taxa. The time bins are: Kimmeridgian-Tithonian (157.3–145 Ma), Berriasian–Valanginian (145–132.9 Ma), Hauterivian–Barremian (132.9–125 Ma), Aptian (125–113 Ma), Albian (113–100.5 Ma), Cenomanian–Turonian (100.5–89.8 Ma), Coniacian–Santonian (89.8–83.6 Ma), Campanian (83.6–72.1 Ma) and Maastrichtian (72.1–66 Ma).

Temporal and taxonomic disparity trends were assessed with two approaches. Firstly, the dental character-taxon matrix was used to calculate pairwise dissimilarity matrices based on generalised Euclidean distances (GED) and maximum observable distances (MOD), in the R package Claddis48. Morphological disparity in temporal and taxonomic bins was then calculated based on within-bin weighted mean pairwise disparity (WMPD), using both distance metrics. WMPD calculations place greater weighting on pairwise dissimilarities derived from more comparable characters, therefore preventing highly incomplete specimens from inflating or underestimating disparity49. In our second approach, we performed disparity calculations based on the positions of taxa in multivariate morphospace. The GED dissimilarity matrix was subjected to PCO to generate multivariate axes of variation, incorporating the Calliez negative eigenvalue correction. Disparity through time and within group bins was then quantified with binned PCO scores from the first 10 axes, using the sum of variances metric. Sensitivity tests were performed using different numbers of axes. In all disparity calculations, 95% confidence intervals were generated based on 10,000 bootstrap replicates.

To explore the morphological diversification of ornithopod dentitions, we examined taxon distribution in an empirical multivariate morphospace. The dental morphospace was based on a bivariate plot of PC axes 1 and 2, derived from the PCO described above. By plotting all taxa in a single pooled morphospace, we examined taxonomic trends and group separation. To explore the phylogenetic branching patterns within the dental morphospace, a pruned phylogeny with branch lengths (more below) was superimposed, producing a phylomorphospace. Temporal morphospaces were generated for the same bins as those used in the disparity analyses and convex hulls denote the overall area of morphospace occupied through time.

Permutation tests were implemented to test for taxonomic and temporal separation in morphopace. One-way PERMANOVA multivariate statistical tests were performed on the PCO scores for both time and taxa-sorted bins, to test for statistical significance in multivariate group means50, with a chosen alpha level of 0.05. The results were adjusted using Benjamini-Hochberg corrections51 to reduce bias from multiple comparisons52.

Evolutionary rates analysis

Rates of dental character evolution were examined with maximum-likelihood methods. All rate calculations were performed using the DiscreteCharacterRate function in the R package Claddis48. Our analyses aimed to determine if rates of dental evolution, based on the same 44 characters, were uniform in ornithopod evolution, or if particular branches were associated with significantly higher or lower rates of character change. In Claddis, rates are assessed by first time-scaling a phylogenetic tree and estimating ancestral character states. Likelihood ratio tests (LRTs) are then used to test the hypothesis that a particular branch shows significantly higher or lower rates than the pooled rate for the rest of the tree, based on the total numbers of character changes and branch durations53,54. An alpha threshold of 0.01 was used to assess significance, with Benjamini-Hochberg false discovery rate correction.

Rates calculations were performed on five random MPTs out of the 190 originally calculated for the supertree. A random number generator was used for this purpose and the selected trees were MPTs 22, 49, 50, 125 and 163. Firstly, the individual MPTs were pruned to remove non-ornithopods and taxa with excessive missing data (as described above for disparity analyses), leaving 89 tips. Stratigraphically calibrated branch lengths were calculated using the equal method29,55 by assigning taxa an age drawn randomly between their first appearance dates (FAD) and last appearance dates (LAD). For each MPT, the dating was repeated 20 times to test for consistency. To summarize rates results for all dating replicates across each MPT, we use pie charts positioned along each branch to illustrate the proportion of iterations that showed significantly high or low rates. To examine rates of character change in a temporal context, we generated time-series “spaghetti” plots using the same bins as the disparity analyses. These graphs illustrate per-lineage-million-years rates of character change in each interval and are calculated by dividing the sum of character changes by the summed duration of the branches within each time bin. This is repeated for each dating replicate of each MPT (100 topologies in total)56. It is important to acknowledge that such time-series rates calculation are influenced by diversity and specimen completeness in each interval – a greater number of complete fossils allows more opportunities to record character changes49.

Diversity metrics

Ornithopod diversity through time was assessed with two metrics. Firstly, we simply plotted taxonomic species diversity in each of the nine time bins used in the disparity calculations. Secondly, we calculated phylogenetic diversity estimates (PDE) incorporating both taxon occurrences and ‘ghost lineages’ inferred from the ornithopod supertree. PDE were made for the same nine time intervals. To account for phylogenetic and dating uncertainties within the supertree, unresolved nodes were randomly resolved and 100 dated topologies were generated. The median of the 100 topologies was plotted along with confidence intervals based on two-tailed 95% lower and upper quantiles.