Craniodental functional evolution in sauropodomorph dinosaurs

Abstract

Sauropodomorpha included the largest known terrestrial vertebrates and was the first dinosaur clade to achieve a global distribution. This success is associated with their early adoption of herbivory, and sauropod gigantism has been hypothesized to be a specialization for bulk feeding and obligate high-fiber herbivory. Here, we apply a combination of biomechanical character analysis and comparative phylogenetic methods with the aim of quantifying the evolutionary mechanics of the sauropodomorph feeding apparatus. We test for the role of convergence to common feeding function and divergence toward functional optima across sauropodomorph evolution, quantify the rate of evolution for functional characters, and test for coincident evolutionary rate shifts in craniodental functional characters and body mass. Results identify a functional shift toward increased cranial robustness, increased bite force, and the onset of static occlusion at the base of the Sauropoda, consistent with a shift toward bulk feeding. Trends toward similarity in functional characters are observed in Diplodocoidea and Titanosauriformes. However, diplodocids and titanosaurs retain significant craniodental functional differences, and evidence for convergent adoption of a common “adaptive zone” between them is weak. Modeling of craniodental character and body-mass evolution demonstrates that these functional shifts were not correlated with evolutionary rate shifts. Instead, a significant correlation between body mass and characters related to bite force and cranial robustness suggests a correlated-progression evolutionary mode, with positive-feedback loops between body mass and dietary specializations fueling sauropod gigantism.

Introduction

Sauropodomorph dinosaurs are represented by approximately 230 genera and include the largest terrestrial animals of all time (Mannion et al. 2011; Sander et al. 2011; Sander 2013). Hence, they are of great interest in understanding the evolution of gigantism and the biophysical constraints acting upon terrestrial life (Clauss 2011; Sander et al. 2011; Sander 2013). This has led to a concerted effort to comprehend the biological factors permitting sauropod gigantism (e.g., Klein et al. 2011; Sander 2013) and the macroevolutionary processes associated with it (e.g., O’Gorman and Hone 2012; Sookias et al. 2012; Benson et al. 2014).

Sauropodomorpha was the first dinosaur clade to become globally abundant, forming a dominant component of terrestrial faunas by the Norian (Galton and Upchurch 2004; Weishampel et al. 2004; Barrett et al. 2011) and remaining globally important until the end of the Cretaceous (Weishampel et al. 2004; Mannion et al. 2011; Sander 2013). This success has been attributed to the early adoption of herbivory by the clade (Barrett et al. 2011; Barrett 2014). Most “prosauropod” (non-sauropod sauropodomorph; Fig. 1) taxa are morphologically and, by assumption, functionally conservative in craniodental morphology (Galton 1985a,b, 1986; Crompton and Attridge 1986; Barrett and Upchurch 2007; Young and Larvan 2010), and some were probably at least facultatively omnivorous (Barrett 2000; Barrett and Upchurch 2007; Barrett et al. 2011). However, derived “prosauropod” taxa demonstrate a stepwise acquisition of craniodental characters hypothesized to be associated with increased body size, cranial robustness, bite force, and food gathering and processing potential that coalesce at the origin of Sauropoda (Barrett and Upchurch 2007; Upchurch et al. 2007; Rauhut et al. 2011; Sander 2013).

Figure 1 Simplified phylogeny of Sauropodomorpha, redrawn after Benson et al. (2014). Clades are numbered as follows: 1. Plateosauria; 2. Sauropodiformes; 3. Sauropoda; 4. Eusauropoda; 5. Neosauropoda; 6. Diplodocoidea; 7. Macronaria; 8. Titanosauriformes; 9. Titanosaura; 10. Lithostrotia.

The hypothesized functional shift toward greater cranial robustness and oral processing has been used to infer specialization toward obligate high-fiber herbivory and bulk feeding in the “prosauropod”–sauropod transition (Upchurch and Barrett 2000; Barrett and Upchurch 2007; Upchurch et al. 2007; Barrett et al. 2011; Barrett 2014). This ecological shift has in turn been hypothesized to have been a driver of gigantism (Hummel and Clauss 2011; O’Gorman and Hone 2012; Sander 2013; Barrett 2014; Benson et al. 2014), potentially with correlated progression or evolutionary cascades facilitating linkages between craniodental and postcranial characters associated with feeding and body mass (Barrett and Upchurch 2007; Sander 2013; Barrett 2014).

By contrast, eusauropods exhibit extensive craniodental morphological disparity (Young and Larvan 2010) and presumed functional disparity, with most taxa falling into two morphofunctional grades named after the distinctive tooth morphology of each (e.g., Fiorillo 1998; Barrett and Upchurch 2005; Chure et al. 2010). “Broad-crowned” taxa have robust skulls and more powerful jaw adductor musculature, whereas “narrow-crowned” taxa exhibit gracile skulls, abbreviated tooth rows, and weaker bites (Button et al. 2014). This variance in the feeding apparatus would have underpinned sauropod feeding ecology (Upchurch and Barrett 2000; Button et al. 2014), and tooth microwear studies have also demonstrated dietary divergence between sympatric “broad-crowned” and “narrow-crowned” taxa (Fiorillo 1991, 1998; Upchurch and Barrett 2000; Whitlock 2011).

Patterns of craniodental character evolution within Sauropodomorpha are complex, with convergent morphological trends observed between diplodocoids and titanosauriforms (Upchurch 1998, 1999; Barrett and Upchurch 2005; Coria and Salgado 2005; Chure et al. 2010; D’Emic et al. 2013; Button et al. 2014; MacLaren et al. 2017). These morphological convergences have been interpreted as functional and ecological convergences (D’Emic et al. 2013; Button et al. 2014) and are invoked in scenarios of sauropod macroevolution (e.g., Barrett and Upchurch 2005; Chure et al. 2010). In particular, the expansion of titanosaurs into “narrow-crowned” morphologies coincident with the decline of diplodocoid taxa has been taken as evidence of either an opportunistic ecological replacement (Coria and Salgado 2005; Chure et al. 2010; de Souza and Santucci 2014) or the possible competitive exclusion of rebbachisaurid diplodocoids by titanosaurs (Barrett and Upchurch 2005).

Morphological and functional disparity can be decoupled in biological systems (Lauder 1995; Alfaro et al. 2004, 2005; Wainwright and Richard 1995; Wainwright 2007). As a result, morphological diversity does not necessarily represent functional diversity, and further analysis is required to fully comprehend functional trends within the sauropodomorph feeding apparatus. The few existing biomechanical studies of sauropodomorph skulls have focused either upon individual taxa (Preuschoft and Witzel 2005; Witzel et al. 2011; Young et al. 2012) or a subset of taxa within Sauropoda (Button et al. 2014) or have been restricted to the mandible alone (MacLaren et al. 2017). The functional disparity of the feeding apparatus in non-sauropod sauropodomorphs has yet to be quantified, hampering evaluation of hypothesized ecological shifts, particularly toward the base of the Sauropodomorpha. Additionally, although macroevolutionary studies on sauropodomorph body mass have been conducted (Sookias et al. 2012; O’Gorman and Hone 2012; Benson et al. 2014; de Souza and Santucci 2014), there has been a lack of quantitative comparisons in considering the links between changes in feeding apparatus and gigantism.

Here, a series of functional properties are used to quantify the evolutionary mechanics of the sauropodomorph feeding apparatus using the largest data set compiled to date and including the first data from basal sauropodomorphs. We quantify how craniodental functional characters evolved on a temporal and phylogenetic scale and test the hypothesis that convergence in craniodental function plays a key role in sauropodomorph evolution. We calculate the rate of evolution of functional characters and, to test the hypothesis that evolutionary rate shifts in characters tied to the function of the sauropodomorph skull are linked to changes in body mass, we compare our data on trends in cranial function and mechanics to trends in body-mass evolution, using phylogenetic comparative methods. This allows, for the first time, quantitative testing of proposed convergence in feeding mechanisms across Sauropodomorpha and of the hypothesized linkages between innovations in the feeding apparatus and changes in body mass that have been frequently cited in models of the evolution of sauropodomorph gigantism (e.g., Barrett and Upchurch 2007; Upchurch et al. 2007; Sander 2013; Barrett 2014).

Materials and Methods

Multivariate Biomechanical Morphospace Analysis

Twenty-nine functional characters were measured from the crania, mandibles, and teeth of 67 sauropodomorph taxa (Table 1), using published images and firsthand observation of original material (see Supplementary Section S1). Similar studies (e.g., Anderson 2009; Anderson et al. 2011, 2013; Stubbs et al. 2013), including those incorporating sauropodomorphs (MacLaren et al. 2017), have often focused on the mandible alone, both to increase taxon coverage (Anderson et al. 2011, 2013) and to avoid potential compromise in the biomechanical signal due to the multiple functional roles of the cranium (Anderson 2009; Anderson et al. 2011, 2013; Stubbs et al. 2013; MacLaren et al. 2017). Characters from both the cranium and mandible were measured here to more fully capture feeding morphology and to increase taxon coverage. The resulting matrix of character scores was 79% complete. These characters are a combination of 20 continuous biomechanical metrics and nine binary characters that describe the functional properties of the feeding apparatus (see Supplementary Section S2 for character descriptions). Continuous measurements were standardized utilizing the z-transformation. These transformed data were then subjected to principal coordinate (PC) analysis, conducted in PAST (Version 2.17; Hammer et al. 2001), utilizing the Gower dissimilarity index (Gower 1971). Gower dissimilarity was chosen because it can be applied to mixed data sets containing both continuous and categorical data. The Mardia (1978) correction was applied to negative eigenvalues. Here, the total variance is taken as the sum of the magnitudes of all PC axes, which is then used to calculate the percentage of variance contributed by each positive axis. This method has the advantage of considering data from positive PC axes only, but it does mean that the total expressed variance will always be less than 100%.

Table 1 The taxa included in this study. The percentage of the biomechanical characters used in this study that could be measured for each taxon is indicated.

Character loading on the first two PC axes was gauged through linear correlations of each character with the PC axes and through the Spearman’s rank correlation coefficient, both calculated in PAST (see Supplementary Section S3), yielding a multivariate biomechanical morphospace. Taxa were grouped according to their distribution in the morphospace, yielding the following grades or clades for plotting: “basalmost Sauropodomorpha” (non-plateosaurian sauropodomorphs), “basal plateosaurians” (remaining non-sauropodiform sauropodomorphs), “basal Sauropodiformes” (non-sauropod sauropodiforms), “basal Sauropoda” (non-neosauropod sauropods and Camarasaurus), Rebbachisauridae, Dicraeosauridae, Diplodocidae, Brachiosauridae, Euhelopodidae, and Titanosauria.

Statistical comparison was made between these groups to test for significant differences in morphospace occupation. Comparisons were based upon scores on the first 31 PC axes (together accounting for 77.6% of the total variance) using a nonparametric multivariate analysis of variance (npMANOVA; Anderson 2001), with 100,000 permutations, conducted in PAST. The Bonferroni correction was applied to p-values for multiple comparisons.

To compare morphospace occupation through time, data were binned at the epoch level, resulting in six time bins (Late Triassic, Early Jurassic, Middle Jurassic, Late Jurassic, Early Cretaceous, and Late Cretaceous). Although the duration of these time bins is variable, these were chosen over finer-scale or equal stratigraphic divisions to avoid bin underpopulation. Sampling of sauropodomorph crania is generally poor but temporally variable: the smallest sample size (six) occurs in the Middle Jurassic, from which few sauropodomorph-bearing collections are known (Mannion et al. 2011; Upchurch et al. 2011). All other bins contain at least 10 taxa. The best-sampled bins are the Early Jurassic (17 taxa) and the Late Jurassic (16 taxa), of which the latter does seem to represent the acme of sauropodomorph diversity (Mannion et al. 2011; Tennant et al. 2016). Statistical comparison was made between each of these bins following an identical procedure to that described earlier. For those taxa with ranges encompassing multiple bins, the comparisons were conducted twice; once with all taxa assigned to the lower of the two time bins, and again with them in the upper bin.

Disparity Analyses

Scores from the first 29 PC axes were used to calculate the craniodental functional disparity of the above-listed groupings, to enable comparisons between non-sauropod sauropodomorphs and sauropods, and between the “prosauropod,” “broad-crowned” sauropod, and “narrow-crowned” sauropod functional grades (Supplementary Table S1). Taxa whose ranges extend through multiple bins were included in each bin. Disparity was calculated using the MDA add-on (Navarro 2003) for Matlab (MathWorks, Version 8.2.0.701, release R2013b). Bootstrapping, with 1000 replicates, was used to calculate 95% confidence intervals, which were used to gauge the significance of differences in disparity between these grades. In order to interrogate disparity patterns through time, disparity was calculated following the same protocol for each of the six time bins described earlier.

Four disparity metrics were computed for each of these analyses: the sum and product of variances and sum and product of ranges. The range-based metrics describe the total morphospace volume occupied by a particular group, whereas the variance-based metrics provide a measurement of mean dissimilarity between taxa. As measures of disparity, these metrics are relatively reliable at up to 25–30% missing data (Ciampaglio et al. 2001), compared with 21% missing data in the matrix analyzed herein. However, range-based measures are highly susceptible to sampling biases, whereas variance-based measures are robust to such biases (Butler et al. 2012). Hypervolume measurements as the products of eigenvalues are also more sensitive to values from axes of negligible variance (Wills et al. 1994; Wills 2001; Ciampaglio et al. 2001) and slightly more sensitive to character completeness (Ciampaglio et al. 2001) than summed metrics. For these reasons the sum of variances was the preferred metric for principal comparison; results for other metrics are presented in Supplementary Section S4.

To dissect these patterns further, partial disparity (Foote 1994) was calculated in each time bin for the following groups: the “prosauropod” and “broad-crowned sauropod” (including Euhelopus) functional grades, and Rebbachisauridae, Dicraeosauridae, Diplodocidae, Brachiosauridae, and Titanosauria. Although comparisons between nonmonophyletic groups can be problematic, the use of the first two groups is justified, as they represent distinct functional grades.

Phylogenetic Biomechanical Morphospace

To place these data in a phylogenetic context, an informal sauropodomorph supertree was constructed from published topologies (see Supplementary Section S5). This was time calibrated based on taxon occurrences at the stage level using the timePaleoPhy function in the ‘paleotree’ package (Bapst 2012) within R (Version 3.2.1; R Core Team 2013).

All taxa were dated to age level. In order to account for this uncertainty in dating, 1000 time-calibrated trees were produced, with random resolution of polytomies. This was performed in R, using the timePaleoPhy function in the ‘paleotree’ package (Bapst 2012) with the “minMax” argument. This draws a single-point date from a distribution drawn between minimum and maximum bounds. Zero-length branches were avoided through the “equal-sharing” method of Brusatte (2008), wherein time is shared equally with that of a preceding non-zero-length branch.

This time-calibrated tree was projected onto the first two PC axes (together accounting for approximately 52% of the total variance), using the phylomorphospace function in the ‘phytools’ R package (Revell 2012), with ancestral states for each node reconstructed under maximum likelihood. This yielded a phylogenetic biomechanical morphospace.

Mode of Craniodental Evolution

Taxon scores on the first two PC axes were mapped onto time-calibrated trees as continuous characters, allowing the mode of feeding apparatus evolution to be investigated through fitting models of trait evolution across dated trees. This modeling focused upon scores on the first two axes, as the functional variance along these axes is relatively well understood, reflecting robustness of the cranium and mandible, relative bite force, and tooth morphology. However, this approach does not account for the impact of functional variation linked to subsequent PC axes. To account for dating uncertainties, 1000 time-calibrated trees, with random resolution of polytomies, were generated following the procedure described earlier. Scores of taxa on the first two PC axis scores were then mapped onto these trees after pruning of taxa not included in the multivariate analysis. To compare the evolution of craniodental biomechanical characters with the results of Benson et al. (2014) on the mode of body-mass evolution, the following models were then fitted across these trees, using the fitContinuous function of the ‘Geiger’ R package (Harmon et al. 2008): Brownian motion, Ornstein-Uhlenbeck (OU), early-burst, Lambda, and Delta (for more information on these models, see Supplementary Section S6). Model selection was performed using the size-corrected Akaike information criterion (AICc).

Testing for Shifts in Evolutionary Rate of Multivariate Data

To test whether clades show a shift in the rate of evolution of craniodental characters as they move into regions of morphospace associated with each of the functional grades, multiple rate models were fit utilizing the transformPhylo.ML function of the ‘Motmot’ package (Thomas and Freckleton 2012), within R. The transformPhylo.ML function can be used to assign different clades different mean evolutionary rates. This permits testing hypotheses of rate shifts at or within a priori selected clades through comparison of the performance of two-rate models to a univariate null model of Brownian motion (Thomas and Freckleton 2012).

Two rate models permitting a rate shift (using the “clade” rate type of transformPhylo.ML) within the following clades were fit to 1000 time-calibrated trees: Sauropodiformes, Sauropoda (sensu McPhee et al. 2014), Diplodocoidea, and Titanosauria. In addition, two models featuring a rate shift concentrated at the branch leading to the Sauropodiformes and Sauropoda, respectively, utilizing the “branch” rate type, were fit to the 1000 dated trees to test for a marked rate shift associated with the onset of bulk herbivory (Barrett and Upchurch 2007; Barrett 2014).

Estimated body masses of 96 sauropodomorph taxa, taken primarily from the data set of Benson et al. (2014) with some additions (see Supplementary Section S7), were then mapped across these 1000 trees and pruned to feature only those taxa for which body-mass data are known. These analyses were then repeated for body-mass data to test whether any significant shifts in craniodental evolution at these nodes are correlated with shifts in the rate of body-mass evolution. This differs from the analysis of Benson et al. (2014), who identified the single most-likely point for a rate shift in body-mass evolution, as opposed to testing for the presence of a shift at a priori selected nodes.

Each two-rate model was compared with a null model of Brownian motion, fitted over the same 1000 time-calibrated trees. Model selection was performed using the AICc. The ΔAICc significance cutoff was calculated through a simulation study (Thomas and Freckleton 2012; Puttick et al. 2014). First, 1000 simulations of Brownian motion were mapped onto a tree. A two-rate model was applied to each and fitted utilizing the transformPhylo.ML function, employing the “tm2” algorithm to find the most likely single shift (Thomas and Freckleton 2012); the single-rate model described earlier was also applied. The 95^th percentile of the AICc improvement of the two-rate model, which equaled 8.49 for the craniodental PC scores and 9.79 for the body-mass data, was then used as the ΔAICc significance cutoff.

SURFACE Modeling of Multivariate Data

An elaboration of the OU process, allowing shifts between multiple local optima, can be used to model adaptive evolution in a Simpsonian landscape (Hansen 1997; Hansen et al. 2000; Blomberg et al. 2003; Butler and King 2004; Ingram and Mahler 2013; Mahler et al. 2013), with each local optimum value representing an adaptive peak (sensu Simpson 1944). Mahler et al. (2013) and Ingram and Mahler (2013) presented the SURFACE method, which uses an incremental approach to fit increasingly complex models featuring a greater number of local optima. This is performed using the AICc, with the sequential addition of a new local optimum to the model until doing so no longer results in a ΔAICc improvement exceeding a user-defined threshold (Ingram and Mahler 2013; Mahler et al. 2013). The resulting best-fit model defines the number of shifts to a new adaptive peak within the clade of interest, allowing the explanation of divergence in terms of distinct selective regimes.

Such a model can also be used to inspect instances of convergence within a clade. Observed morphological convergence between lineages may be the result of common selective pressures due to the adoption of similar ecologies but can also be the result of random drift (Stayton 2006; Ingram and Mahler 2013). SURFACE can be used to distinguish these scenarios in a second, “backward” iterative phase, which tests whether two lineages show shifts to occupy a common adaptive peak. This involves pairwise testing of the distinct selective regimes identified in the forward phase, fitting a model in which each pair is collapsed into a single selective regime (Ingram and Mahler 2013). Again, model selection is performed according to a user-defined ΔAICc cutoff. The result is to identify a new best-fit model that may contain independent convergence toward the same regime multiple times within the tree.

This procedure, implemented via the ‘surface’ package (Ingram and Mahler 2013) in R, was used to test whether the assembly of novel craniodental functional grades is associated with shifts to a new adaptive peak and to test for the presence of functional convergence between the diplodocoid and titanosaur lineages. Data for the first two PC axes (with only the first two axes used, as SURFACE is best used only on axes for which the biological significance is already well clarified; Ingram and Mahler 2013) were mapped onto 100 randomly resolved time-calibrated trees pruned of taxa lacking PC axis data. Additionally, the taxon Tornieria africana, whose position in the biomechanical morphospace was deemed a potential artifact of missing data (see “Discussion”), was found in preliminary analyses to be associated commonly with shifts that were ambiguous. Consequently, it was omitted from further analyses. A ΔAICc cutoff of 5 (the “conservative” threshold of Ingram and Mahler [2013]) was employed for the forward phase, with a cutoff of 0 used for the backward phase to yield a conservative estimate of the number of adaptive peaks (Ingram and Mahler 2013). Common shifts between sets of closely related nodes were compared between trees, and individual results were compared on the basis of AICc scores. Additionally, comparison was made to a single-shift OU model and a null model of Brownian motion (performed using the startingModel function in ‘surface’: Ingram and Mahler 2013) on the basis of AICc scores.

A second set of analyses was performed with the omission of taxa known from highly incomplete remains, with levels of missing data similar to Tornieria. These are detailed in Supplementary Section S8. SURFACE analysis was not conducted on the body-mass data, as the method is poorly suited to univariate data (Ingram and Mahler 2013).

Results

Multivariate Biomechanical Morphospace

PC axis variance scores are given in Table 2. PC axes 1 and 2 together account for ≈52% of the observed variance (Fig. 2A), PC3 accounts for ≈5%, after which variance scores rapidly tail off, so that PC axes above 11 account for <1% each (for other PC axes, character loading, and the positions of individual taxa, see Supplementary Section S3). PC1 is associated primarily with characters describing skull size, the relative robustness of the skull and jaw (average and maximum jaw depth, mandibular symphysis size and angle, external mandibular fenestra area, presence/absence of the lateral plate, snout breadth), relative bite force (anterior mechanical advantage, supratemporal and mandibular fenestra size, adductor muscle angle), degree of articular expansion, relative tooth row length, and tooth morphology and occlusal pattern. More positive values of PC1 refer to less robust skulls, lower anterior mechanical advantage of the mandible, and heterodont dentitions lacking tooth–tooth wear facets, whereas negative values describe larger and relatively more robust skulls and jaws with enlarged adductor chambers. Extreme negative values of PC1 describe narrow, precisely occluding dentitions. PC2 is similarly dominated by characters associated with relative cranial and mandibular strength (mandibular depth, development of lateral plates, snout morphology), relative bite force (anterior and posterior mechanical advantage, adductor fossa size), and tooth morphology. Positive values of PC2 are occupied by taxa exhibiting larger and relatively robust jaws, of more favorable mechanical advantage and with teeth occluding in an interdigitating manner, whereas negative values describe more gracile and less mechanically efficient jaws. The region of negative PC1 values and extremely negative PC2 values (<−0.15) is occupied by diplodocid taxa with highly procumbent dentitions.

Figure 2 A, Biomechanical morphospace plot of PC1 and PC2. Convex hulls illustrate the distribution of the three main functional grades: “prosauropod” (right), “broad-crowned” sauropods (top), and “narrow-crowned” sauropods (left). Taxa are plotted into the following grades: “basalmost Sauropodomorpha,” “basal plateosaurians,” “basal Sauropodiformes,” “basal Sauropoda”; and clades: Rebbachisauridae, Dicraeosauridae, Diplodocidae, Brachiosauridae, Euhelopodidae, and Titanosauria. Positions of exemplar taxa are labeled; skull illustrations accompany taxa listed in bold. Skulls are redrawn after the following sources: Eoraptor (Sereno et al. 2013); Plateosaurus (Button et al. 2016); Riojasaurus (Bonaparte and Pumares 1995); Shunosaurus (Chatterjee and Zheng 2002); Camarasaurus (Button et al. 2014); Brachiosaurus (Carpenter and Tidwell 1998); Rapetosaurus (Curry Rogers and Forster 2004); Nigersaurus (Sereno et al. 2007); and Kaatedocus (Tschopp and Mateus 2013). B, Phylogenetic biomechanical morphospace plot, featuring the supertree topology used herein projected onto the first two PC axes to illustrate functional trends. Strong shifts are observed between basalmost sauropodomorphs and plateosaurians; within non-sauropodiform sauropodomorphs, toward the base of Sauropoda; and in both Diplodocoidea and Titanosauriformes.

Table 2 Summary statistics of the first 31 PC axis scores.

Biomechanical Morphospace Group Separation

Global npMANOVA testing for group separation in morphospace yields significant results (p=9.99E-05). Pairwise comparisons indicate that “prosauropod” (non-sauropod sauropodomorph), “broad-crowned,” and “narrow-crowned” sauropods occupy significantly different regions of biomechanical morphospace (Table 3), and brachiosaurids can also be distinguished from each of these grades (Tables 2, 3). “Prosauropod” taxa are restricted to negative values of PC1 and strongly negative to weakly positive values of PC2, whereas “broad-crowned” sauropods occupy intermediate values of PC1 and greater values of PC2. Titanosaurs and diplodocoids cannot be distinguished from each other (Table 4), inhabiting a common “narrow-crowned” region of morphospace at negative values of PC1 and weakly positive to strongly negative values of PC2.

Table 3 Results of npMANOVA testing of multivariate biomechanical morphospace separation of the sauropodomorph functional grades, using scores on the first 31 PC axes. Bonferroni-corrected p-values for each pairwise comparison are given; significant results (<0.05) are highlighted in bold.

Table 4 Results of npMANOVA testing of multivariate biomechanical morphospace separation of the sauropodomorph grades and higher-level clades, using scores on the first 31 PC axes. Bonferroni-corrected p-values for each pairwise comparison are given; significant results (<0.05) are highlighted in bold.

However, if diplodocoids are subdivided to the family level (Table 5), then diplodocids can be distinguished from titanosaur taxa and from “broad-crowned” and “prosauropod” taxa. However, the separation between diplodocids alone and brachiosaurids is not significant. Rebbachisaurids and dicraeosaurids, however, show only significant separation from “broad-crowned” taxa, and so remain indistinguishable from all other “narrow-crowned” groups (Table 5). It should be noted that some of these results may reflect low statistical power associated with the low numbers of rebbachisaurids (n=3), dicraeosaurids (n=3), brachiosaurids (n=4), and diplodocids (n=5) used in this analysis.

Table 5 npMANOVA results of separation on the first 31 PC axes, with subdivision of the Diplodocoidea into separate clades. Bonferroni-corrected p-values for each pairwise comparison are given; significant results (<0.05) are highlighted in bold.

Results of statistical comparison of morphospace occupation through time show no significant difference between analyses utilizing the upper or lower ranges of taxa (Table 6). Global comparison for separation between time bins yields significant results (p=1E-05). Pairwise comparisons between time bins shows that there is no significant difference in overall biomechanical morphospace occupation between the Late Triassic and the Early Jurassic; between the Middle Jurassic and the Late Jurassic; and between the Late Jurassic, Early Cretaceous, and Late Cretaceous. All other pairwise comparisons were found to be significantly different.

Table 6 Results of npMANOVA test of significance of separation of multivariate biomechanical morphospace occupation through time. LT, Late Triassic; EJ, Early Jurassic; MJ, Middle Jurassic; LJ, Late Jurassic; EK, Early Cretaceous; LK, Late Cretaceous. Where uncertainty in dating led to taxon ranges crossing two time bins, two sets of pairwise comparisons were performed. For the first, these taxa were only included in the lower bin: Bonferroni-corrected p-values for these comparisons are given above in each cell. The second set of comparisons only included these taxa in the upper of the two time bins: results for these are given below, italicized and marked by an asterisk (*). Significant results are highlighted in bold.

Functional Disparity Comparison

Comparison of “Prosauropoda” to Sauropoda indicates that each exhibits comparable levels of overall functional disparity (Fig. 3A) regardless of the disparity metric used (see Supplementary Section S4). Comparison of the three sauropodomorph functional grades demonstrates a similar pattern. “Prosauropods” exhibit disparity values similar to, but slightly greater than, those of “narrow-crowned” sauropods, with “broad-crowned” forms the least functionally disparate (Fig. 3B). However, significant differences are only observed between “prosauropods” and “broad-crowned” sauropods under variance-based metrics (Supplementary Fig. S17).

Figure 3 Product of variances results for craniodental functional disparity, calculated from scores on the first 29 PC axes. Bars refer to 95% confidence intervals as calculated from bootstrapping with 1000 replicates. A, Comparison of total craniodental disparity of “prosauropod” and sauropod taxa. B, Comparison of total craniodental disparity exhibited by taxa of the “prosauropod,” “broad-crowned” sauropod, and “narrow-crowned” sauropod functional grades. Brachiosaurids were omitted from this analysis. C, Total sauropodomorph craniodental functional disparity through time. Results for other disparity metrics are given in Supplementary Section S4.

Total sauropodomorph craniodental functional disparity remains approximately constant from the Late Triassic to the Late Cretaceous (Fig. 3C). Although a decrease in disparity is observed in the Middle Jurassic for all plots, this result is significant in product of variances results only (Fig. 3C); no other significant differences are observed between any time bins for all disparity metrics used (Supplementary Fig. S18).

Dissecting these temporal trends through comparison of the partial disparity contributed by each group demonstrates that functional disparity in the first two time bins is dominated by “prosauropods,” but then rapidly drops off with the extinction of non-sauropod groups by the end of the Early Jurassic (Fig. 4A). During the Late Jurassic, non-neosauropods, diplodocoids, and macronarians exhibit similar levels of functional disparity. However, the expansion of the Macronaria into titanosaur morphotypes during the Cretaceous results in a continual increase in total observed macronarian disparity, at the expense of both basal sauropods (which are lost by the Late Cretaceous) and diplodocoids (which are extinct before the end of the Cretaceous) (Fig. 4A).

Figure 4 A, Partial craniodental functional disparity plotted through time for the following groups: “Prosauropoda” (non-sauropod sauropodomorphs), “broad-crowned” sauropods (non-neosauropods, Camarasaurus, and Euhelopus), rebbachisaurids, dicraeosaurids, diplodocids, brachiosaurids, and titanosaurs. Timescale given at the bottom. B, Scatter plot of body mass data used for this analysis through time. Body mass data were taken mostly from Benson et al. (2014), with some additions (see Supplementary Data), and plotted using the ‘Strap’ package (Bell and Lloyd 2015) in R. Points were plotted in the following groupings: “Prosauropoda: (non-sauropod sauropodomorphs), “Basal Sauropoda” (non-neosauropod sauropods), Diplodocoidea, and Macronaria.

Phylomorphospace

Superposition of sauropodomorph phylogeny onto the first two PC axes (Fig. 2B) demonstrates a strong trend within non-sauropod sauropodiforms toward decreased values of PC1 and increased values of PC2, culminating at the base of the Sauropoda. The other two most prominent functional shifts are observed in the seemingly convergent phylogenetic trajectories of Diplodocoidea and Titanosauriformes into negative values of PC1.

Modeling the Mode of PC Axis Score Evolution

Summary statistics of the fitted models are given in Table 7. For both PC1 and PC2 OU, delta, and early-burst models perform similarly to the Brownian motion model, in terms of both the average and total spread of AICc values. The average delta value for PC2 is >1, indicating that diversification in associated traits was more concentrated toward the distal parts of the tree. In contrast, the value for PC1 is ≈1, indicating a relatively constant rate of evolution through time. Lambda models perform better, and values of lambda approach 1 for both PC axes, indicating significant phylogenetic signal (as is also apparent from the phylomorphospace plot in Fig. 2B).

Table 7 Summary statistics of Brownian motion (BM), Ornstein-Uhlenbeck (OU), early-burst (EB), Lambda, and Delta models fitted to scores on PC axes 1 and 2, which were treated as continuous characters, across 1000 trees. The metric reported in each case refers to: BM, Brownian variance; OU, α; Lambda, λ; Delta, δ; EB, the rate parameter, r. ΔAICc values refer to the mean difference between AICc scores for each model and the null Brownian motion model (so that a positive ΔAICc value refers to a more negative AICc score than that for Brownian motion) across the 1000 trees.

Evolutionary Rates Analysis

Modeling of PC axis score evolution indicates that significant shifts in the rate of evolution of craniodental characters are seen only in a minority of trees within the a priori hypothesized clades (Table 8). The exception to this is Diplodocoidea, in which ≈68% of trees exhibit an increase in rate to approximately five times the background level. Evidence for a shift concentrated at the base of Sauropoda was weak.

Table 8 Summary of the rate-shift analyses calculated across 1000 dated trees. The mean maximum likelihood (ML) rate-shift estimation across all results is given, as is the percentage of trees in which a significant rate-shift signal (exceeding the simulated AICc threshold vs. a null single-shift model) is observed.

Significant shifts in the evolutionary rate of body mass were also sparse within any of the a priori tested clades, and evidence for a shift concentrated at the base of Sauropoda was weak (Table 8). The strongest evidence for a shift is located within the Neosauropoda in 60% of trees, and more specifically within Diplodocidea in 56% of trees.

SURFACE Results

SURFACE modeling of trees results in a mean improvement in AICc score of 107.14 versus a Brownian motion “starting model,” and of 104.74 versus a single-shift OU process (see Table 9 for summary statistics). Although ΔAICc scores are variable between trees, all show an improvement relative to OU or Brownian motion models. Results are variable between trees due to sensitivity to taxon dating, but common trends are apparent. All trees show at least one shift in local optimum between the base of Sauropodiformes and the base of Sauropoda (with many demonstrating multiple shifts), and 94% demonstrate a regime shift at the node Atlasaurus + Neosauropoda. Regime shifts are also commonly observed in Riojasaurus (79% of trees), Diplodocidae (61%), Camarasaurus (50%), and Titanosauria (43%) (Fig. 5). Results for the single best-performing SURFACE run (ΔAICc relative to single-α OU model=248.24, ΔAICc relative to Brownian motion model=304.29) are shown in Figures 6 and 7. Full results are given in Supplementary Section S8 and the supporting data. Sensitivity analyses returned similar results (see Supplementary Section S9).

Figure 5 Supertree topology with the position of shifts observed in >20% of trees indicated. Nodes associated with a shift are numbered, with branches in each regime colored. Numbers in boxes refer to the proportion of trees in which a shift at that point is observed. Numbered regime shifts are: 1. Riojasaurus; 2. Sauropodiformes; 3. Yunnanosaurus; 4. [Aardonyx+Sauropoda]; 5. Sauropoda; 6. [Atlasaurus+Neosauropoda]; 7. Diplodocoidea; 8. Diplodocidae; 9. Kaatedocus; 10. Macronaria; 11. Camarasaurus; 12. Titanosauria. For full results of SURFACE analysis, see Supplementary Data.

Figure 6 Best-performing tree under SURFACE analysis (ΔOU AICc=248.24). Evolutionary regimes are numbered; convergence to a common regime (here seen between Yunnanosaurus and Diplodocidae) indicated by reuse of the same number. Regime shifts are located at: 1. Sauropodomorpha; 2. Riojasaurus; 3. Yunnanosaurus; 4. Diplodocoidea [Aardonyx + Sauropoda]; 5. Sauropoda; 6. [Atlasaurus+Neosauropoda]; 7. Macronaria; 8. Titanosauria.

Figure 7 SURFACE results plotted onto PC1 and PC2. Large ovals refer to the positions of local optima in phylogenetic biomechanical morphospace. Regime 2 only contains the taxon Riojasaurus; hence it also represents a local optimum. Evolutionary regimes are numbered as in Fig. 6, with taxa colored according to regime.

Table 9 Summary statistics of SURFACE modeling across 100 dated trees. ΔBM AICc, improvement in AICc scores relative to a null Brownian motion “starting model”; ΔOU AICc, improvement relative to a single α OU process; total shifts, total number of shifts, with convergent shifts toward the same local optimum counted separately; conv. shifts, number of separate local optima resolved. The discrepancy between total shifts and conv. shifts is a measure of convergence; e.g., here three of the 11 identified shifts are convergent with respect to another regime, so only eight separate local optima are resolved. See Supplementary Material and online Supporting Data for the full results.

Multiple shifts toward convergent adaptive peaks within Sauropodomorpha are observed in 71% of trees, but this signal is dominated by individual genera. Average ΔAICc improvement versus a single optimum OU model for trees exhibiting convergence is 100.69 versus 854 for trees exhibiting no convergence; however, the number of convergent shifts does not seem to influence the comparative AICc improvement of different trees.

Camarasaurus demonstrates convergence with the group containing non-neosauropod sauropods in 47% of trees. The Diplodocidae or Diplodocinae show convergence with Riojasaurus in 19% of trees and with Yunnanosaurus in 19% of trees. As Riojasaurus and Yunnanosaurus often exhibit convergence (28% of trees), this results in at least some members of Diplodocidae converging with at least one of these taxa in 26% of trees. Additionally, Riojasaurus shows a convergent shift to the same local optimum as Atlasaurus, Rebbachisauridae, and Dicraeosauridae in a further 17% of trees. As a result, at least some members of the Diplodocoidea occupy a common regime to Riojasaurus in 36% of total trees. Meanwhile, diplodocoid and titanosaur taxa only occupy the same regime in 27% of trees. Apart from two instances of specific taxon convergence, this is as a result of common retention of the plesiomorphic neosauropod (or more inclusive taxon) regime, rather than convergent shifts. A full list of observed convergences across the 100 trees is given in Supplementary Section S8.

Discussion

Phylogenetic Biomechanical Morphospace and Functional Disparity

The three broad morphological grades observed in sauropodomorph skulls are functionally distinct, with each occupying significantly different regions of function space (Fig. 2; Table 3; see also Button et al. 2014; MacLaren et al. 2017). Phylogenetic signal in the data is strong (Table 7), but none of these grades are monophyletic (Fig. 2B). Whereas the variation observed herein can be characterized according to these grades it should be noted that sauropodomorph crania are poorly sampled: only around 30% of known sauropodomorph taxa exhibit sufficiently complete cranial remains to be included in these analyses. The rebbachisaurids, dicraeosaurids, and titanosaurs are particularly poorly sampled (Table 1); future discoveries of these taxa may result in an expansion of the overall recognized functional variety of Sauropodomorpha.

Non-sauropod (“prosauropod”) taxa possess nonoccluding dentitions and gracile skulls and mandibles. Although relatively conservative morphologically (Barrett and Upchurch 2007; Young and Larvan 2010) “prosauropods” are highly variable functionally, displaying a similar range of functional disparity to Sauropoda as a whole (Fig. 3A) and greater disparity than either of the two sauropod functional grades (Fig. 3B). However, this difference is only significant in comparison to “broad-crowned” taxa for variance-based metrics. In addition to the osteological characters considered herein, the purported presence of a rhamphotheca in some taxa (e.g., Bonaparte and Pumares 1995; Martínez 2009) may have further increased functional disparity, although the evidence for this structure remains contentious (Barrett and Upchurch 2007). This disparity in characters known to vary with diet in tetrapods (e.g., Reisz and Sues 2000; Sues 2000), such as the recurvature, serration, and imbrication of the teeth, degree of ventral jaw offset, and the efficiency versus speed of biting, probably reflects variance along the omnivory–herbivory spectrum. Plateosaurids, massospondylids, and some other taxa (here informally termed “core prosauropods”) cluster in a relatively restricted region of biomechanical morphospace. The high functional disparity within “prosauropods” is instead driven primarily by a relatively low number of eccentric taxa: basalmost sauropodomorphs, Riojasaurus, and non-sauropod sauropodiforms.

Basalmost Sauropodomorpha/Basal Saurischia

Panphagia, Eoraptor, and Pampadromaeus represent either basalmost sauropodomorphs or closely related basal theropods/basalmost saurischians (see discussions in Cabreira et al. [2011] and Martínez et al. [2011, 2013]). Whatever their relationships, they are informative with respect to the primitive condition for Sauropodomorpha and cluster in a unique region of the biomechanical morphospace. This divergence is due to traits commonly associated with faunivory—the possession of recurved teeth and very low anterior mechanical advantage of the jaw, which results in a weak, fast bite. Despite this, the dentitions of these taxa exhibit some features usually correlated with herbivory. The observed imbricate/subimbicrate tooth arrangement in Panphagia (Martínez and Alcober 2009) and Pampadromaeus (Cabreira et al. 2011) and the swollen crown bases and relatively coarse, obliquely inclined denticles of these taxa (Martínez and Alcober 2009; Cabreira et al. 2011; Sereno et al. 2013) and Saturnalia (Langer et al. 1999) more closely resemble those of other basal sauropodomorphs and ornithischians than the fine perpendicular serrations typical of carnivorous theropod teeth (Sereno et al. 2013). However, it is notable that these correlates are based on extant facultatively omnivorous iguanid lizards (Barrett 2000), rather than strictly herbivorous taxa. This, coupled with the variation along the tooth row in these taxa from anterior, conical, recurved “more carnivorous” teeth to posterior, less recurved, lanceolate, more coarsely denticulate “more herbivorous” teeth, and the mechanical characters identified herein, suggest that these taxa are probably best considered as omnivorous (Barrett et al. 2011).

“Core Prosauropods”

The majority of more derived “prosauropods” (Thecodontosaurus, Pantydraco, Arcusaurus, Efraasia, Plateosauridae, and Massospondylidae) cluster relatively tightly in a region of morphospace described by strongly positive PC1 scores and low positive and negative values in PC2. Relative to more basal taxa, they demonstrate the general absence of recurved teeth, an overlapping imbricate arrangement of the lanceolate maxillary and opposing dentary teeth (resulting in a more continuous cutting edge), and increased mechanical advantage of the jaw (which is increased in massospondylid vs. more basal plateosaurian taxa). Increased mechanical advantage is typical of herbivores (e.g., Janis 1995; Stayton 2006), where jaw closure speed is not important for prey capture.

These traits imply a shift to greater dependence on herbivory, corroborating evidence from both craniodental and postcranial morphology (e.g., Galton 1984 1985a,b, 1986; Galton and Upchurch 2004; Barrett and Upchurch 2007) and the relative abundance of “prosauropods” in Late Triassic and Early Jurassic faunas (Benton 1983; Galton 1984, 1985a,b, 1986; Barrett 2000). The absence of the cranial specializations seen in other herbivores, such as complex jaw actions, reflects an emphasis on gut processing in these taxa (Farlow 1987; Barrett 2014) rather than an absence of herbivory (contra Cooper 1981). However, other aspects of the craniodental apparatus, particularly the unspecialized dentition, still imply some degree of faunivory. All of these taxa maintain prominent heterodonty, with the premaxillary teeth and opposing dentary teeth being conical with reduced ornamentation and with more posterior teeth being lanceolate and highly denticulate. Tooth–tooth wear facets are absent, with the upper dentition lingually overlapping the lower dentition during jaw closure, indicating puncture-crushing of fodder by individual tooth crowns. The anterior teeth, as well as being used for plucking foliage, could equally be used for obtaining and dismembering animal matter, as seen in extant iguanids with a similar dentition (Barrett 2000). The biomechanical characters also reflect this. Anterior mechanical advantage is still low in these taxa, and they are highly variable with regard to characters like the degree of ventral offset of the jaw joint, a well-known herbivorous adaptation (Reisz and Sues 2000; Sues 2000) that both increases the in-lever of the jaw adductor musculature and permits a more even bite along the tooth row (Janis 1995). As a result, “core prosauropods” are best considered to have been relatively unspecialized herbivores, with some facultative omnivory likely (see also Barrett 2000, 2014; Barrett and Upchurch 2007).

Riojasaurus is strongly separated from the other “core prosauropods,” occupying more negative values of PC1 and lying closer to sauropod taxa in this axis (Fig. 2A). In addition to the osteological and dental characters considered herein, Riojasaurus has also previously been considered to have been distinct from other “prosauropod” taxa in the inferred presence of a premaxillary rhamphotheca (Bonaparte and Pumares 1995), for which the evidence is stronger in Riojasaurus than any other sauropodomorph (Barrett and Upchurch 2007). This functional separation, potential rhamphotheca, and its potentially atypical tooth morphology (Bonaparte and Pumares 1995; but see Barrett and Upchurch 2007) indicate that Riojasaurus would have been a more specialized herbivore than many other contemporary “prosauropod” taxa.

Basal Sauropodiformes

Sauropodiformes are highly functionally disparate (Fig. 2A). Jingshanosaurus and Yunnanosaurus, the most basal sauropodiform taxa considered herein, are both outliers in comparison to most “prosauropod” taxa. The dentition of Jingshanosaurus is relatively homodont, with many of the teeth being highly recurved (Zhang and Yang 1994: Fig. 9). Compared with phylogenetically proximate taxa (Massospondylidae, Yunnanosaurus) Jingshanosaurus exhibits a relatively low mechanical advantage and a gracile mandible with an elongated retroarticular process. These traits suggest a greater importance of faunivory in Jingshanosaurus. Yunnanosaurus differs primarily in its tooth morphology. It exhibits a homodont dentition of apicobasally elongate but mesiodistally narrow teeth, resulting in limited overlap between tooth crowns (Barrett et al. 2007). Moreover, denticles are absent, although incipient crenulations are apparent at tooth apices (Barrett et al. 2007). Although lacking well-developed wear facets resulting from tooth–tooth occlusion (Barrett et al. 2007) this dentition exhibits some convergences with sauropod teeth, suggesting more specialized herbivory in Yunnanosaurus.

Chuxiongosaurus, Anchisaurus, and Mussaurus are functionally similar to “core prosauropod” taxa, clustering tightly with them in PC1 and PC2. A prominent trend is observed through the more derived sauropodiforms toward more robust mandibles (Fig. 2B), greater mechanical advantage of the mandible, more procumbent dentitions, and reduced tooth rows. This culminates at the base of the Sauropoda with the onset of static, shearing occlusion between homodont teeth, which are clearly those of obligate herbivores (Calvo 1994; Christiansen 2000; Upchurch and Barrett 2000). As in the continuous metrics, this trend is also seen through the piecemeal acquisition of discrete characters, such as the lateral plates. The observed trends toward greater cranial and mandibular robusticity, increased bite force, and increased processing power of the dentition are all consistent with a shift toward bulk feeding on coarse, fibrous plant material, providing quantitative biomechanical evidence for the presence of such an ecological shift at the base of the Sauropoda and supporting previous suggestions based on morphology (Barrett and Upchurch 2007; Upchurch et al. 2007; Rauhut et al. 2011; Sander 2013; Barrett 2014).

Sauropoda

“Broad-crowned” sauropods (non-neosauropods, Camarasaurus, and Euhelopus) show relatively low disparity (Fig. 1A), remaining in a restricted region of the biomechanical morphospace associated with relatively high bite forces (due both to large adductor musculature and high mechanical advantage of the jaw), high cranial and mandibular robustness, and robust teeth meeting in interdigitating occlusion. This character combination would permit a broad diet, with cropping of tough or even woody material possible. These results are consistent with those obtained from previous biomechanical disparity analyses (Button et al. 2014; MacLaren et al. 2017) and from tooth microwear (Fiorillo 1991, 1998; Upchurch and Barrett 2000; Whitlock 2011) and myological reconstruction and finite-element modeling of sauropod skulls, which indicate the exertion and accommodation of greater bite forces in the more robust crania of “broad-crowned” rather than “narrow-crowned” sauropod taxa (Button et al. 2014).

Diplodocoids and titanosauriformes each show trajectories into overlapping regions of biomechanical morphospace associated with relatively gracile mandibles and reduced and inclined adductor chambers (Fig. 1A,B). These characters, coupled with the development of specialized cropping mechanisms, reduced length of the tooth rows, and slender teeth, suggest further emphasis on mesial cropping with little, if any, oral processing. Interestingly, diplodocoids and titanosaurs also exhibit convergence in characters not considered herein. All develop high tooth-replacement rates (Chure et al. 2010; D’Emic et al. 2013), probably in response to a highly abrasive/gritty diet (Chure et al. 2010; Whitlock 2011; D’Emic et al. 2013). Also, on the basis of both craniodental and postcranial evidence, rebbachisaurids (Barrett and Upchurch 2005; Sereno and Wilson 2005; Sereno et al. 2007; Whitlock 2011) and dicraeosaurids (Upchurch and Barrett 2000; Christian 2002; Barrett and Upchurch 2005; but see Whitlock 2011) are usually considered to have been specialized low browsers, and diplodocids (Upchurch and Barrett 2000; Christian 2002; Barrett and Upchurch 2005; Dzemski and Christian 2007; Christian and Dzemski 2011; Whitlock 2011) and probably some titanosaurs (Barrett and Upchurch 2005) would have indulged in low browsing for at least parts of the foraging cycle. However, diplodocids occupy a region of biomechanical morphospace separated from that of titanosaur taxa. The particularly low relative bite force of diplodocids, the loss of tooth–tooth wear facets, expanded articular, and the low-angle arrangement of the adductor muscles indicate a specialized feeding mechanism with an emphasis on propaliny, as has been previously suggested in Diplodocus (Barrett and Upchurch 1994; Upchurch and Barrett 2000; Young et al. 2012). The position of Tornieria as closer to the dicraeosaurids in PC2 is suspect, as it lacks data for the specialized characters loading in PC2 as described earlier. Kaatedocus occupies more negative values of PC2 and fewer positive values of PC1 than other diplodocids. Additionally, a shift within diplodocinae is detected within a minority of trees, resulting in Kaatedocus occupying its own independent regime (Fig. 5). This is a result of the highly elongate and gracile skull and mandible of Kaatedocus and its more rounded snout relative to adult specimens of other diplodocid taxa. Craniofacial remodeling, including squaring-off of the snout, has been suggested to have resulted in functional divergence between ontogenetic stages of Diplodocus (Whitlock et al. 2010). The separation of Kaatedocus from other diplodocid taxa seen here lends some additional support to this, although this is difficult to test in the absence of adult Kaatedocus cranial material.

Dicraeosaurids and rebbachisaurids show greater functional overlap with titanosaurs, with these groups following parallel functional trends toward a reduced tooth row, slender teeth, a reduction in the size of the supratemporal fenestra, and the development of a dentition with high-angled shear between the upper and lower tooth rows (Upchurch and Barrett 2000; Barrett and Upchurch 2005). This guillotine-like “precision-shear” bite would be suited to “nipping” off of foliage, with no further oral processing (Upchurch and Barrett 2000; Barrett and Upchurch 2005). Dicraeosaurids cluster most closely with the titanosaurs Antarctosaurus, Bonitasaura, and Brasilotitan (Supplementary Figure S15), all sharing a highly abbreviated tooth row and, at least with Bonitasaura, an extreme reduction of the supratemporal fenestra. These three titanosaur taxa further converge with rebbachisaurids in showing the common development of a highly squared-off snout.

Although neither rebbachisaurids nor dicraeosaurids show significant separation from titanosaurs in morphospace occupation, they cannot be distinguished from diplodocids, brachiosaurids, or “prosauropods” either. This reflects low statistical power due to the small number of rebbachisaurid and dicraeosaurid taxa sampled rather than functional overlap. Sampled dicraeosaurids exhibit greater inclination of the adductor chamber than titanosaurs, which would result in lower relative bite forces. However, the crania of dicraeosaurids remain poorly known (Table 1), limiting the functional inferences that can be drawn for these taxa. The same is generally true for rebbachisaurids (Table 1), with the exception of Nigersaurus. Nigersaurus is unique among sauropodomorphs in possessing teeth with asymmetrical enamel distribution (Sereno and Wilson 2005; Sereno et al. 2007; D’Emic et al. 2013) arranged in a complex, self-supporting and sharpening “tooth battery” (Sereno and Wilson 2005; Sereno et al. 2007) and in apparently having closed the supratemporal fenestra, necessitating a dramatic rearrangement of the adductor musculature (Sereno et al. 2007). The poor sampling of rebbachisaurid and dicraeosaurid crania means that the taxonomic distribution of specific functional complexes remains unknown, and the disparity present within these forms may be underestimated. Consequently, results for these clades, and functional comparisons with other taxa, must be treated with caution.

Titanosaurs are highly disparate (Fig. 2A), with basal taxa such as Malawisaurus and some more derived forms such as Rapetosaurus showing more functional similarity to brachiosaurids than to diplodocoids (see also MacLaren et al. 2017). This is not surprising, given the high diversity and ecological importance of titanosaurs during the Cretaceous (Barrett and Upchurch 2005; Upchurch and Barrett 2005; Mannion et al. 2011), especially in South America (Salgado and Bonaparte 2007; de Souza and Santucci 2014; García et al. 2015). The inference of a keratinous, rhamphotheca-like, cutting blade covering the posterior region of the dentary of Bonitasaura (Apesteguía 2004) might expand this disparity further. However, the existence of such a structure is considered unlikely here: Nemegtosaurus and Quaesitosaurus exhibit a morphologically similar coronoid to that of Bonitasaura (Wilson 2005) that is devoid of the dense neurovascular foramina whose presence was used to argue for a rhamphotheca in the latter taxon (Apesteguía 2004). Additionally, the position of such a blade behind the teeth would be of limited effectiveness.

Disparity Patterns through Time

The great functional breadth of basal taxa results in sauropodomorphs attaining maximum disparity early in their history, by the Early Jurassic (Fig. 3C), after which disparity plateaus. Indeed, with the exception of the product of variances results for the Middle Jurassic (see below), no significant differences in disparity are observed between any successive time bins (Fig. 3C; see also Table 6). Although much of the sauropodomorph fossil record is both patchily sampled and poorly dated (Mannion et al. 2011), particularly with respect to cranial remains, variance-based disparity metrics are generally relatively robust to sampling biases (Butler et al. 2012), and even the appearance of highly unusual sauropod taxa such as Nigersaurus (Sereno and Wilson 2005; Sereno et al. 2007) do not result in a significant expansion of total disparity. Craniodental remains too poorly known or too problematic to be included within our analyses (e.g., Saturnalia Langer et al., 1999, Ignavusaurus Knoll, 2010, Yimenosaurus Bai et al., 1990, Bellusaurus Dong, 1990, “Astrodon” [Carpenter and Tidwell 2005], Balochisaurus Malkani, 2003, Phuwiangosaurus [Suteethorn et al. 2009]) all conform to one of the functional grades identified within. This conformity further suggests that a genuine biological signal is represented. The apparent (although still largely nonsignificant) decline in the Middle Jurassic is a probable exception to this, as it coincides with the minimum sampled interval. The retrieval of a significant result here for the product of ranges results may be a reflection of the sensitivity of hypervolume calculations to axes of negligible variance (Wills et al. 1994; Ciampaglio et al. 2001), given that all positive axes were used to calculate disparity herein. The plateauing of total disparity is all the more surprising given the strong functional shifts observed within Sauropodomorpha during the Mesozoic, with significant differences in biomechanical morphospace occupation being observed accompanying the radiation of Sauropoda in the Middle Jurassic and of Neosauropoda in the Late Jurassic. This pattern of turnover, as opposed to expansion, through time is broadly similar to the pattern observed from biomechanical analysis of the mandible alone (MacLaren et al. 2017). However, mandibular disparity peaks in the Late Jurassic (MacLaren et al. 2017), differing from the early plateauing of whole-skull biomechanical disparity resolved here. This is partially due to the inclusion of cranial characters allowing increased sampling of titanosaurs herein. It may also reflect a decoupling of the disparity of the cranium and mandible within sauropodomorphs as a result of the greater number of the potentially conflicting roles performed by the cranium.

Dissection of these trends indicates that partial disparity (Fig. 4A) closely follows the diversity patterns of each group (see Barrett and Upchurch 2005; Upchurch and Barrett 2005). During the Early Jurassic, “prosauropods” represent >90% of observed craniodental disparity, but this rapidly tails off with the extinction of non-sauropods by the Middle Jurassic.

During the Late Jurassic, “broad-crowned” and diplodocoid sauropods show roughly equal levels of disparity, with brachiosaurids making a relatively minor contribution (Fig. 4A). However, through the Cretaceous, the expansion of titanosaur disparity occurs at the expense of other groups, through to the extinction of all other groups in the Late Cretaceous. This signal is consistent with either opportunistic or competitive interactions between titanosaurs and diplodocoids, but is too coarse to allow us to distinguish between these two hypotheses.

Functional Grades and Evolutionary Rate Shifts

Despite the strong functional separation between sauropodomorph craniodental grades, evidence for associated shifts in evolutionary rate are weak; of the a priori selected clades and particular branches tested here, only Diplodocoidea demonstrated significant shifts in both craniodental PC scores in a majority of trees (Table 8). Even then, significant shifts are detected in <70% of trees. Such a shift in Diplodocoidea is potentially attributable to sampling artifacts, with the first cranial remains for the clade represented by the simultaneous occurrence of multiple taxa in the Kimmeridgian–Tithonian Morrison Formation and Tendaguru Beds. However, the projected Bathonian age of the Diplodocoidea root used here compares favorably with the oldest known putative diplodocoid remains from the late Bathonian (Upchurch and Martin 2003; Upchurch et al. 2004; Mannion et al. 2011), providing some support for a genuine rate shift within Diplodocoidea.

Evidence for shifts in body-mass evolution associated with the observed functional shifts is similarly scant (Table 8). This corresponds with the results of Benson et al. (2014), who found a general pattern of reduction in rate over time, according to either an early burst or OU model. This pattern within Sauropodomorpha (and Dinosauria more generally) was attributed to the early filling of ecological niches (Benson et al. 2014). This would appear to corroborate the early plateauing of disparity observed herein. However, both early burst and single-shift OU models fit functional craniodental evolution relatively poorly (Table 7). They are strongly outperformed by models permitting multiple shifts in local optimum (Table 9), suggesting continued ecological innovation within sauropodomorph evolution.

Adaptive Landscape Shifts in Sauropodomorph Craniodental Evolution

Over macroevolutionary timescales, each local optimum recovered in SURFACE analysis can be interpreted as a peak in an adaptive landscape sensu Simpson (1944) (Ingram and Mahler 2013; Mahler et al. 2013). Shifts between adaptive zones will then be the result of circumvention of previous functional constraints through ecological or morphological innovation (Hunt and Carrano 2010).

However, interpretation of results is complicated by strong sensitivity to dating, particularly when inspecting potentially convergent shifts in relatively small clades. Nevertheless, consistent retrieval of the adoption of a novel adaptive zone proximate to the base of Sauropoda (Fig. 5) corroborates the results from the functional multivariate analysis indicating a significant ecological shift coincident with the evolution of sauropods.

A second consistent shift is observed at the node uniting Atlasaurus and Neosauropoda; this is a result of the sharp increase in disparity occurring at the basal split between diplodocoids and macronarians. There is also relatively good support for the adoption of independent local optima by Riojasaurus and Diplodocidae, each of which is notable for its functional dissimilarity to related taxa. The independent occurrence of some functional convergences between the skulls of Riojasaurus and diplodocids (although body-size disparity still suggests divergences in ecology) further suggests that specialized forms of herbivory had been adopted by sauropodomorphs as early as the Late Triassic.

Diplodocid cranial material is difficult to diagnose to fine taxonomic levels, with Tschopp et al. (2015) suggesting that the majority of previously referred cranial material cannot be assigned to genus level. However, this would make little difference to the results presented herein, as the only taxon found to occasionally differ from other diplodocids in functional measures—Kaatedocus—is associated with diagnosable postcranial remains (Tschopp and Mateus 2013; Tschopp et al. 2015).

Functional Convergence in Sauropodomorph Evolution

Although there is a strong signal for divergence, evidence for the convergent occupation of common adaptive zones by larger sauropodomorph clades is weak. In particular, despite the observed parallel functional trends between rebbachisaurids, dicraeosaurids, and titanosaurs (Fig. 2A), convergent shifts in these clades are not recovered here. Titanosaurs and diplodocoids only cluster around a common adaptive peak in a minority of trees (Fig. 5), and even in those cases it is typically the result of inheriting a plesiomorphic shift close to the base of Neosauropoda. Instead, Camarasaurus often shows convergence back toward the same regime occupied by neosauropods less derived than Atlasaurus (Fig. 5). It is possible that this might reflect a poorly reconstructed condition at the base of the Neosauropoda, with Camarasaurus inheriting the ancestral non-neosauropod state and titanosauriforms and diplodocoids then adopting the same regime convergently. This problem is exacerbated by the generally poor sampling of sauropodomorph crania; basal neosauropods are poorly represented, and the skull of Atlasaurus, the sister-taxon to Neosauropoda, is poorly known (with around 45% missing data). Nevertheless, the strong similarity observed here between Atlasaurus and Camarasaurus (Figs. 2, 7) would be expected to result in a conservative test of this pattern of regime shifts within the Neosauropoda. However, omission of Atlasaurus makes little difference to the results (see Supplementary Information Section S5). More problematically, the low numbers of dicraeosaurids and rebbachisaurids present in the analysis will make the detection of trends within each group difficult: it should be noted that removing taxa with <50% data would leave only one rebbachisaurid and one dicraeosaurid (Nigersaurus and Dicraeosaurus, respectively) in the analysis.

Still, this result would suggest that the parallel functional trends observed among rebbachisaurids, dicraeosaurids, and titanosaurs may result from random processes rather than selective forces due to occupation of a common adaptive peak (Ingram and Mahler 2013). This, coupled with the presence of unique functional complexes in taxa such as Nigersaurus, suggests that competitive replacement of diplodocoids by titanosaurs during the Cretaceous would have been unlikely. Interestingly, the titanosaur taxa showing the greatest levels of functional similarity to rebbachisaurids—Antarctosaurus, Bonitasaura, and Brasilotitan (Supplementary Fig. S15)—all postdate the extinction of diplodocoids in the Turonian. This corresponds with signals from body-mass evolution, with titanosaurs only expanding into the smaller-size niches occupied by rebbachisaurids and dicraeosaurids following the extinction of these diplodocoid groups (de Souza and Santucci 2014). This suggests opportunistic ecological expansion of titanosaur taxa following the extinction of diplodocoids, as opposed to competitive exclusion of dicraeosaurid and rebbachisaurid taxa during the Cretaceous (Coria and Salgado 2005; de Souza and Santucci 2014; contra Barrett and Upchurch 2005). However, the poor sampling of rebbachisaurid, dicraeosaurid, and, to a lesser extent, titanosaur crania make the elucidation of functional trends in these clades highly problematic. More complete dicraeosaurid and rebbachisaurid crania are required to adequately test for functional convergence between these taxa and titanosaurs.

Bulk Feeding and the Evolution of Sauropod Gigantism

The evolutionary modeling carried out herein demonstrates the adoption of a novel adaptive zone close to the base of Sauropoda. Multivariate functional analysis indicates that this manifests in characters associated with cranial robusticity, increased bite force, and increased food-processing potential, consistent with the hypothesized shift toward obligate herbivory and bulk feeding at the base of Sauropoda (Barrett and Upchurch 2007; Upchurch et al. 2007; Barrett et al. 2011; Sander 2013; Barrett 2014). It is also consistent with suggestions of ecological distinction between “prosauropods” and sauropods on the basis of postcranial anatomy (Barrett and Upchurch 2007; McPhee et al. 2015). This shift was unilateral: “core prosauropod” regions of biomechanical morphospace were not reoccupied by sauropod taxa (Fig. 2A), and total morphospace occupation remained significantly different from the Middle Jurassic onward after the extinction of “prosauropods” (Tables 3, 4). The relatively constant level of observed functional disparity (Fig. 3) appears to have been maintained despite episodic functional shifts, through turnover, with the vacation of previous niches. This signal is also seen in sauropod body mass; whereas “prosauropods” range over three orders of magnitude (from <10 kg to >2000 kg), even the earliest sauropods commonly reached masses of ≈5000 kg, and all were over 1000 kg (Fig. 4B), with the exception of the insular dwarves Europasaurus (Sander et al. 2006) and Magyarosaurus (Stein et al. 2010).

The observed functional shift within Sauropodiformes coincides with several trends in postcranial anatomy such as an increase in relative neck length (Barrett and Upchurch 2007) and the adoption of a graviportal quadrupedal stance (Barrett and Upchurch 2007; Yates and Kitching 2010; Yates et al. 2010; McPhee et al. 2014), although the evolution of the latter character was complex (McPhee et al. 2014). A shift toward bulk feeding has been suggested to have driven the evolution of these cranial and postcranial characters and body size through correlated progression (Barrett and Upchurch 2007; Upchurch et al. 2007; Barrett 2014), wherein integration of multiple functional complexes leads to correlated change through positive-feedback loops (e.g., Kemp 2006). Elongation of the neck and increasing body size would have expanded the feeding envelope (Barrett and Upchurch 2007; Sander 2013). Increasing body size would also result in reduction of mass-specific metabolic rate and expansion of gut capacity (Demment and van Soest 1985; Farlow 1987), enabling increased overall intake and fiber digestibility (Clauss and Hummel 2005; Clauss et al. 2013), if not increased total digestibility (Clauss et al. 2013). These would both be advantageous in bulk feeding on large quantities of low-quality, highly fibrous plant matter. The development of quadrupedalism would have facilitated all of these modifications of the postcranial skeleton (Barrett and Upchurch 2007).

The consistent trends observed in aspects of the sauropodomorph craniodental system observed herein suggest correlated progression; however, further investigation of the correlation of postcranial and cranial traits across sauropodomorph phylogeny is required to test this further. Correlated progression of these traits would explain the unilateral nature of this shift at the base of Sauropoda, with the apparent failure of sauropods to expand into small and medium-sized herbivorous niches following the extinction of “prosauropods.” The tight positive-feedback loops between multiple character complexes would be difficult to break (Kemp 2006; Barrett 2014), making reversions difficult, and some of these niches may have instead been filled by ornithischian and theropod herbivores in the Jurassic and Cretaceous. A basal specialization of sauropods toward bulk feeding may have then constrained them to large, megaherbivorous niches, promoting gigantism, as enabled by the unique character combination of the sauropod Bauplan (Sander et al. 2011; Sander 2013).

Such a scenario is consistent with the results from analysis of functional characters presented herein by “broad-crowned” taxa. However, the relationships between form, function, and diet are complicated (Lauder 1995; Wainwright 2007; Lautenschlager et al. 2016), and more direct dietary evidence, particularly from “prosauropod,” basal sauropod, and titanosaur taxa, are necessary to test such a scenario further.

Conclusions

Multivariate analysis of multiple functional traits of the sauropodomorph craniodental system indicates that variance can be characterized into multiple functional grades. Functional disparity remained approximately constant from the Late Triassic until the end of the Mesozoic, despite taxonomic turnover. Basal “prosauropod” taxa, despite being relatively morphologically conservative, are highly functionally disparate, possibly relating to variance along the omnivory–herbivory spectrum. A prominent shift in functional metrics of the feeding apparatus is observed starting within Sauropodiformes and culminating at the base of Sauropoda in the “broad-crowned” functional grade. Modeling of trait evolution indicates that this is associated with a significant shift into a new adaptive zone. This shift, toward greater cranial robusticity, increased bite forces, and increased processing potential, is consistent with the adoption of obligate high-fiber herbivory. This demonstration of an ecological shift provides quantitative evidence for models that posit the adoption of bulk feeding as a key factor in driving various character complexes, including sauropod gigantism through correlated progression.

Although Diplodocoidea and Titanosauria show multiple morphological and functional convergences in their skulls and teeth, each clade remains highly disparate. Diplodocids are found to be significantly different from titanosaur taxa, preventing characterization of a single “narrow-crowned” grade. Although titanosaur, dicraeosaurid, and rebbachisaurid taxa are functionally similar, evidence for competitive exclusion of diplodocoids by titanosaurs is weak. However, poor sampling of rebbachisaurid and dicraeosaurid lineages make the elucidation of trends within these clades difficult: discovery of more complete crania from these clades is imperative in order to test this hypothesis further.