Multifaceted disparity approach reveals dinosaur herbivory flourished before the end-Cretaceous mass extinction

Abstract

Understanding temporal patterns in biodiversity is an enduring question in paleontology. Compared with studies of taxonomic diversity, long-term perspectives on ecological diversity are rare, particularly in terrestrial systems. Yet ecological diversity is critical for the maintenance of biodiversity, especially during times of major perturbations. Here, we explore the ecological diversity of Cretaceous herbivorous dinosaurs leading up to the K-Pg extinction, using dental and jaw morphological disparity as a proxy. We test the hypothesis that a decline in ecological diversity could have facilitated their rapid extinction 66 Ma. We apply three disparity metrics that together capture different aspects of morphospace occupation and show how this approach is key to understanding patterns of morphological evolution. We find no evidence of declining disparity in herbivorous dinosaurs as a whole—suggesting that dinosaur ecological diversity remained high during the last 10 Myr of their existence. Clades show different disparity trends through the Cretaceous, but none except sauropods exhibits a long-term decline. Herbivorous dinosaurs show two disparity peaks characterized by different processes; in the Early Cretaceous by expansion in morphospace and in the Campanian by morphospace packing. These trends were only revealed by using a combination of disparity metrics, demonstrating how this approach can offer novel insights into macroevolutionary processes underlying patterns of disparity and ecological diversity.

Introduction

Studies of biodiversity in deep time are often focused on taxonomic diversity (Jackson and Johnson 2001). In contrast, how taxonomic diversity relates to a second aspect of biodiversity—ecological diversity—is little understood (Swenson 2011). Ecological diversity is broadly defined as the diversity of ecotypes, or niches, occupied (Magurran 1988). It is likely that species diversity is positively correlated with ecological diversity—but it is not obvious that this relationship will be linear. High diversity in species can be accommodated either through increased packing within occupied ecological space (niche partitioning) or through expansion into new niches (niche expansion [MacArthur 1965]). A recent study exploring both taxonomic and ecological diversity of bivalves during extinction events in the Phanerozoic found that great losses in taxonomic diversity did not equal losses in ecological diversity (Edie et al. 2018), supporting the notion that these processes can be decoupled. This is important, because ecological diversity may increase the robustness and resilience of ecosystems, even as species diversity declines. Resilience is defined as the ability of an ecosystem to recover after major perturbations (Elmqvist et al. 2003; Hooper et al. 2005). Ecosystems that have high resilience can adapt better and regain species diversity that was lost during times of environmental stress, while low resilience facilitates extinction cascades (Walker 1995; Elmqvist et al. 2003; Folke et al. 2004; Lindegren et al. 2016). Ecological diversity is therefore likely to be important for both the maintenance and recovery of biodiversity.

These topics are of great interest to paleontologists, yet studies with a long-term perspective on ecological diversity have been rare, especially in terrestrial systems. While many studies have explored morphological disparity patterns in clades (Erwin 2007; Hughes et al. 2013), only a subset has focused on ecologically important traits (such as teeth and jaws). Many studies also limit the inferences that can be made by calculating only a single metric, usually variance-based disparity. This is unfortunate, as using several metrics describing different aspects of morphological disparity is key to identifying patterns in ecological diversity (Foote 1999; Bush and Novack-Gottshall 2012; Novack-Gottshall 2016a, b). Novack-Gottshall (2016a) developed a conceptual framework describing the disparity trends we might expect during different processes of “ecospace” utilization, including contraction (loss of ecological diversity), redundancy, and niche partitioning. A key message is that different modes of evolution, such as loss of morphospace occupation versus packing, cannot be distinguished using only metrics describing variance or mean distance, for example, sum of variances (SoV) and mean pairwise distances (Fig. 1; Ciampaglio et al. 2001; Novack-Gottshall 2016a). Packing refers to how densely taxa occupy a given region of a morphospace. Close packing may be caused by redundancy in functional morphology, where the same morphology is shared among many taxa (Novack-Gottshall 2016a; Pigot et al. 2016). Yet these alternatives will lead to radically different interpretations (extinction on the one hand and packing on the other).

Figure 1 Different disparity metrics capture different aspects of disparity. A, Original morphospace occupation (shaded). A decrease in distance-based metrics (mean pairwise disparity [MPD] or sum of variances [SoV]) could be either due to an increase in packing in morphospace (B) or a decrease in the volume occupied in morphospace (C). Measuring morphovolume in addition to distance-based metrics allows us to separate these events, as morphovolume will not be affected by morphospace packing.

To gain a better understanding of how ecological diversity evolves, particularly in the face of environmental change and how this may relate to mass extinctions, there is a need for a long-term and detailed study of ecological diversity in terrestrial systems. Such approaches have previously been successful in uncovering the evolutionary dynamics of clades through time in marine biotas (Foote 1999; Korn et al. 2013; Hopkins and Smith 2015; Marx and Fordyce 2015; Stubbs and Benton 2016; Smithwick and Stubbs 2018). To this purpose, we explore the ecological diversity of herbivorous dinosaurs during the last 80 Myr of their reign, ending with the K-Pg mass extinction 66 Ma. We use dental and jaw morphological disparity as a proxy for diversity in feeding habits and, by extension, ecotypes. The rationale for this is that morphological variation in the jaws and dentition is expected to have a strong ecomorphological signal, with major adaptive innovations and modifications being intrinsically related to diet and feeding mode (Norman and Weishampel 1985; Barrett and Rayfield 2006; Cantalapiedra et al. 2013; Button et al. 2014; Gill et al. 2014; Stubbs and Benton 2016).

The K-Pg mass extinction event was one of the most severe in the history of life, leading to the complete extinction of nonavian dinosaurs, among many other groups. There is now abundant evidence that an asteroid impact coincided with the extinction (Alvarez et al. 1980; Schulte et al. 2010), and it is also concurrent with a large igneous event, the eruption of the Deccan Traps (Renne et al. 2015). The mass extinction was preceded by a period of extreme climate changes (Friedrich et al. 2012)—possibly weakening ecosystems before the final extinction. The extinction of the dinosaurs is therefore an excellent test case of whether lowered ecosystem resilience can catalyze major extinction events, as has been suggested for other mass extinctions (Roopnarine et al. 2007; Dick and Maxwell 2015). The lead-up to the extinction of the dinosaurs is of added interest, because it has been suggested that dinosaurs showed a broad-scale global decline in speciation rates through the last 40 Myr of the Cretaceous (Sakamoto et al. 2016), challenging the general consensus that there was no decline in diversity in their last 10 Myr (Fastovsky et al. 2004; Wang and Dodson 2006; Barrett et al. 2009; Lloyd 2012; Mitchell et al. 2012; Brusatte et al. 2015; Starrfelt and Liow 2016). Herbivorous dinosaurs are of particular interest, as they were hugely abundant and positioned at the bottom of the vertebrate food chain, and therefore changes in their diversity and ecology had the potential to create cascading effects through the ecosystem (Pringle et al. 2007; Vila et al. 2016). Moreover, herbivores are expected to encompass more species with dietary specializations than carnivores (Bernays and Graham 1988; Codron et al. 2016), so herbivorous dinosaurs should have been particularly sensitive to environmental perturbations (Davies et al. 2004; Clavel et al. 2011).

Several studies have investigated dinosaur disparity leading up to the K-Pg extinction, three of which exclusively focused on ecomorphology (Larson et al. 2016; Strickson et al. 2016; MacLaren et al. 2017). Larson et al. (2016) used a geometric morphometric approach, finding that the dental morphospace occupation of different maniraptoran clades remained stable throughout the Santonian–Maastrichtian interval. Strickson et al. (2016) used discrete dental characters of ornithopods and found that their disparity peaked in the Coniacian–Santonian, with a trough in the Campanian, followed by an increase in the Maastrichtian. MacLaren et al. (2017) explored disparity using geometric morphometrics and biomechanical characters of lower jaws from herbivorous dinosaurs and showed a steady increase in disparity throughout the Cretaceous, with no decline in the Maastrichtian. Based on discrete characters from the whole skeleton, Brusatte et al. (2012) detected declines in the disparity of ceratopsians and hadrosauroids, but not among smaller ornithopods, theropods, and sauropodomorphs, from the Campanian to the Maastrichtian. This result both contrasts with (Strickson et al. 2016) and agrees with (Larson et al. 2016) previously mentioned studies.

These recent studies have been valuable in exploring both patterns of dinosaurian adaptation and different data and methods, but it is hard to make direct comparisons between studies and build an overarching understanding of how the disparity of all herbivorous dinosaurs was changing during the Cretaceous. Conflicting results for the same clades (hadrosaurs and theropods) suggest that different approaches may yield different results (e.g., different sets of taxa, different sets of morphological characters, different disparity metrics). Further, it is reasonable to believe that disparity will vary depending on which part of the organism one measures, and whether geometric morphometrics or discrete characters are used (Mongiardino Koch et al. 2017). Time binning and the temporal scope of the analysis also matter, with some analyses having investigated only the latest Cretaceous (Brusatte et al. 2012), and others the whole of the Mesozoic, but with larger time bins (MacLaren et al. 2017). Timescales matter: peaks and troughs might be part of long-term fluctuations, which would be missed if the temporal focus is too narrow, while using large time bins may obscure short-term events and trends.

This study extends, and complements, earlier disparity studies on dinosaurs by exploring this theme using a larger data set, multiple disparity metrics, smaller time bins, including the great majority of herbivorous dinosaur taxa, and investigating global, regional, and clade-specific trends. By doing so, we uncover new detail and previously unaddressed aspects of dinosaur disparity dynamics preceding their extinction and demonstrate how a multifaceted approach to disparity is critical to understand the dynamics of morphological evolution.

Materials and Methods

Taxon Sampling

A list of all herbivorous dinosaur taxa (excluding theropods) during the Cretaceous period was created using the Paleobiology Database (http://paleodb.org) in conjunction with the literature. All valid species, according to Weishampel et al. (2004), plus new species described after that paper’s publication (referenced in Supplementary Data) with any cranial and/or dental material were included. After the character matrix was constructed, taxa with >70% missing data were excluded. This cutoff was chosen to exclude as few taxa as possible (26), yet remove extreme outliers. This was checked visually by inspecting morphospace plots for different cutoff values. The final matrix used for analysis consisted of 194 species, with 34% missing data overall (Supplementary Data 2).

We excluded herbivorous theropods for four reasons. (1) There is some debate about which theropods were herbivorous and which might have been omnivorous. (2) Character choice: our characters are focused on tooth morphology, which would be a poor estimate of beak disparity. Adjusting the characters to accommodate beak morphologies would have drastically increased missing data in the data set, both from nonapplicable dental characters in beaked taxa and the addition of nonapplicable beak-specific characters in toothed taxa. High percentages of missing data affect disparity analyses negatively; therefore, to maintain good quality we would need to cut a larger proportion of taxa, ultimately reducing the sample size of the study. (3) Rarity: herbivorous theropods constitute ~5% of herbivorous dinosaur species diversity in the Cretaceous (Butler and Barrett 2008; Butler et al. 2009). (4) Minimal impact: all key herbivorous theropod clades appear in the Early Cretaceous (therizinosaurs, ornithomimids, and oviraptorids) and continue through the Cretaceous, so if they were included, they would likely increase morphospace occupation, and increase expansion in the Early Cretaceous, but not change the main trends during the Cretaceous. Other methodological approaches would be more suitable for including beaked forms, such as geometric morphometrics measuring the shape of the jaw (although theropods were also excluded from the largest such study, by MacLaren et al. [2017] on the basis of lack of complete specimens).

Morphological Characters and Scoring

Eighty-eight morphological characters, encompassing dental and jaw morphology, were used in the analysis (Supplementary Data 1). Forty-two of the dental characters were adopted from Strickson et al. (2016), with minor adjustments to the scoring of ornithopods. Another 46 were collected from sources chosen to cover morphological variation in all groups included and to add characters describing jaw morphology. We used photographs and descriptions of specimens in the literature for scoring of characters (listed and referenced in Supplementary Data 2).

Timescale and Time Bins

The Cretaceous timescale used in this study is based on Gradstein et al. (2012) and the updates in the 2016 International Chronostratigraphic Chart. Using this timescale, dinosaur-bearing geological formations were given start and end dates based on Weishampel et al. (2004) or, if available, more recent publications (Supplementary Data 2). The first and last appearance dates of each taxon were then assigned based on these formation ages. We conducted a sensitivity test to evaluate the effect of different time-bin divisions (geological substages, geological stages and equal-length time bins of 10, 5, and 3 Myr; Supplementary Figs. 1, 2). By evaluating the effects of these divisions, we constructed a modified series in which long time bins were divided into substages (Albian and Campanian) and two short bins with particularly few taxa, Coniacian and Turonian, were combined. This division allowed for time bins that were most equal in sample size and duration (∼5 Myr), without either increasing error bars dramatically or losing temporal information in the data.

Disparity Analyses

All analyses were conducted in R v. 3.2.2. (R Core Team 2015). For constructing a distance matrix, we calculated generalized Euclidean distances (GED) (Wills 1998, 2001) and maximum observable distances (MORD) using the package Claddis (Lloyd 2016). GED has the benefit of always returning a complete matrix (required for ordination) by filling in a mean distance (based on calculable distances) into missing fields. MORD rescales the calculable distances based on the information available, by dividing each distance by the maximum distance calculated from the observed characters (Lloyd 2016). Using both distance matrices, we then calculated preordination (weighted mean pairwise disparity [WMPD]) and postordination (SoV, principal coordinates analysis [PCO] convex-hull volume) disparity metrics.

Mean pairwise disparity measures how similar a given number of taxa are, that is, how closely packed they are in morphospace. The WMPD is different, in that it gives more weight to pairwise distances based on more comparable characters (Close et al. 2015). To test for effects of sample size and produce 90% confidence intervals, we applied a bootstrap analysis (1000 replicates with replacement). Confidence intervals were based on percentiles of the bootstrap distribution.

For postordination metrics (SoV, PCO volume) we applied a PCO to the GED and MORD distance matrices. A common problem when PCO is applied to a distance matrix based on discrete characters is the creation of negative eigenvalues. Negative eigenvalues in a PCO can be corrected by adding a constant of equal size to the largest negative eigenvalue to all distances (e.g., Cailliez correction; Cailliez 1983). We explored how well “corrected” and “uncorrected” scores represent the original distance matrix, and how this affected the SoV statistic (Supplementary Note 2, Supplementary Figs. 3, 6). We note that the variation explained by PC 1 and PC 2 is notably higher when no correction is used. Negative eigenvalues need not be a problem for the representation of the first few PC axes if the greatest absolute negative value is smaller than the greatest positive value (Legendre and Legendre 1998). This condition is met in our data set. We discuss the fit of our ordination data to the original distances in the Supplementary Material. In summary, we found a high correlation (95%) between distances after ordination axes and original distances when no correction was applied. This suggests that results from SoV and WMPD should be comparable.

The SoV statistic measures the spread of taxa in morphospace. In this aspect, it is like WMPD. SoV can be measured for any number of axes, and how many are included varies between studies. We included 30 axes based on the R ² values of the correlation between distances in the PCO space with the original distances in the GED or MORD matrices (Supplementary Fig. 3). At the point where the curve plateaus, additional axes do not add to the variation explained, and this can hence be used as the cutoff point. A bootstrap method (1000 replicates) was applied to estimate 90% confidence intervals. Following Kotrc and Knoll (2015a), we also quantified the strength of association between all discrete characters and the PCO axes using Cramér’s V coefficient (Supplementary Note 1, Supplementary Fig. 4).

We quantified the total amount of morphospace occupied through time by measuring the PCO convex-hull volume. PCO volume based on a convex hull is sensitive to outliers (Kotrc and Knoll 2015b), which can create large peaks in volume from single data points. This can be adjusted by changing the alpha shape parameter, which allows large empty spaces (with few data points) to be removed. The alpha shape parameter regulates the radius by which such spaces are “scooped out.” We calculated volumes for three dimensions (PC 1–3) using different values of alpha. The range of values to be tested was determined by visually inspecting morphospace plots for different alpha values. We also calculated raw morphospace volumes. For volumes and alpha shapes, we calculated 90% confidence intervals based on percentiles of a bootstrap distribution (1000 replicates).

To measure the amount of morphospace packing and expansion for each time bin, we used the niche packing “flexible” metric (NP flexible) created by Pigot et al. (2016). This metric is calculated by estimating the number of species in a volume A1 that can be packed within another volume A2. A1 and A2 are sequential time bins, so, for example, to calculate the percentage packing in the Campanian, A1 is the Campanian and A2 the Santonian. A greedy algorithm sequentially removes species from A1, each time choosing the data point that incurs maximum loss of morphovolume. Once the modified A1 volume reaches the same volume as A2, the taxa remaining in A1 are considered to contribute to packing, and the removed taxa are considered to expand trait space. In this way, percentages for expansion and packing can be calculated. The benefit of this algorithm is that it compares packing in the absolute volumes of A1 and A2, irrespective of their positions in morphospace.

Results

Distance-based Disparity Trends

Distance-based metrics (WMPD and SoV) highlight three key trends in herbivorous dinosaur disparity, revealing an Aptian peak, a Campanian low, and a stable/increasing trend in the Maastrichtian (Fig. 2A,B,D,E). As expected, both distance-based metrics, within-bin WMPD (Fig. 2A,D) and the SoV (Fig. 2B,E), show almost identical trends. This confirms that our ordination has not distorted the original distances in any major way (we explore this more in depth in the Supplementary Material). They are also unchanged when using two approaches for calculating morphological distances from the discrete character data, the GED (Fig. 2D,E) and MORD (Fig. 2A,B). There is a major discrepancy between the two distance measures in the time from the Albian to the Turonian–Coniacian, when the MORD metric records a sharp rise, but the GED metric suggests unchanged or decreasing disparity. The Aptian peak is not pronounced using the MORD metric, but it does appear using smaller time bins (Supplementary Figs. 1, 2). Apart from these key features, distance-based disparity remains relatively stable during the Cretaceous.

Figure 2 Disparity of herbivorous dinosaur teeth and jaws through the Cretaceous. A, Weighted mean pairwise disparity (WMPD), maximum observable distances (MORD)-based disparity; B, sum of variances (SoV), MORD-based disparity; C, morphovolume (principal coordinates analysis [PCO] convex-hull volume) based on MORD: upper graph raw volumes, lower graph alpha shapes; D, WMPD, generalized Euclidean distances (GED)-based disparity; E, SoV, GED-based disparity; F, morphovolume based on GED: upper graph raw volumes, lower graph alpha shapes. Shaded envelopes represent 90% confidence intervals based on 1000 bootstrap replicates. Stratigraphic abbreviations: Be, Berriasian; V, Valanginian; H, Hauterivian; B, Barremian; Ap, Aptian; Al, Albian; Ce, Cenomanian; T, Turonian; C, Coniacian; S, Santonian; Cm, Campanian; Ma, Maastrichtian.

MorphospaceOccupation

Phenotypic similar-ity between taxa can be visualized in a morphospace based on the jaw and dental character data (Fig. 3). The first two axes represent 32% of the total variation, high for a morphospace built on a discrete character data set (e.g., Brusatte et al. 2012; Halliday and Goswami 2016), but low compared with a landmark-based data set. However, a strong correlation (R=0.86) between PC scores on the first and second axes and original distances in the distance matrix confirms that this is a good representation of the relative placement of taxa in morphospace (Supplementary Fig. 3A). Herbivorous dinosaur clades generally form clusters within the morphospace, with some groups occupying distinct regions, such as hadrosauroids and the highly divergent sauropods (Fig. 3). The separation between clades in morphospace is significant between all groups except non-hadrosauroid ornithopods and neoceratopsians, using nonparametric multivariate analysis of variance (NPMANOVA) for pairwise comparison (Supplementary Table 1). We note, however, that NPMANOVA can be misleading when dispersion differs greatly in groups: a significant result can be confounded by difference in spread, rather than solely by difference in position (Anderson 2001). In PC 1–PC 2 morphospace, ornithischians with less derived teeth, such as early ornithopods, ankylosaurs, psittacosaurids, and pachycephalosaurids, converge, forming a group. More derived ceratopsians and ornithopods diverge and approach the area occupied by hadrosauroids. Statistical tests, based on Cramér’s V coefficient, confirm that most characters contribute to the positioning of taxa in morphospace, and there is no bias for certain individual characters or categories of characters (e.g., dentary and maxillary teeth or jaw-related features; Supplementary Fig. 4).

Figure 3 Jaw and dental morphospace occupation for groups of herbivorous dinosaurs. PC 1–PC 2 explains 32% of variance, and the correlation between pairwise distances from the first two PC axes and squared pairwise distances from the original matrices is 0.86 (Supplementary Note 1). A, Total morphospace occupation with convex hulls drawn around the major clades. B, Temporal trends in morphospace occupation through the Cretaceous. Both plots are based on the generalized Euclidean distances metric, and expanding taxa (shaded in B) are based on a two-dimensional packing analysis. Stratigraphic abbreviations: Ber, Berriasian; Val, Valanginian; Hau, Hauterivian; Bar, Barremian; Apt, Aptian; Alb, Albian; Cen, Cenomanian; Tur, Turonian; Con, Coniacian; San, Santonian; Cmp, Campanian; Maa, Maastrichtian. Silhouettes from PhyloPic.org: upper left image by Mathew Wedel, under Creative Commons Attribution 3.0 unported license (https://creativecommons.org/licenses/by/3.0/); upper middle image by Raven Amos, under Creative Commons Attribution 3.0 unported license (https://creativecommons.org/licenses/by/3.0/); upper right image and lower image licensed under Public Domain Dedication 1.0 license (https://creativecommons.org/publicdomain/zero/1.0/).

Morphovolume Disparity Trends

Morphovolume (PCO convex-hull volume) values during the Cretaceous show both similarities and marked differences compared with distance-based metrics (Fig. 2C–F). The volume of morphospace occupied is high in the Aptian, in agreement with distance-based metrics, but it is also high in the Campanian, when distance-based metrics give low values. These patterns are unchanged using different distance metrics. Morphospace volumes are sensitive to outliers, and we have accounted for this by changing the alpha shape parameter, which removes large areas of empty space, and by applying a bootstrap method. Alpha volumes clearly reduce the variability in calculated volumes, suggesting this approach is effective in removing the impact of outliers. Even with confidence intervals and for low alpha values, the peaks in the Aptian–Albian and Campanian remain (Fig. 2C–F), suggesting they are not caused by outliers. Patterns are similar for both distance metrics, but the MORD metric shows high variability in the Turonian–Coniacian bin, probably reflecting the peak recorded by distance-based metrics. However, this increase does not persist at decreasing alpha values. Similar patterns are seen when six morphospace axes are included (hypervolumes; Supplementary Fig. 5).

Morphospace Packing and Expansion

We define morphospace packing as an increase in species within the same volume as recorded from the previous time bin (which will by necessity increase the density of taxa). If sample size does not increase, species can still expand, or lack expansion, but the latter does not exemplify packing. Our morphospace packing analysis shows that the Aptian and Campanian peaks in morphovolume are characterized by different processes. Morphospace packing was elevated in the Campanian, while the Aptian exhibits a combination of both packing and expansion (Fig. 4). Packing is high in the Campanian, both in absolute numbers and as a percentage compared with expansion. Apart from the Aptian, expansion is high in the Valanginian, Barremian, and Santonian. The Hauterivian, Albian, Cenomanian, Turonian–Coniacian, and Maastrichtian all show drops in sample size and a lack of morphospace expansion.

Figure 4 Morphospace packing through time. Each value is calculated as the volume in the younger bin accommodated in previous bin. A, Proportion of unique taxa contributing to packing/expansion for two dimensions; B, number of unique taxa contributing to packing/expansion for two dimensions. “No packing/expansion” refers to bins where sample size decreases, therefore they do not constitute packing but a lack of expansion in morphospace. All volume estimates are based on the generalized Euclidean distances metric. Stratigraphic abbreviations as in Fig. 2.

Regional Trends

Our data set allows division between the most well sampled geographical regions, the Eurasian and North American continents (Fig. 5). The distance-based metrics for Eurasia lack the Aptian high disparity but show the same low Campanian disparity. The North American trend shows both high disparity in the Aptian and low disparity in the Campanian, in agreement with the global trend. In contrast to the global trend, however, there is an increasing trend from the Albian to the Turonian–Coniacian. Morphovolumes for both continents retain the Aptian–Albian and Campanian peaks in disparity, mirroring the global trend (similar trends are recorded for MORD-based distances; see Supplementary Fig. 8).

Figure 5 Regional disparity trends. Disparity trends for Eurasian taxa: A, Weighted mean pairwise disparity (WMPD), maximum observable distances (MORD) based; B, WMPD, generalized Euclidean distances (GED) based; C, principal coordinates analysis (PCO) volume, GED based. Disparity trends for North American taxa: D, WMPD, MORD based; E, WMPD, GED based; F, PCO volume, GED based. Shaded envelopes represent 90% confidence intervals based on 1000 bootstrap replicates. Stratigraphic abbreviations as in Fig. 2.

Clade-Specific Trends

Herbivorous dinosaur clades show distinct disparity dynamics (Fig. 6). Both distance metrics generally show identical patterns; therefore, only results based on the GED distance metric are discussed here (see Supplementary Information for MORD analyses). Neoceratopsians exhibit a decrease in distance-based disparity in the Santonian and then increase throughout the rest of the Late Cretaceous, with morphovolumes increasing in tandem in the Campanian, then decreasing in the Maastrichtian (Fig. 6A). Hadrosauroids dramatically decrease in distance-based disparity in the middle Campanian, and then recover previous levels of disparity in the Maastrichtian (Fig. 6B). Hadrosauroid morphovolumes, in contrast, increase through the Campanian, peaking in the late Campanian, with a decrease in the Maastrichtian. Non-hadrosauroid ornithopods increase in disparity in the Barremian, then remain on a stable level of disparity through the Aptian–Albian and Maastrichtian (Fig. 6C). High morphovolume is recorded in the Aptian. Pachycephalosaurids show stable distance-based disparity in the Campanian and Maastrichtian and increasing morphovolumes during the same interval (Fig. 6D). Sauropods exhibit high distance-based disparity values in the Aptian and middle to late Albian, which then plummet toward the Maastrichtian (Fig. 6E). Similarly, sauropod morphovolumes peak in the Aptian, then decrease through the rest of the Cretaceous. Ankylosaurs show the lowest levels of distance-based disparity in the Santonian, compared with their relatively stable and higher disparity in the Early Cretaceous, but show an increase in the Maastrichtian (Fig. 6F). Ankylosaur morphovolume also remains stable through the Cretaceous, with a slightly increasing trend toward the late Campanian. Finally, psittacosaurids record the highest levels of distance-based disparity in the Hauterivian, a decrease in the Aptian, but recovery in the Albian, albeit with slightly lower than Hauterivian values (Fig. 6G). In contrast, psittacosaurid morphovolume peaks in the Aptian.

Figure 6 Clade-specific disparity trends. Trends in principal coordinates analysis (PCO) volume (darker shading) and weighted mean pairwise disparity (WMPD) (lighter shading) through the Cretaceous for: A, neoceratopsians; B, hadrosauroids; C, non-hadrosauroid ornithopods; D, pachycephalosaurids; E, sauropods; F, ankylosaurids; G, psittacosaurids. All values based on the generalized Euclidean distances metric (see Supplementary Information for maximum observable distances–based disparity). Shaded envelopes represent 90% confidence intervals based on 1000 bootstrap replicates. Stratigraphic abbreviations as in Fig. 2. Silhouettes from PhyloPic.org: A, By Raven Amos, under Creative Commons Attribution 3.0 unported license (https://creativecommons.org/licenses/by/3.0/); B and E, under Public Domain Dedication 1.0 license (https://creativecommons.org/publicdomain/zero/1.0/); C, by Jaime Headden, under Creative Commons Attribution 3.0 unported license (https://creativecommons.org/licenses/by/3.0/); D, by Emily Willoughby, under Creative Commons Attribution-ShareAlike 3.0 unported license (https://creativecommons.org/licenses/by-sa/3.0/); F, by Scott Hartman, under Creative Commons Attribution 3.0 unported license (https://creativecommons.org/licenses/by/3.0/); G, by Pete Buchholz, licensed under Creative Commons Attribution-ShareAlike 3.0 unported license (https://creativecommons.org/licenses/by-sa/3.0/).

Discussion

Disparity Trends in the Cretaceous

We find three key results: (1) an Aptian peak in disparity characterized by expansion and packing; (2) a second disparity peak in the Campanian characterized by intense morphospace packing; and (3) a thinning out of morphospace in the Maastrichtian. These results paint a picture of a largely resilient dinosaurian ecosystem that was able to adapt and flourish in the face of climate-driven environmental changes.

Morphospace expansion in the Aptian suggests that this was a time of diversification of dental and jaw morphologies, particularly among sauropods and early hadrosauroids (Fig. 3). This includes Nigersaurus taqueti, which is the only sauropod known to have evolved a specialized “dental battery” (Sereno and Wilson 2005), and early hadrosauroids (Altirhinus kurzanovi, Jinzhousaurus yangi) with primitive versions of dental batteries. Dental batteries would later become more advanced in hadrosaurids and were likely a key factor in the success of this group in the Late Cretaceous (Norman and Weishampel 1985; Strickson et al. 2016).

Morphospace saturation coupled with high morphovolume in the Campanian indicates that this was not a time when dinosaurs were declining in morphological disparity. On the contrary, dinosaurs were flourishing, with the Late Cretaceous seeing the radiation of successful groups such as hadrosaurids and neoceratopsians. This exemplifies how low disparity values from distance-based metrics can be caused by increased packing and high morphodensity, rather than a loss of morphovolume (Novack-Gottshall 2016a). New species filled in and saturated the morphospace already occupied, a process that was particularly seen in hadrosauroids, which became very diverse and abundant in the latest Cretaceous but did not show wide variation in jaw and dental morphology (Strickson et al. 2016). Very low distance-based disparity coupled with high morphovolumes in the hadrosauroids confirms this (Fig. 6B).

The decreased morphodensity and morphovolume in the Maastrichtian represents an effect opposite to what happened in the Campanian: a thinning in morphospace. Whether this is a true trend or an edge effect is unclear. Either way, our results do not indicate any long-term decline in dinosaur disparity preceding the end-Cretaceous extinction, in agreement with previous disparity studies (Brusatte et al. 2012; Larson et al. 2016; Strickson et al. 2016; MacLaren et al. 2017).

An exception can be seen in sauropods, which appear to have truly been in long-term decline. A combined decrease in distance-based disparity and volume suggests declining morphological diversity (as opposed to morphospace packing) (Fig. 6E), in agreement with the known disappearance of broad-crowned sauropods in the Late Cretaceous (Chure et al. 2010). Sauropods, large herbivores with slow generation times, might have been more sensitive to extreme environmental changes during the mid-Cretaceous (Gaston and Blackburn 1995; Purvis et al. 2000; Cardillo et al. 2005).

In sum, our disparity results show a dynamic pattern of waxing and waning of clades—sauropods flourishing in the Aptian–Albian then diminishing in the Late Cretaceous, while other clades rose to prominence (e.g., hadrosaurids, ceratopsians). Collectively, herbivorous dinosaurian ecosystems show no signs of “weakening” preceding the Cretaceous mass extinction.

Contrasting Diversity and Disparity in Hadrosauroids

Empirical studies suggest that ecosystem species richness is positively correlated with specialization, suggesting that ecosystems with high species diversity mediate the increase in species numbers mainly through niche partitioning rather than niche expansion (Case 1981; Belmaker et al. 2012; Pigot et al. 2016; Pellissier et al. 2018). Our results show that a similar pattern might be true for one of the most species-rich groups of dinosaurs, the hadrosauroids. Once advanced dental batteries had evolved in early hadrosauroids, these structures required little change to ensure huge success for the clade. Our measures are currently limited to estimating morphospace packing, but an interesting further avenue of research is to test whether the taxa filling morphospace are doing so in an even manner or are more densely populating certain areas. Novack-Gottshall (2016a) offers a framework to test these hypotheses, by measuring the degree of redundancy versus partitioning. Redundancy refers to the occupancy of niches identical to those previously occupied by existing species.

Though suggestive, a high degree of redundancy/partitioning in our morphospace does not necessarily equal niche partitioning. “Niche space” has many more dimensions than differing feeding modes. For example, differences in body size incur differences in predation pressures (with higher predation pressures on small animals). In modern African savanna ecosystems, such differences can lead to differing habitat preference (open/vegetated) and, consequently, differing resource use among herbivores (Riginos and Grace 2008). If hadrosauroids are disproportionally spread out along some other axis of niche space, it is possible that they might not have exhibited greater niche partitioning than other groups. A more detailed analysis should also include a spatial dimension, as hadrosaurs on different continents could inhabit similar niches without competing.

Differences to Previous Analyses

Our analyses reveal differences in temporal, regional, and clade-specific trends compared with earlier studies. MacLaren et al. (2017) did not identify a disparity peak in the Aptian, as in this study, but instead their distance-based disparity value rises steadily throughout the Cretaceous. Neither did they recover a trough in the Campanian. Their study used 10 Myr time bins and a smaller sample, which might be responsible for the discrepancy. Large time bins can alter disparity curves substantially, as is evident from adjusting time bins for our own data, and this effect has also been demonstrated previously (Deline and Ausich 2011). Applying 10 Myr bins, both distance-based metrics show a steadily decreasing trend starting in the Early Cretaceous (Supplementary Figs. 1, 2). Another reason might be the difference in method and morphology measured, as MacLaren et al. (2017) used geometric morphometrics and functional traits relating to lower jaw shape, while our study used discrete characters with a focus on both dental and jaw morphology. Brusatte et al. (2012) found global decreases in hadrosaur, neoceratopsian, pachycephalosaur, and ankylosaur disparity from the Campanian to the Maastrichtian, and an increase in sauropod disparity. This is the exact opposite of the results reported here: sauropods show a decrease and all other groups record an increase in distance-based disparity from the Campanian to the Maastrichtian (Fig. 6). This is intriguing and suggests that dental and jaw disparity, linked more tightly to feeding ecology, might show a different pattern from the disparity of the entire skeleton. Further studies should explore the disparity of the entire skeleton, to clarify whether and how overall disparity dynamics differ from those of dental and jaw morphology alone. Also in contrast to Brusatte et al. (2012), our main trends are shared on the Eurasian and North American continents (Fig. 5). This further strengthens the argument that herbivorous dinosaurs worldwide were not experiencing long-term declines (Brusatte et al. 2015) and suggests that the trends recovered herein are not restricted to regional patterns.

The overall message of this study is different from that of Sakamoto et al. (2016), who found declines in speciation rates, but the results need not be contradictory. Sakamoto et al. (2016) did find increasing speciation rates in neoceratopsians and hadrosaurids, indicating radiations of these groups, which agrees with the morphospace packing results presented here. Furthermore, we present data suggesting sauropods were the only group declining in disparity, apparently coinciding with declining speciation rates. This is not to say that speciation rates and disparity must be coupled, but it is interesting to note that in the case of nonavian dinosaurs, both trends seem to co-occur, with reduced speciation rates matched by low morphological diversity in sauropods, and high speciation rates matched by high volume-based disparity (but low distance-based disparity due to packing) in hadrosauroids. Therefore, the overall picture of dinosaur evolutionary dynamics preceding the extinction is one of faunal turnover, with some clades prospering and others declining, rather than a collective decline. In addition, the focus here on morphospace volume and packing provides different insights than a study of coupled speciation and lineage extinction rates.

Potential Biases

When using fossil data, the relative proportion of bias and true signal can be difficult to estimate. Contrary to MacLaren et al. (2017), a decrease in disparity coinciding with an increase in sample size does not prove that disparity measures are not biased by sample size. A decrease in distance-based disparity could be due to either reduced morphospace occupation or extensive packing (high morphodensity) (Fig. 1). Packing in morphospace is likely to be higher when a larger sample is included, and thus increased sample size can also cause drops in disparity. Here, we apply a bootstrap method to simulate the effect of random sampling, as done elsewhere (e.g., Halliday and Goswami 2016; Strickson et al. 2016; MacLaren et al. 2017). We also use alpha shape volumes, following Kotrc and Knoll (2015a), to remove large areas of empty morphospace in volume estimates. The volume analyses show that applying alpha shapes strongly reduces the confidence interval, suggesting it is a useful method to remove large peaks caused by outliers. We also use two dissimilarity metrics (MORD and GED) to investigate how robust our results are to different methodological approaches.

How effective are these methods for identifying spurious peaks? The Turonian–Coniacian peak using the MORD metric is significantly higher than any other point. However, it is not recorded by the GED metric. This suggests that this peak in disparity is not a stable result, as it is not reproducible using two alternative metrics. Similarly, the same Turonian–Coniacian peak does not persist at low alpha levels in the volume estimates. The sample size in this bin is low, and upon inspection of morphospace occupation through time, it appears that the MORD metric estimates higher disparity between neoceratopsians, sauropods, and hadrosauroids than does the GED metric. This suggests that the Turonian–Coniacian peak is caused by a few very disparate taxa. Large areas of empty morphospace are often a sign of bias caused by missing data, and the ability of the MORD metric to handle missing data that are not randomly distributed (Lloyd 2016) should be investigated. The implications of this are that we should be cautious to interpret the Turonian–Coniacian peak (and trough) as a real phenomenon—it is probably a result of poor sampling in this bin (Supplementary Fig. 7). However, the Aptian and Campanian patterns, which we have focused on as our main result, persist with MORD and GED metrics, both for low alpha values and with the bootstrap confidence interval. Our bootstrap and alpha shape analysis do reveal, however, that the very high volume in the latest Campanian is unstable. We cannot be certain that the Campanian had a higher morphovolume than the Aptian—but high volumes in the Aptian and the Campanian, as compared with other time bins, are well supported. We suggest our data set is reliable enough to illuminate larger-scale patterns.

The clade-specific analyses, particularly for morphovolume, are influenced by low sample size. Fine temporal trends can therefore not be discerned, and peaks should be interpreted with caution. Higher sample sizes are achieved for neoceratopsians and hadrosauroids, and these are the only clades that can be interpreted with certainty. Both clades are well sampled, and it is interesting to note that in the Campanian bin, ceratopsians exhibit coupled increase in volume and variance-based disparity (expansion), while hadrosauroids show an increase in volume but decrease in the variance metric.

Ultimately, disparity analyses can only record the disparity of the known fossil record. Therefore, although we have shown that some trends are robust to random sampling, we cannot account for the fact that some geological intervals are relatively poorly sampled and will underestimate the true disparity more than other bins. PCO volumes will naturally be closely linked to sample size when sample size is low, as each species will add/remove some volume (as opposed to measuring variance). For very large data sets, this effect is likely to be less pronounced, because the effect of sample size tapers out (Ciampaglio et al. 2001). The extent and effect of sampling bias on the dinosaur fossil record is a subject of much debate, and a number of methods and models have been developed to estimate bias (Barrett et al. 2009; Brusatte et al. 2015; Sakamoto et al. 2016; Starrfelt and Liow 2016). Subsampling by standardizing to very low sample sizes likely causes much of the real signal to be wiped out, resulting in a nearly flat diversity curve (Close et al. 2018). The phylogenetic diversity estimate accounts for ghost lineages and, when applied to the dinosaur fossil record, mostly follows that of the raw taxonomic count; therefore, depending on position, it either confirms the major diversity trends or demonstrates that the phylogenetic diversity estimate cannot remove biases effectively. The majority of recent diversity studies have been based on the residuals model (Barrett et al. 2009; Upchurch et al. 2011; Lloyd 2012; Mannion et al. 2012; Brusatte et al. 2015), which has since been shown to be based on incorrect statistical assumptions and is therefore flawed (Sakamoto et al. 2017; Dunhill et al. 2018). A different approach is to quantify the sampling rate of each bin to estimate true richness, the TRiPS method, as applied by Starrfelt and Liow (2016). Their results confirm diversity peaks in the Aptian and Campanian and suggest that the Aptian–Albian and Campanian–Maastrichtian are particularly well sampled (Starrfelt and Liow 2016). It is therefore reasonable to consider the Albian–Aptian and Campanian–Maastrichtian disparity values as good estimates, while the other bins should be treated with more caution, particularly for volume estimates. These two peaks are also the most interesting from the morphospace packing/expansion perspective, with the combined expansion/packing in the Albian and intense packing in the Campanian. Therefore, the pattern of morphospace packing in the Campanian is relatively well supported. Note, however, that the TRiPS method has been criticized (Close et al. 2018), because it can only truly standardize diversity data once sampling is sufficient in each time bin (i.e., has reached the asymptote in the species discovery curve), and when sampling is not sufficient, it overestimates binomial sampling probabilities.

Different clades also show differing degrees of preservation: in particular, sauropods do not often preserve skull material. Is the steady decline in sauropods just a sampling artifact? This possibility cannot be excluded, and only increased sampling can resolve this. Nevertheless, the available fossil record of sauropods suggest broad-crowned forms were lost in the Late Cretaceous, which saw the radiation of the narrow-crowned titanosaurs (Chure et al. 2010). Therefore, as far as we can tell from the known fossil record, there is a real loss in sauropod morphological diversity in the Late Cretaceous, which fits with our disparity curve—a coupled loss in volume and distance-based disparity. Hadrosaurs and ceratopsians, on the other hand, are relatively well sampled in the Campanian–Maastrichtian, and the disparity trends recovered in this period for these clades are better supported.

Disparity Is Multifaceted

“Disparity” is sometimes discussed in terms of a uniform measure of morphological variety. The preferred metric in studies of the fossil record is distance-based metrics, as this measure is least sensitive to sample-size biases. Such biases are particularly a concern for the vertebrate fossil record. But our data show explicitly that disparity has two aspects, variance and volume, and different metrics are needed to capture these aspects (Foote 1999; Ciampaglio et al. 2001; Wills 2001; Kotrc and Knoll 2015a, b; Novack-Gottshall 2016a, b). Pairwise dissimilarity is a very robust measure of disparity, but it has been demonstrated that it “does not give an adequate estimation of the amount of morphospace occupied” (Ciampaglio et al. 2001: p. 711). Appreciating this can change the way we interpret disparity trends. The distance-based metric on its own suggests “disparity” was low in the Campanian. But was it low because the number of morphologies present was low, or because new clades were radiating into the same morphospace? This is an important distinction to make: the first scenario suggests decline and extinction, the other the diversification of successful clades.

The contrast between disparity trends in hadrosauroids and sauropods also demonstrates this. Both clades show approximately the same decrease in distance-based disparity from the Early to Late Cretaceous. Just considering this metric, we could conclude they were decreasing equally in “disparity.” When morphovolumes are measured, however, hadrosauroid disparity increases in volume in the Late Cretaceous, while sauropod disparity decreases. These distinct trends show that the two clades show radically different patterns: hadrosauroids flourished, while sauropods were in decline.

Our study explicitly underlines the importance of considering several aspects of disparity. The measures used here are not exhaustive but can be extended further to get a more in-depth understanding of disparity patterns (e.g., measures of clustering/dispersion in morphospace). Using modeling to identify particular modes of evolution in morphospace is also a promising approach and would be an interesting way to further our analysis (Novack-Gottshall 2016a, b). We show that measuring morphospace volume and packing, a relatively easy to implement extension of typical methods employed in disparity studies, can substantially increase the depth of the analysis.

There might be cases in which a measure robust to sample-size differences, such as mean pairwise distances or SoV, is the only feasible option, particularly where sample size is very small in some bins (as in our clade-specific analyses). A conservative approach is never to apply volume-based metrics, as measures might be so distorted by biases that no “real” signal will be visible, even on relatively good data sets. This is possible, but does not change the fact that distance-based metrics and volume-based metrics describe two different aspects of disparity. If we can only measure one, our inferences of the evolutionary processes underlying this trend must equally be limited.

Conclusion

Our results suggest that dinosaur ecological diversity, as far as it is reflected by dental and jaw disparity, was not declining in the Late Cretaceous, rejecting the hypothesis that dinosaur ecosystems were weakened preceding the extinction. Instead, herbivorous dinosaurs were flourishing in the Campanian, with no significant loss in disparity in the Maastrichtian. Clade-specific analyses reveal how different clades rose and fell in disparity during the Cretaceous, notably a long-term decline in sauropod disparity and morphospace packing in hadrosaurids, a group that was radiating in the Campanian. These trends reveal complex relationships among disparity, diversity, and diversity dynamics. Sauropods increased in diversity but decreased in speciation rates and disparity, and hadrosauroids, while diversifying, did not expand proportionally in morphospace, but rather retained a specific tooth/jaw morphology. We show how disparity trends are multifaceted and cannot be described in simple terms of increases and decreases, or by just a single metric. Thus, using a combination of metrics, each describing different aspects of disparity, is key to understanding the dynamics of morphological diversity through time.