Introduction

Theropod dinosaurs were one of the most remarkable lineages of terrestrial vertebrates in the Mesozoic Era. They attained a high level of taxonomic and ecological diversity (Weishampel et al. 2004), and represent the only dinosaur clade that survived the mass extinction event at the end of the Cretaceous, in the form of birds (Dingus and Rowe 1997). Mesozoic theropod species occupied the mass spectrum from a few hundred grams to more than six tonnes (Christiansen and Fariña 2004; Turner et al. 2007), and showed a huge diversity in skull morphologies (Fig. 1; Weishampel et al. 2004) and feeding strategies (Barrett 2005; Zanno and Makovicky 2011). Numerous papers have been published on the phylogenetic relationships of non-avian theropods and basal birds (e.g., Gauthier 1986; Sereno 1999; Clark et al. 2002a; Rauhut 2003; Smith et al. 2007; Choiniere et al. 2010). These analyses largely agree in the general interrelationships of major groups, but the phylogenetic position and validity of several clades (e.g., Alvarezsauridae, Ceratosauria, Compsognathidae, Therizinosauridae) and the detailed positions of many species are still controversial (see Rauhut 2003; Choiniere et al. 2010; Zanno 2010; Xu et al. 2011). In contrast to this rather high number of phylogenetic analyses, relatively few studies have investigated the morphofunctional evolution of theropod character complexes, or have addressed macroevolutionary questions, such as the importance of heterochrony or biomechanical constraints in theropod evolution. Those studies that have addressed such questions have overwhelmingly concentrated on the evolution of the limbs (e.g., Gatesy 1990; Wagner and Gauthier 1999; Middleton and Gatesy 2000; Hutchinson 2001a, b; Dececchi and Larsson 2011), growth patterns as indicated by bone histology (e.g., Erickson et al. 2001, 2009; Padian et al. 2001), body size, breathing and physiology (e.g., Schweitzer and Marshall 2001; Tuner et al. 2007; Benson et al. 2012) or, most recently, the evolution of feathers (e.g., Xu and Guo 2009) and the variety of diets (Barrett et al. 2011; Zanno and Makovicky 2011). However, H. Dromaeosaurid Velociraptor. I. Basal bird (avialian) Archaeopteryx. Modified from Rauhut (2003). apart from the works of Rayfield (2005, 2011), Barrett (2005), Barrett and Rayfield (2006), Sakamoto (2010), Zanno and Makovicky (2011), and Brusatte et al. (2012), the relationships between cranial diversity, functional constraints, diet, and evolutionary processes have received surprisingly little attention so far.

Fig. 1.

Diversity of skull shapes in theropod dinosaurs. A. Basal theropod Herrerasaurus. B. Coelophysid Syntarsus. C. Generalized spinosaurid. D. Abelisaurid Carnotaurus. E. Allosaurid Allosaurus. F. Ornithomimosaur Ornithomimus. G. Unnamed oviraptorid (originally referred to Oviraptor).

In recent years, geometric morphometrics has been used increasingly in palaeontology. The most comprehensive study focusing on the morphometrics of archosaurian skulls is Marugán-Lobón and Buscalioni (2003), though this investigation was based on simple distance measurements for three homologous units of the skull (braincase, orbit, and rostrum). Geometric morphometric studies have often been carried out for ornithischian dinosaurs (Chapman et al. 1981; Chapman 1990; Chapman and Brett-Surman 1990; Goodwin 1990; Dodson 1993), where they were mainly used for taxonomic purposes. Further studies deal with geometric morphometrics of the skulls of sauropods (Young and Larvan 2010) and of single or small numbers of taxa of theropods, such as Allosaurus and Tyrannosaurus (Chapman 1990), and Carnotaurus and Ceratosaurus (Mazzetta et al. 2000), as well as with isolated theropod teeth (D'Amore 2009). More comprehensive analyses were published by Marugán-Lobán and Buscalioni (2004, 2006) for extant birds. Recently, Brusatte et al. (2012) investigated the cranial diversity of nonavian theropods, using two data sets (small data set: 26 taxa and 24 landmarks; large data set: 36 taxa and 13 landmarks). These authors concluded that cranial shape of theropods is highly correlated with phylogeny, but only weakly with functional biting behaviour and thus that phylogeny was the major determinant of theropod skull shape. Their result challenges previous studies, which suggested marked functional constraints on the evolution of theropod skull shape (Henderson 2002; Rayfield 2005).

The goal of this study, which was initiated in parallel with that of Brusatte et al. (2012), is therefore to evaluate theropod cranial diversity and its relation to phylogeny, ecology and function, using geometric morphometrics. We used an independent subsample of taxa and a different combination of landmarks, as well as different proxies for cranial function. Thus, this study helps to test the results obtained by Brusatte et al. (2012), using a different set of data. Further, we investigated how shape variation of different skull regions is correlated with function, and how different dietary patterns affect skull shape.

Abbreviations.—AMS, average maximum stress; CVA, Canonical Variate Analyses; GPA, Generalized Procrustes Analyses; PC, principal component; PCA, Principal Component Analysis; PCO, Principal Coordinates Analysis; PIC, phylogenetic independent contrast; SSI, skull strength indicator; 2B-PLS, two-block partial least squares analysis.

Material and methods

Taxon sampling.—The cranium of 35 non-avian theropod species (+ four outgroup species and two Avialae) was analysed, using a two-dimensional geometric morphometric approach. The data set is based on published reconstructions of adult (or nearly adult) fossil material in lateral view (see supplementary information). For the majority of the ∼270 valid non-avian theropods described so far (see Butler et al. 2011), skull material is incomplete, juvenile or missing; therefore, the present data set includes only a small fraction (∼13%) of “real” theropod cranial diversity. However, usable reconstructions were available for at least one species from all major “family”-level clades of theropods. Thus, it is likely that the present data set successfully captures broad phylogenetic and functional patterns of cranial morphological diversity across theropods, even if it is not possible for us to document detailed variation within “family”-level clades.

Because the skulls of basal birds are extremely fragile, their preservation is usually poor. Furthermore, the vast majority of basal avialan taxa come from Konservat-Lagerstätten, such as the Solnhofen limestone of southern Germany or the famous Jehol beds of China, in which specimens are usually rather two-dimensionally preserved. Thus, good reconstructions of the skulls of basal avialan taxa are extremely rare. Because of the highly derived morphology of some taxa we were not able to place all landmarks on all specimens, so we created a second, smaller data set, including Pengornis and Shenquiornis. In addition to the paravian taxa from the larger data set, we included only a few nonparavian coelurosaurs (Compsognathus, Dilong, Ornithomimus, Erlikosaurus, Shuvuuia, and Conchoraptor) as outgroup taxa. We were able to include the alvarezsaurid Shuvuuia only in the small data set because of its bird-like skull shape, which includes the loss of the postorbital process of the jugal. All taxa are listed in SOM: table S1 (see Supplementary Online Material at  http://app.pan.pl/SOM/app58-Foth_Rauhut_SOM.pdf).

Geometric morphometrics.—Geometric morphometric approaches are commonly used to quantify and study interspecific or intraspecific shape variation across a number of specimens based on outline or landmark data that capture the shape of the specimens in question (Adams et al. 2004; Zelditch et al. 2004). The advantage of this approach is that landmark coordinates capture more information about shape than can be obtained from traditional morphometric measurements (linear distances, ratios and angular measures), which are often insufficient to capture the whole geometry of the original object (Hammer and Harper 2006). Geometric morphometric approaches have been reviewed and discussed by many authors including Bookstein (1991), Elewa (2004), Zelditch et al. (2004), and Adams et al. (2004).

In the large data set, cranial geometry was captured using 20 homologous landmarks, which were plotted on the reconstructed skulls using the program tpsDig2 (Rohlf 2005), which outputs a tps (thin plate spline) file with two-dimensional landmark coordinates and scale (size) data for each specimen. The chosen landmarks are of type 1 (good evidence for anatomical homology, such as points where two bone sutures meet) and type 2 (good evidence for geometric homology, such as points of maximal curvature or extremities) in the terminology of Bookstein (1991). The landmark data set includes the outer shape of the whole cranium (excluding nasal crests or horns on the skull roof), maxilla, antorbital fenestra, orbit, lateral temporal fenestra, jugal, quadratojugal, postorbital, and the posterior part of the skull roof (parietal and squamosal). For comparison, the data set shares eleven landmarks (55%) with the 26-taxon data set of Brusatte et al. (2012) and only five (25%) with the 36-taxon data set. In the smaller data set we used only 15 landmarks, owing to the fusion or loss of various skull elements in some taxa (Fig. 2; see SOM for the anatomical description of the landmarks).

The landmark coordinates were superimposed using Generalized Procrustes Analyses (GPA) in tpsRelW (Rohlf 2003) and PAST 2.09 (Hammer et al. 2001), which serves to minimize non-shape variation between species, such as that caused by size, location, orientation, and rotation (Zelditch et al. 2004). Afterwards, the Procrustes coordinates were subjected to Principal Component Analysis (PCA) using PAST and tpsRelW. This procedure assimilates data from all landmarks and reduces it to a set of PC scores that summarize the skull shape of each taxon and describe maximal shape variation in a morphospace (Hammer and Harper 2006).

Fig. 2.

Position of the landmarks on theropod skulls. A. 20 landmarks used for the large data set plotted on the skull of Ceratosaurus (modified from Sampson and Witmer 2007). B. 15 landmarks used for the small data set plotted on the skull of Anchiornis (modified from Hu et al. 2009). The description of landmark positions is given in the SOM at  http://app.pan.pl/SOM/app58-Foth_Rauhut_SOM.pdf.

Character evolution.–One main question is how far is theropod skull shape correlated with phylogeny? This is important for reconstructing skull shape changes using character-mapping approaches. If a strong phylogenetic signal is present, closely related taxa should occur closer to one another in morphospace than more distantly related taxa (Klingenberg and Gidaszewski 2010). The degree of this correlation can be calculated by mapping the phylogeny into the morphospace. This requires an ancestral state reconstruction of the morphometric data for each internal node on the tree using squared change parsimony (Maddison 1991; Klingenberg 2011). This algorithm collates the sum of squared changes of continuous characters (here PC scores) along all branches of the tree and calculates parsimonious ancestral states by minimizing the total sum of squared change across the phylogeny. The most optimal configuration of ancestral PC scores is mathematically described by a squared length value, which can be used to determine the strength of the correlation between shape and phylogeny (see below).

For this approach an informal supertree (sensu Butler and Goswami 2008) was created. The tree is based on several of the most recent phylogenetic analyses of theropods: basal theropods, including coelophysoids (Sues et al. 2011), Ceratosauria (sensu Rauhut 2003; Smith et al. 2007; Xu et al. 2009a), basal tetanurans (Benson et al. 2010), Coelurosauria (Hu et al. 2009), Tyrannosauroidea (Brusatte et al. 2010), and Dromaeosauridae (Csiki et al. 2010). In general, most trees show strong similarities with regard to the higher-level relationships of theropod dinosaurs, but may disagree on the positions of individual lineages. Euparkeria, Lesothosaurus, Massospondylus, and Plateosaurus were taken as outgroup taxa (SOM: fig. S3; the topology showing the interrelationships of the coelurosaurian taxa used in the small data set is shown in SOM: fig. S4). As discussed by Hunt and Carrano (2010), models of phenotypic evolution require information about time. To include this information, branch lengths were scaled to the present topology, using stratigraphic ages of taxa obtained from Weishampel et al. (2004) or from original literature (SOM: table S1).

In the next step, the original landmark data from the tps file and supertree were imported into MorphoJ (Klingenberg 2011). The landmark data were superimposed (GPA) and converted into a covariance matrix and subjected to Principal Component Analysis (PCA). Subsequently, the tree was mapped into the morphospace (see above). Furthermore, a permutation test was performed in MorphoJ in which the topology is held constant and the PC scores for each taxon are randomly permuted across the tree 10 000 times (Laurin 2004; Klingenberg and Gidaszewski 2010). For this approach both the tps file and the supertree were imported into MorphoJ. If the squared length of the supertree is less than occurs in at least 95% of the randomly generated trees then the phylogenetic signal may be deemed significant.

For a more detailed description of cranial shape changes through time, a similar character mapping approach was performed using the software package Mesquite 2.72 (Maddison and Maddison 2009). For this approach, a Nexus file containing the PC scores taken from PAST was produced. The data were used as continuous characters and mapped on the supertree using squared change parsimony (see above). The shape changes along the tree were visualised by plotting the ancestral state values in the morphospace within the visualization window of tpsRelW.

Shape versus function.—Another main question is whether cranial shape is correlated with skull function? As proxies for function we used the skull strength indicator (SSI; after Henderson 2002) and the average maximum stress (after Rayfield 2011) (SOM: table S2), which are both size-related parameters. Originally, Henderson (2002) calculated the skull strength at the longitudinal position of the orbital midpoint by treating the skull as a cantilevered beam, with the posterior region held immobile while a vertical force was applied at a point on the ventral edge of the snout. In the first step, we tested the correlation between SSI and shape using only those taxa that were also used by Henderson (2002). However, because SSI is strongly correlated with skull depth in the orbital region (SOM: fig. S1, table S2), this distance was measured in order to estimate SSI for all taxa included in this study. As a second proxy, we used the average maximum stress (AMS). The estimation of AMS is based on a two-dimensional finite element approach by calculating the value of maximum stress per element, which was carried out for the crania of various allosauroids and megalosauroids (Rayfield 2011). For those taxa, AMS shows a significant correlation with the skull length (see Rayfield 2011). We therefore used this relationship to estimate AMS for all taxa in this study (SOM: fig. S2, table S2).

The estimated values of SSI and AMS were logarithmically transformed to normalize for data distribution (Freckleton et al. 2002). To evaluate the correlation between shape and function we performed a two-block partial least squares analysis (2B-PLS; see Rohlf and Corti 2000) in MorphoJ using the Procrustes coordinates from the GPA (shape) and both functional proxies (functional coefficients). This method explores the pattern of covariation between two sets of variables by constructing pairs of variables that are linear combinations of the variables within each of the two data sets, and accounts for as much as possible of the covariation between the two original data sets. The strength of correlation is given by the RV coefficient and a p-value, generated by 10 000 permutations. Additionally, we divided the landmark data set into several modules (preorbital region, postorbital region, antorbital fenestra, orbit, and lateral temporal fenestra; see supplementary information), and reran the analyses. Using this approach, we were able to test how the shape variation of specific skull regions and skull openings is correlated with functional proxies. If different skull regions show different degrees of correlation with the functional proxies it is also possible that shape variations that occur in these regions are independent of each other. To test this we performed 2B-PLS for the preorbital and postorbital regions, as well as PLS in one configuration. The difference between both approaches is that the former is based on a separate Procrustes fit, testing the covariation between the shapes of the parts of each considered separately, whereas the latter is based on a joint Procrustes fit, testing the covariation between parts within the context of the structure as a whole (see Klingenberg 2009).

We performed phylogenetic independent contrast (PICs) analyses on the first three PC axes and both functional proxies including the small data set, which was based on the original data from Henderson (2002). First, the correlation of SSI and AMS with phylogeny was tested, by loading both proxies into Mesquite as continuous characters and mapping them on the supertree using squared change parsimony. Then, a permutation test was performed as outlined above. Assuming that a correlation of shape and function with phylogeny is present, and the terminal scores of both factors are non-independent, the scores have to be transformed into PICs (Felsenstein 1985). This was done using the PDTREE package for Mesquite (Midford et al. 2005). This procedure considers the relationships of species to each other and calculates contrasts that are statistically independent by assuming that character evolution can be modelled as a random walk (brownian motion model) and that characters change at a uniform rate per unit branch length across all branches. To produce “standardised” branch lengths, the original branch lengths were log_e-transformed prior to analysis, as recommended by Garland et al. (1992). To test if the data fulfil these assumptions the absolute values of the standardized PICs were plotted against: (i) their standard deviations (Garland et al. 1992), (ii) their estimated nodal values (ancestral PIC values), and (iii) the corrected age of their base nodes (Purvis and Rambaut 1995). Finally, the estimated nodal values were plotted against the corrected node ages. The assumption is justified if no significant correlation is present in all plots (Garland et al. 1992; Purvis and Rambaut 1995; Midford et al. 2005). Finally, if all quantities fulfil the four criteria, the resulting contrasts of PC scores were plotted against the contrasts of SSI and AMS. If any signal is present between shape and function then a statistical correlation should be detectable. This relationship was evaluated by calculating the coefficient of determination (R²) and the corresponding p-value.

Because both functional parameters were originally size related, we tested if the shape variation (after the landmarks were superimposed) is correlated with centroid size (log transformed) as a proxy for total size, which is defined as the square root of the sum of squared differences between landmark coordinates and centroid coordinates for any dimension. Because original size was previously removed from the data by performing GPA, a significant correlation between centroid size and shape variation could indicate an allometric trend (see e.g., Piras et al. 2011). The test was performed for all Procrustes coordinates and for the first three PC axes.

Shape versus ecology.—To test if feeding ecology is correlated with skull shape we categorized taxa based on dietary preference and feeding style, using the following characters: character 1, carnivorous vs. non-carnivorous; character 2, carnivorous vs. omnivorous vs. herbivorous; and character 3, weak biting vs. medium biting vs. strong biting. The subdivision of the second character is based on Barrett and Rayfield (2006), supplemented by data from Zanno and Makovicky (2011). In addition, Euparkeria, Daemonosaurus, and Zupaysaurus were coded as carnivorous, Eoraptor, Archaeopteryx, Anchiornis, Confuciusornis, and oviraptorids as omnivorous, and Limusaurus as herbivorous. The subdivision of the third character is based on Sakamoto (2010) (SOM: table S2). The characters were originally coded as discrete characters. To test the correlation between shape and diet, we performed 2B-PLS in MorphoJ. Scorings of characters 2 and 3 (see above) were transformed into covariates via a Principal Coordinates Analysis (PCA) with Euclidean distances and the transformation exponent c = 2. This analysis was repeated for the preorbital and postorbital data sets, to test whether dietary patterns have an influence on the shape variation of different skull regions, and for the small data set, which was mainly focused on skull variation in basal birds.

To test whether taxa with different dietary preferences (carnivorous vs. non-carnivorous and carnivory vs. omnivory vs. herbivory) occupied different regions within the morphospace, we additionally performed a NPMANOVA test (nonparametric multivariate analysis of variance) in PAST with 10 000 permutations, Euclidean distances and a Canonical Variate Analyses (CVA), using the first six PC axes (large data set), which describe over 75% of the total variance in theropod cranial shape (all outgroup taxa were excluded from the analysis, to avoid false-positive or negative signals). The NPMANOVA estimates whether the distribution of the three groups shows significant differences in morphospace (see Anderson 2001). One of the strengths of this approach is that it does not assume or require normality of the multivariate data. The test computes an F statistic and a p-value, pointing to a significant difference between dietary preferences if the F value is high and p-value less than 0.05. The p-values were Bonferroni corrected, which set the significance level lower than the overall significance to avoid false positive signals in a data set comparing more than two groups (Hammer and Harper 2006). In a second approach we excluded the oviraptorid taxa from the data set, to evaluate the influence of their aberrant skull morphology on previous results.

Shape versus phylogeny, function, and ecology.—Finally, as already addressed by Brusatte et al. (2012), we wanted to evaluate whether skull shape in theropods is better explained by function, ecology or phylogeny. To evaluate this, we calculated various agglomerative, hierarchical cluster topologies based on the Procrustes coordinates, the SSI, ASM, SSI, and ASM, as well as feeding ecology. In the latter case we combined the covariates from the PCO (see above) with SSI and ASM values, which therefore represent a diet-function cluster. The cluster analyses (UPGMA and Ward's method) were performed in PAST. In general, these kinds of cluster algorithms search for the two most similar objects and join them into a cluster. Then, the next most similar object is identified and joined and this continues until all objects are joined into one supercluster (Hammer and Harper 2006). In UPGMA, the clusters are joined based on the average distance between all members in the two groups, and in Ward's method the clusters are joined such that increase in withingroup variance is minimized (Hammer 2009). Subsequently, we loaded all topologies into MorphoJ, mapped them onto the tree, and performed a permutation test. By comparing the squared length and the p-value with that of the supertree topology, we were able to estimate which of the parameters best explained skull shape, based on a parsimony criterion. All topologies are shown in the supplementary information (SOM: figs. S5–S14).

Results

Morphospace and evolutionary trends in skull shape.— The whole data set (large data set) is summarized by 40 PC axes. Most shape variation in theropod skulls (including outgroup taxa) is captured by the first three PC axes (34.4%, 17.1% and 11.4% of total variance), describing over 60% of total variance in all. The first PC describes the relative length and depth of the snout (premaxilla and maxilla) and the anteroposterior dimensions of the antorbital fenestra and lateral temporal fenestra. The second PC describes the angle of the premaxillary body, the relative length and depth of the anterior extension of the jugal, the dorsoventral dimension of the antorbital fenestra, the anteroposterior dimension of the orbit, as indicated by the length of the suborbital body of the jugal, the depth of the postorbital region, and the position of the jaw joint. The third PC describes the relative length of the posterior extension of the maxilla, the inclination and depth of the lacrimal and the influence that this has on the relative size and shape of the antorbital fenestra and orbit, and the relative height of the quadratojugal. Most taxa plotted on the positive side of the first PC axis and are equally distributed along the second PC axis. Oviraptorids plot far from the other taxa because of high negative PC 1 scores, reflecting a shortened snout but enlarged lateral temporal fenestrae (Fig. 3).

Fig. 3.

Two-dimensional morphospace of the theropod skull shape (large data set) with distribution of the different theropod clades, based on the first two PC axes. The most derived taxa with respect to the first two components are visualized with relative warps. Shape changes along PC 1 reflect changes in the length and height of the snout and the anteroposterior dimensions of the lateral temporal fenestra. Changes along PC 2 reflect the height of the postorbital region, the size of the orbit and the position of the jaw joint.

The permutation test indicates that theropod cranial form is significantly correlated with phylogeny (tree length weighted by branch lengths = 0.78293449, p < 0.0001). Based on the ancestral state reconstruction, PC 1 remains almost constant from the hypothetical ancestor of basal archosauriforms to that of basal deinonychosaurs, indicating that this component was very uniform and close to the consensus shape (Fig. 4A). A significant deviation in this value is seen in the hypothetical ancestor of Neotheropoda, with a rapid shift to extreme snout elongation in the hypothetical ancestor of Coelophysidae. Similar deviations, and thus trends of snout elongation, can be seen in the hypothetical ancestors of Spinosauridae, Tyrannosauroidea, Ornithomimosauria, and Dromaeosauridae. Opposite deviations are present in the hypothetical ancestors of ceratosaurs and oviraptorids, which have markedly shorter snouts; this trend is especially pronounced in oviraptorids. In contrast, according to this component (PC 1) the skull shape of the hypothetical ancestor of basal birds does not differ greatly from that of basal deinonychosaurs (Fig. 4A). PC 2 (Fig. 4B) shows a continuous decrease from the hypothetical ancestor of Archosauriformes to that of Theropoda, thus demonstrating a trend of relative increase of the anteroposterior length of the orbit (in relation to orbit height) and decrease of the dorsoventral depth of the infratemporal fenestra. From the hypothetical ancestor of theropods to that of neotetanurans the value is relatively constant, but decreases again from this point to the hypothetical ancestor of deinonychosaurs. However, this value shows a number of deviations along the main phylogenetic axis of theropods. Whereas the hypothetical ancestors of Ceratosauria, Tyrannosauroidea, Ornithomimosauria, and Avialae show a further decrease of PC 2, values increase in the hypothetical ancestor of Abelisauridae, Allosauroidea, Megalosauroidea, and Tyrannosauridae, indicating a shift to a skull shape with a deep temporal fenestra, a dorsoventrally elongated orbit and a deep suborbital body of the jugal. In all of these taxa, these changes are related to a marked increase in body size (Fig. 4B). The third PC (Fig. 4C) shows a continuous increase from the hypothetical ancestor of archosauriforms to that of coelurosaurs, indicating a shift in the orientation of the lacrimal, resulting in an increase of the relative size of the antorbital fenestra and a decrease of the orbit. From the hypothetical ancestor of coelurosaurians to that of deinonychosaurs the component decreases only slightly. However, with the exceptions of the hypothetical ancestors of coelophysids, basal ceratosaurs and avialans, all other groups show a further increase in their respective lines. The skull shape of the hypothetical ancestor of Avialae fits with that of basal Deinonychosauria (Fig. 4C).

Fig. 4.

A—C. One-dimensional morphospaces of the most important theropod groups based on the ancestral state reconstruction of the first three principal components illustrating the major shape changes of the respective principal components. The X-axis represents the clade rank. D. Morphospace of the smaller data set showing the cranial diversity of paravians. Deinonychosaurs and basal coelurosaurs (with the exception of Conchoraptor) plot in a small area possessing skulls with a rectangular shape and large antorbital fenestrae, whereas basal birds occupy a large area possessing skulls with a more triangular shape, but with great variation in the beak shape and the relative size of the antorbital fenestra. Note that the symbols in D differ from that of A–C.

Table 1.

Results of the two-block partial least squares, showing the degree of correlation of Procrustes Coordinates with biomechanical coefficients (skull strength indicator and average maximum stress, both log transformed) and the diet (RV coefficient and p-value) with the whole skull and the preorbital and the postorbital regions. Abbreviations: AMS, average maximum stress; SSI, skull strength indicator.

The morphospace of the second, smaller data set is summarized with 13 PC axes. Here, the first three PC axes describe more than 70% of total variation (PC 1 = 34.6%, PC 2 = 21.9%, PC 3 = 14.1%). The paravian taxa are well separated from other coelurosaurian taxa, in which the close paravian outgroup taxa Velociraptor and Sinornithosaurus plot in the same morphospace as the more basal coelurosaurs Compsognathus, Dilong, and Ornithomimus. This is also true for the alvarezsaurid Shuvuuia. However, Erlikosaurus and Anchiornis lie closer to Archaeopteryx. Confuciusornis represents the greatest outlier within Avialae, whereas the skull shape of enantiornithine birds plots close to Archaeopteryx (Fig. 4D). Performing the permutation test for the smaller data set, the tree length is 0.45221156 (p = 0.0124), indicating a significant phylogenetic signal. The following trends of shape changes can be described within the avialan lineage: a decrease in the angle of the anterior margin of the premaxillary body accompanied by an elongation of the nasal process of the premaxilla, both leading to a markedly triangular skull shape, a decrease in the height of the maxillary body, a decrease in the size of the antorbital fenestra, a decrease in the depth of the jugal body, and a decrease in the size of the temporal fenestra.

Shape versus function.—The 2B-PLS reveals significant correlation between Procrustes coordinates (shape) and both functional parameters (SSI and AMS). Interestingly, this correlation is stronger in the postorbital than in the preorbital region (Table 1), which indicates that skull function has a stronger influence on shape variation in the postorbital region than in the preorbital region. Despite this difference, shape variation in both regions is significantly correlated (2B-PLS: RV coefficient= 0.261, p < 0.0001; PLS: RV coefficient = 0.494, p < 0.0001). Comparing the results of the shape variation for the three lateral skull openings, the orbit showed the strongest correlation with both functional proxies. The lateral temporal fenestra, as part of the postorbital region, also shows a significant correlation, whereas the shape of the antorbital fenestra does not. This result supports Henderson's (2002) proposal that orbit shape is an important indicator of the mechanical performance of a theropod skull (see below).

Based on the PIC analyses of the large data set, the skull strength indicator is significantly correlated with PC 2 and 3, whereas the average maximum stress is significantly correlated with PC 1 and 3 (Table 2). Thus, shape variation described by the first three PC axes seems to include morphofunctional information in respect to average maximum stress and skull strength. All three components at least partially describe the shape of the postorbital region (PC 1, anteroposterior dimension of lateral temporal fenestra; PC 2, depth of postorbital region and relative position of jaw joint; PC 3, relative height of quadratojugal), which, given the strong correlation between the shape variation in this region and functional indicators, might explain this correlation. Furthermore, PC 2 and 3 both concern the shape of the orbit and surrounding structures, supporting the result described above. In contrast, the PIC analysis based on the original data set from Henderson (2002) shows no significant correlation between shape and function (Table 2). However, this contradictory outcome could be the result of the small sample size of the original data set, and fails to be a diagnostic test for the PIC analysis (SOM: table S4).

The overall skull shape represented by the Procrustes coordinates is strongly correlated with centroid size in the large data set. Furthermore, all three PC axes show a significant correlation, with the second axis showing the best fit (SOM: table S3). This result differs from that of Brusatte et al. (2012), where only the second PC axis showed a significant relationship with centroid size. This indicates that all three PC axes still retain some information on size. A correlation between centroid size and Procrustes coordinates is also found for the small data set, thus also indicating a correlation between skull shape and skull size. However, in contrast to the results for the large data set, there is no significant correlation of centroid size with each of the first three PC axes.

Table 2.

Results of the PIC analyses, showing the degree of correlation between PC scores and biomechanical coefficients (R² and p-value). SSI* represents the small data set which includes only the original data from Henderson et al. (2002).

Table 3.

Results of the NPMANOVA (with and without oviraptorids) verifying the differences of the skull shape based on different diets (F value / p-value).

Table 4.

Results of the permutation test showing the correlation of the morphospace with the supertree and various cluster topologies (squared length and p-value).

Shape versus ecology.—The 2B-PLS analysis reveals a significant correlation between overall skull shape and dietary parameters. As was the case for the functional proxies, this correlation is stronger in the postorbital region (Table 1). A comparison between only carnivorous and non-carnivorous taxa shows a significant difference between both groups. When looking at non-carnivorous taxa in more detail (i.e., distinguishing between omnivorous and herbivorous forms), we find that, within the morphospace of the first three PC axes, herbivorous and omnivorous taxa overlap in a large area, but only slightly with the carnivorous taxa (Fig. 5A, B). This is also supported by the distribution of taxa in the CVA plot (Fig. 5C, D). A pairwise comparison of all three groups after performing a NPMANOVA test supports this observation, i.e. carnivorous taxa differ significantly from omnivorous and herbivorous taxa, whereas omnivores do not differ significantly from herbivores (Table 3). The exclusion of the highly aberrant oviraptorid taxa in both analyses does not affect these results. Performing the 2B-PLS analysis for the small data set recovers no significant signal, which might, however, be a result of the significantly smaller sample size.

Shape versus phylogeny, function, and ecology.—Based on the squared length, the topologies of both diet-function clusters are the most parsimonious explanation for the distribution of taxa within the morphospace, followed by the Ward's Cluster that combines the values of the average maximum stress and the skull strength indicator, and the supertree topology. The other functional clusters are less parsimonious than the supertree phylogeny (Table 4). Based on this result it seems that feeding ecology (as a combination of diet and biting performance) explains skull shape in theropod dinosaurs better than phylogeny. However, similar to the results of Brusatte et al. (2012), phylogeny seems to have a larger influence on skull shape variance than any single functional proxy.

Discussion

Major patterns in theropod cranial morphospace.—Analysing theropod cranial diversity via geometric morphometrics helps to quantify variation in shape. The data show that snout length and the length of the postorbital region are inversely correlated with each other (as already found by Marugán Lobón and Buscalioni 2003), whereas snout length is weakly positively correlated with orbit size (PC 1). The depth of the postorbital region is correlated with the relative position of the The theropod group with the most aberrant skull morphology is the Oviraptoridae, which shows an extreme negative PC 1, indicating a short and high snout and large postorbital region (see also Clark et al. 2002b; Osmólska et al. 2004). This conclusion is supported by the cluster analyses (SOM: figs. S5, S6) where all three oviraptorid taxa cluster together, but form a branch to the exclusion of all other theropods. Brusatte et al. (2012) also found oviraptorids to have the most unusual skull shape within theropods. Further outliers are the abelisaurid Carnotaurus (high positive PC 2), which has an extremely high skull with a small, dorsoventrally elongated orbit, and the ornithomimosaur Gallimimus (high negative PC 2), which has a flat skull with an enlarged orbit.

Fig. 5.

Two-dimensional morphospace of the theropod skull shape and CVA plot with distribution of the carnivorous (black triangle), omnivorous (grey squares) and herbivorous (white diamonds) taxa (non-theropod taxa are not shown). A. PC 1 vs. PC 2. B. PC 2 vs. PC 3. In both diagrams herbivorous and omnivorous taxa overlap one another strongly, whereas carnivorous taxa overlap only marginally with herbivorous forms, but moderately with omnivorous forms. C. CVA plot with oviraptorids. D. CVA plot with oviraptorids excluded.

Previous studies on recent birds demonstrated that the greatest morphological variation is also found in the rostrum (primarily in the shape of the premaxilla; Marugán-Lobón and Buscalioni 2006), which is correlated with a great variety of feeding strategies ( ; Smith 1993). In birds, the length of the rostrum is independent from that of the orbital and postorbital region, whereas the orientation of the rostrum has an influence on the shape of the posterior skull parts (Marugán-Lobón and Buscalioni 2006). However, this morphological variation seems to be controlled by only a small set of signal molecules (β-catenin, BMP, calmudulin, Dkk3, TGFβllr; Abzhanov et al. 2004, 2006; Wu et al. 2004, 2006; Mallarino et al. 2011). At this stage, it is of course speculative to hypothesize similar gene regulatory networks in theropods. However, by investigating the molecular control of the development of the rostrum shape in different bird and crocodile groups, it might be possible to create an extant phylogenetic bracket (sensu Witmer 1995) for cranial development in Archosauria. If it is possible to correlate different expression patterns of signal molecules with the morphological variation of the rostrum in birds and crocodiles in a mathematical way (see Johnson and O'Higgins 1996; Campás et al. 2010), it might be possible to use geometric morphometric methods to investigate the genetic control of cranial diversity in fossil archosaurs (see also Klingenberg 2010).

Shape versus phylogeny.—The skull shape of theropods is significantly correlated with higher-level phylogeny in both data sets. A similar result was found by Brusatte et al. (2012). However, the correlation in the smaller data set was weaker at a lower taxonomic level (see Jones and Goswami 2010), which could be the result of incomplete sampling.

A problem with this correlation is certainly the highly incomplete sampling of theropod phylogeny, because of incompleteness of the fossil record, most importantly the lack of good cranial material for most taxa. Thus, only 13% of known global theropod diversity could be included in the analyses. Cranial data from Megalosauridae (e.g., Eustreptospondylus, Torvosaurus), Alvarezsauroidea (included in the small data set), basal Oviraptorosauria (e.g., Caudipteryx, Incisivosaurus), basal Therizinosauroidea (e.g., Beipiaosaurus) and some basal birds are missing in the current analysis. Furthermore, many data are derived from skull reconstructions, the accuracy of which might be questionable. This problem is exemplified by the spinosaurid skull used, which is based on a combination of information on skull morphology from three different taxa (see SOM and Rauhut 2003 for details). Additionally, there are large time gaps, especially in the Early and Middle Jurassic as well as in the late Early Cretaceous to early Late Cretaceous, which might influence the results of the squared changed parsimony and PIC analyses. The time gaps mainly concern the basal radiation of the clades Averostra, Abelisauridae, Tyrannosauridae, Ornithomimosauria, Therizinosauridae, Oviraptoridae, and Dromaeosauridae.

Interestingly, the skull shape of basal birds does not differ significantly from that of basal deinonychosaurians, which are represented by Anchiornis. On the one hand this might mean that the skulls of basal birds and troodontids were very similar, from close relationship or similar diet preference. Alternatively, this observation might reflect phylogenetic uncertainty, and it is possible that Anchiornis was instead a basal member of Avialae, as originally described (see Xu et al. 2009b). This has to be tested in future phylogenetic analyses, as is the case with the currently published hypothesis that Anchiornis forms a clade with Archaeopteryx at the base of Deinonychosauria, outside of Avialae (Xu et al. 2011), in which case the similarities in shape between these two taxa might represent a low-level taxonomic signal. In contrast, the skull of the alvarezsaurid Shuvuuia plots within the coelurosaurs, close to Compsognathus, but far away from the avialan taxa. This is surprising, because the skull of Shuvuuia is extremely bird-like and possesses several characters otherwise unknown in non-avian theropods, which was one of the factors that led to the original identification of this animal as a basal avialian (see e.g., Chiappe et al. 2002). However, in contrast to basal birds, Shuvuuia possesses an enlarged antorbital fenestra and an extremely elongated maxillary body, which is similar to basal coelurosaurs. Thus, the different positions of Shuvuuia and basal birds within the morphospace are mainly based upon the shape of the snout. However, the shape analysis might lend further support to the hypothesis that alvarezsaurids are more basal coelurosaurs outside Paraves (Clark et al. 2002a; Novas and Pol 2002; Choiniere et al. 2010). Within Avialae, the skull shapes of Archaeopteryx, Confuciusornis, and Enantiornithes differ greatly from each other. Mapping the phylogeny onto the morphospace further demonstrates that the stem species of Avialae, Pygostylia, and Enantiornithes are well separated from each other, indicating a rapid diversification of skull morphology, probably in connection with the phylogenetic and ecological radiation of this group in the Early Cretaceous (Zhou and Zhang 2003; You et al. 2006; Li et al. 2010; O'Connor and Chiappe 2011). This result is further supported by the cluster analyses carried out for the large dataset, which includes only two avialan species. However, here both Archaeopteryx and Confuciusornis plot in very different positions in the morphological clusters, indicating strong dissimilarities in cranial shape (SOM: figs. S5, S6). Taking the skull morphology of the long-headed Zhongjianornis and Longipterygidae or the short-headed Sapeornithidae into account, the actual morphospace of basal birds is probably much greater than estimated.

Shape versus function.—As mentioned above, the study presented by Brusatte et al. (2012) attempted to correlate geometric morphometric data of the theropod cranium with biting performance, based on a mechanical advantage approach (see Sakamoto 2010). These authors found a significant, but weak correlation between shape and function and concluded that phylogeny was a more important determinant of skull shape than function. By comparing the RV coefficients of the 2B-PLS presented in both studies, the correlation between shape and function is about three times higher in the current data set, but the results of the PIC analyses are almost the same in both studies. The average maximum stress explains only about 6.6% and the skull strength indicator about 9.5% of total cranial shape variation in theropods. However, an interesting result of the present study is that the shape of the postorbital region is better correlated with function than that of the preorbital region, and that the shape of the orbit shows the strongest correlation, whereas the shape of the antorbital fenestra shows no significant correlation. These results are also supported by the correlation of the skull strength indicator with the second and third PC axes (see PIC analysis), which include aspects that describe the shape and depth of the orbital and postorbital regions (e.g., shape and size of the orbit, depth of the suborbital body of the jugal, the position of the jaw joint). These results support previous morphofunctional studies that demonstrated that the orbital and postorbital regions of theropod skulls seem to be more important for an understanding of skull biomechanics than the snout. According to Henderson (2002) there is an inverse relationship between the size of the theropod orbit and the resistance of the skull to bending in the sagittal plane, and the narrowness of the orbit shows a positive correlation with skull strength. Henderson (2002) interpreted the correspondence between orbit size and shape to relate to increases in the amount of bone comprising the skull, and so its strength. He furthermore proposed that the shape of the orbit in theropods with strong skulls is governed by the requirements of the posterior half of the skull to resist muscle-generated forces associated with prey capture and/or dismemberment. Furthermore, based on a finite element approach, Rayfield (2005) demonstrated that high tensile and shear stresses especially affect the jugal bone, the shape of which in turn influences the anteroposterior dimension of the orbit. According to Rayfield (2011), the postorbital part of the skull (especially the squamosal, quadrate and quadratojugal) experiences most of the Von Mises stress associated with biting in theropods. It is thus to be expected that these stresses result in stronger mechanical constraints on the postorbital region of the skull than on the preorbital region, which is confirmed by the findings presented here. Furthermore, since overall stress acting on the skull (for which the functional data used here are a proxy) is necessarily correlated with skull size, it is not surprising that the shape of the postorbital region is also correlated with centroid size (as proxy for skull size). This is also in accordance with the finding of Rauhut (2007) that phylogenetic characters in the braincases of theropods seem to be influenced by size. Thus, the results of biomechanical studies (e.g., Rayfield 2011), phylogenetic data (Rauhut 2007) and morphometric analyses (present study) converge on the result that the postorbital region of the skull is particularly strongly influenced by biomechanical forces and related aspects of size. It might be worth noting, however, that the correlation between functional proxies and the preorbital region, found in both the 2B-PLS and in the PIC analysis (average maximum stress vs. PC 1), indicates that snout shape is also functionally constrained, but to a weaker degree.

The lack of correlation between the shape of the antorbital fenestra and the proxies for cranial mechanics used here might be surprising at first glance. Witmer (1997) suggested that the size and shape of the antorbital fenestra were largely determined by biomechanical factors, in that bone was resorbed opportunistically by pneumatic diverticula if it was not biomechanically necessary. Thus, our results seem to contradict this hypothesis. However, it is necessary to keep in mind that the shape of the antorbital fenestra is only captured by three landmarks and that the proxies for cranial function used here are proxies for the overall stress acting on the skull, not for the stress distribution within the cranium. Thus, although the overall stress experienced by the skull during biting might be high, this does not contradict the idea that particular parts of the cranium experience stresses that are below the threshold for the formation of bone, as has been demonstrated by Finite Element Structural Synthesis for other dinosaurs (Witzel and Preuschoft 2005; Witzel et al. 2011). In this respect, the hypothesis that the antorbital fenestra is associated with a region of low stress might even be supported by the weaker correlation between the shape of the preorbital region with the functional proxies for overall stress in the cranium, since the latter indicates that there might be less overall stress acting on this part of the skull (which, in turn, is in accordance with the results obtained by Rayfield [2011]).

It must be emphasised that the current approach is much more simple than that used by Brusatte et al. (2012), and both functional parameters used here are necessarily strongly correlated with size, since larger taxa are expected to experience higher total stresses than smaller taxa (see Henderson 2002; Rayfield 2011). In contrast, Brusatte et al. (2012) used a metric of functional biting profiles, which are more complex and independent from size. However, the impact of both the functional parameters (AMS and SSI) on skull shape is rather small, as previously hypothesized by Brusatte et al. (2012). Nevertheless, it would be worthwhile to test the results further by using other functional parameters or metrics (e.g., finite element analyses) or by using a different landmark combination, different skull views, or a three-dimensional approach (see also Brusatte et al. 2012).

Shape versus ecology.—As is the case with phylogeny and function, dietary patterns are also correlated with skull shape variation. Interestingly, we found a higher correlation for the postorbital region than for the preorbital region. This might be regarded as a surprising result, because one might expect that dietary preferences would have an equal or even stronger influence on the shape of the snout, the main organ used to obtain and process food. However, the result probably reflects the generalised subdivision in diet preference (carnivory vs. omnivory vs. herbivory) and that the dietary patterns also contain information on biting behaviour (see character 3), and thus biomechanical aspects (see above). Furthermore, different tooth morphologies, which are not captured in the current approach, might better reflect more specific dietary patterns than overall shape of snout (see also Smith 1993; Barrett et al. 2011; Zanno and Makovicky 2011).

Carnivorous, omnivorous, and herbivorous theropods occupy large areas of morphospace. Based on disparity analyses, Brusatte et al. (2012) demonstrated that non-carnivorous theropods (i.e., herbivores and omnivores) display a higher cranial disparity than carnivores. This result was largely influenced by the aberrant skull shape of oviraptorid dinosaurs, and their exclusion significantly reduced the difference in disparity between carnivorous and non-carnivorous forms. However, the current results indicate that both omnivorous and herbivorous taxa overlap strongly in morphospace, but only slightly with carnivorous theropod taxa (Fig. 5, Table 3). The exclusion of oviraptorids does not affect this result. This might indicate that changes between omnivory and herbivory had only small effects on skull shape, whereas changes between carnivory and non-carnivory (i.e., omnivory or herbivory) affect skull shape more strongly. Alternatively, this result might also indicate that the classification of omnivorous and herbivorous taxa on the basis of skull shape in fossil taxa is rather poor, particularly because among extant vertebrates the boundary between herbivory and omnivory is gradual (Zanno and Makovicky 2011). Further information should be incorporated, such as pelvis morphology, fossilized gut contents and coprolites, or the presence of gastroliths (Zhou et al. 2004; Barrett and Rayfield 2006; Zanno and Makovicky 2011), i.e., evidence that is independent from cranial morphology. Subdividing the taxa into carnivorous and non-carnivorous forms leads to a significant difference as well, supporting the assumption that a change between carnivory and non-carnivory had a marked effect on cranial shape. However, the large area of morphospace occupied by non-carnivorous taxa indicates further that an omnivorous or herbivorous diet was an important resource for many small-bodied theropods and may have been a fundamental driver of theropod evolution, especially within the coelurosaurian (Zanno and Makovicky 2011; Brusatte et al. 2012) and avialan clades (O'Connor and Chiappe 2011).

Conclusions

This study investigates the cranial shape variation of various non-avian theropod dinosaurs and some basal birds in a broad scale approach, using two-dimensional geometric morphometrics. The results indicate that most variation in the theropod cranium occurs in the shape of the snout (PC 1), the shape and size of the orbit and the shape of the postorbital region (PC 2 and PC 3). Interestingly, especially in the first principal component, there is surprisingly little change in the ancestral node reconstructions from basal archosauriforms all the way to birds. This might indicate that, in respect to snout shape, there is a generalist archosauriform condition, from which different clades deviate when specializing for certain ecological niches. Oviraptorids had the most aberrant skull shape in the theropod data set, but we further found that the skull shapes of abelisaurids and spinosaurids, for instance, differ greatly from those of other large bodied predators, with the former being characterized by an unusually deep and short skull and the latter by the other extreme, a low and long skull. Skull shape is strongly correlated with phylogeny, but also feeding ecology. Interestingly, the skull shape of non-carnivorous taxa differs significantly from that of carnivorous forms, which might indicate that a change between both diet preferences strongly affected skull shape. In sum, non-carnivorous taxa occupy large areas of morphospace, indicating that a diverse diet might have been a fundamental driver of the evolution of the morphological diversity of theropod skulls. We further found that skull shape is also correlated with dietary patterns, average maximum stress and skull strength, indicating that the cranium of theropods (especially the shapes of the orbital and postorbital regions) was constrained by ecology and function, especially biomechanics.

Using a different subsample of taxa (including some outgroup taxa and basal birds) and landmarks, we were able to support most of the results found by Brusatte et al. (2012) in their independently conducted study of cranial shape in theropods. These include major shape changes along the first two PC axes, i.e., the relative length of pre- and postorbital regions (PC 1), and the depth of the postorbital region and the size and shape of the orbit (PC 2), leading to a very similar distribution of taxa within the morphospace. We also found that skull shape is significantly correlated with phylogeny, but also with function. Comparing the squared lengths of the permutation tests, the functional clusters of the average maximum stress and the skull strength indicator taken separately are less parsimonious, which might indicate that phylogeny has a stronger influence on skull shape than function, as already hypothesized by Brusatte et al. (2012). However, a single Ward cluster, which includes both functional proxies, was found to be more parsimonious, challenging the previous result, whereas the best match was found for those clusters that include dietary patterns and both functional proxies. Nevertheless, given the consensus in several important points between Brusatte et al. (2012) and the current study, the other results mentioned above can be considered as strongly supported.

Based on the current results, we prefer not to speculate about which factor had the largest influence on theropod skull shape, as this might vary within different groups and also depends on the different proxies used. Further investigations into the relationship between function and cranial shape are needed, preferably using more specific data not only on overall stress, but also on stress distribution within the cranium, such as data from finite element analyses. To further test the current results it would also be worthwhile to use other subsamples of landmarks, incorporating data from other views of the skull, broadening the taxonomic sample, and, as far as possible, using three-dimensional data. Furthermore, it would also be interesting to evaluate the variation of skull shape of one species based on different reconstructions or different specimens, or differences between closely related species, and how this might influence the results of the PCA.