Structure-function covariation with nonfeeding ecological variables influences evolution of feeding specialization in Carnivora

See allHide authors and affiliations

Science Advances  07 Feb 2018:
Vol. 4, no. 2, eaao5441
DOI: 10.1126/sciadv.aao5441


Skull shape convergence is pervasive among vertebrates. Although this is frequently inferred to indicate similar functional underpinnings, neither the specific structure-function linkages nor the selective environments in which the supposed functional adaptations arose are commonly identified and tested. We demonstrate that nonfeeding factors relating to sexual maturity and precipitation-related arboreality also can generate structure-function relationships in the skulls of carnivorans (dogs, cats, seals, and relatives) through covariation with masticatory performance. We estimated measures of masticatory performance related to ecological variables that covary with cranial shape in the mammalian order Carnivora, integrating geometric morphometrics and finite element analyses. Even after accounting for phylogenetic autocorrelation, cranial shapes are significantly correlated to both feeding and nonfeeding ecological variables, and covariation with both variable types generated significant masticatory performance gradients. This suggests that mechanisms of obligate shape covariation with nonfeeding variables can produce performance changes resembling those arising from feeding adaptations in Carnivora.


Evolutionary convergence has been extensively documented across an array of organismal clades and on multiple levels of biological organization (1). Morphological or phenotypic convergence underlies many hypotheses of organismal adaptation in the biological sciences, because selection acts directly on the performance of phenotypes through time. The repeated, or iterative, evolution of a limited number of feeding morphotypes [also called feeding ecomorphologies (2)] is a well-demonstrated pattern in many groups of vertebrates, including throughout the more than 60-million-year evolutionary history of Carnivoramorpha (living and extinct Carnivora and their stem relatives) (3). Such a pattern of repeated craniodental shape convergence toward a few dominant ecomorphologies associated with feeding specializations is hypothesized to have significantly influenced extinction risk as the clade diversified during the Cenozoic Era, making this clade an excellent model system for studying diversity- and environment-mediated macroevolution (37). A well-demonstrated process, termed the “macroevolutionary ratchet,” may explain irreversible convergent trends toward increasingly specialized feeding adaptations that, in turn, reduce the evolutionary flexibility of those specialized species and lineages in the face of subsequent environmental changes (for example, sabertoothed predators) (3, 7). The major force thought to drive such trends is selection for craniodental shapes with biomechanical attributes that permit ecological specialization for feeding on vertebrate tissues, or hypercarnivory.

Although incidences of morphological convergence are commonly cited as evidence of similar adaptation under similar selective conditions, such presumed linkages between convergence and adaptation require additional scrutiny and testing to rule out possibilities of selection of correlated traits unrelated to the adaptation or function of interest in producing the convergent characteristics observed, or coopting of traits that either were originally selected for other functions or have evolved from drift (8, 9). If selection for biomechanical specialization is indeed the principal driver of cranial shape evolution in carnivorans, we predict the following observations: Hypothesis 1, cranial shape should be more significantly associated with feeding ecology (diet breadth and trophic level) and locomotor and activity patterns (activity cycle and terrestriality/arboreality) that directly affect prey acquisition (hunting) and processing (mastication) and less closely correlated with life history or environmental variables; and Hypothesis 2, covariations between cranial shape and feeding ecology in turn should be significantly associated with variations in engineering-based measures of performance that distinguish between trophic levels (carnivory, herbivory, and omnivory) and between feeding specialists (those with narrow diet breadth) and generalists (those with broad diet breadth). These predicted observations constitute evidence that cranial biomechanical specialization provides functional and adaptational advantage for dietary specialization.

To assess whether selection for biomechanical specialization is indeed the principal driver of cranial shape evolution in carnivorans, we use engineering-based performance measures that have been correlated to ecology and biomechanical patterns in diverse mammal clades (1012), combined with theoretical shape morphing along ecological variable axes based on principles of geometric morphometrics (GMM) (1315), to investigate cranial shape changes correlated with ecological variables. We also conduct simulation-based tests of the biomechanical consequences of any shape-ecology correlations within a phylogenetic context. The performance variables assessed include cranial stiffness (measured by strain energy density), safety factor (measured by maximum von Mises stress), and mechanical efficiency (measured by ratio between output bite force at the canine teeth and input muscle force at the temporalis or masseter region), all previously shown to be important parameters in vertebrate musculoskeletal biomechanics (10, 16, 17). We developed and analyzed a high-resolution computed tomography (CT) digital cranial morphology data set across all living carnivoran families relative to nine well-sampled (known for >80% of species analyzed) ecological traits (table S3). Synthesizing multiple sources of species and environmental data, as a realization of the theoretical framework of Arnold (18), and applying current quantitative methods of shape and performance analyses (17), this paper describes functionally significant covariation between nonfeeding ecological variables and cranial shape in a major clade of mammals—the Carnivora (dogs, cats, seals, and relatives)—and proposes an additional mechanism that could explain evolutionary shifts in skull shape and biomechanical feeding specialization.


Covariances between skull shape and ecological variables were analyzed within a phylogenetic framework. In general, the same significance relationships between shape and ecological variables are observed in centroid size (CS) and feeding ecology using uniform branch lengths (main analysis) and in various sensitivity analyses incorporating fossil calibrations (for example, sensitivity analysis 6) or molecular clock node estimates (for example, sensitivity analysis 8) (see Supplementary Materials).

The principal phylogeny-associated results below are reported from the uniform branch length configuration of the main analysis, but differences in statistical analyses in other topological and branch length configurations are noted, with associated data presented in Supplementary tables (table S4) and figures (figs. S2 to S11). Comparing the main analysis to the sensitivity analyses, several additional significant associations between skull shape and ecological/environmental variables were gained (for example, in sensitivity analyses 2, 3, 6, and 8).

Size allometric trends of cranial shape variation are statistically significant for the all-Carnivora and feliform-only data sets, respectively (GMM data, G1. CS in Table 1). Cranial shape is significantly correlated with size (measured as CS; P = 0.03) in all Carnivora, but with the correlation driven mainly by feliform species (cats, hyenas, mongooses, etc.) with P = 0.04 and a nonsignificant correlation in analysis of caniforms only (Table 1). Among feliform species, size-driven trends in cranial shape became nonsignificant after removing the largest feliforms in the data set (Crocuta crocuta and Parahyaena brunnea); these two species appear to drive the all-Carnivora trend as well (feliform only, P = 0.23; all Carnivora, P = 0.07). Both species are large hyaenids that incorporate a substantial amount of hard, brittle foods (bone) in their diets (19, 20).

Table 1 Statistical significance (P values) of PGLS regressions based on resampling permutations (n = 1000).

Statistically significant (at the P = 0.05 level) regressions shown in boldfaced font. For combined regression analyses of multiple variables, none are statistically significant; P values of the individual factors are shown in those cases. For sensitivity analyses, see table S4.

View this table:

Species with larger CS have an expanded occipital crest and narrower cranial vault, as well as a broader rostrum (Fig. 1 and Table 2). This size allometry corresponds to a relative increase in the space between the sagittal crest and the zygomatic arches, the region occupied by the temporalis muscles. Size (CS), in turn, is significantly correlated with age at sexual maturity [Table 3; see the “Nonfeeding ecology (life history and environmental variables)” section below]. All subsequent phylogenetic regression and analysis of variance analyses were conducted between post–CS-regressed shape variables and the ecological variable of interest, followed by finite element (FE) simulations conducted on the data partitions with the most significant P values in each ecological variable category (or on the all-Carnivora data set if P values are equal).

Fig. 1 Cranial shape deviations associated with ecological traits.

Background arrows indicate direction of change depicted. Heat maps show range of negative deviations (cool colors) to positive deviations (hot colors). Deviation heat maps and shape deformations were performed on the consensus cranial shape.

Table 2 Shape changes by cranial region in each of the ecological gradients significantly associated with cranial shape covariation.

Category abbreviations as in Table 1.

View this table:
Table 3 Results of linear regression model fits of ecological variables that covary significantly with cranial shape, against other ecological variables.

df, degrees of freedom; Adj. SS, adjusted sum of squares; Adj. MS, adjusted mean square; F, F statistic; P, significance level.

View this table:

Feeding ecology, locomotor, and activity patterns

Of the nine ecological and life history traits analyzed, four are significantly correlated with cranial shape. Narrow diet breadth is correlated with a shortened rostrum and expanded frontoparietal region (especially an expanded postorbital constriction) (P = 0.01 for all Carnivora and P = 0.003 for feliforms in Table 1; Fig. 1 and Table 2); this comparison is only significant between species with narrow and broad diet breadths but not with species of intermediate diet breadth (also see fig. S13A). Trophic levels in the all-Carnivora data set and the feliform-only partition correlate with shifts in skull shapes (P = 0.03 for both). For herbivores, there are cranial shortening and widening from the consensus cranial shape for Carnivora; for omnivores, an elongation of the rostrum, narrowing of the postorbital constriction, and widening of the cranial vault are observed relative to the Carnivora consensus shape; and carnivores exhibit a shortening of the rostrum, expansion of the postorbital constriction, narrowing of the cranial vault, and widening of the zygomatic arches relative to the consensus cranial shape (P = 0.03; Fig. 1 and Table 2).

Cranial shape covariation between trophic levels is only significant for herbivores versus carnivores or omnivores; omnivores relative to carnivores are not significant. The omnivore-carnivore similarity is further evidenced in three-dimensional (3D) plots of the first three principal components (PC) axes, which show extensive overlap between the two groups in the morphospace (fig. S13B). None of the locomotor (terrestriality/arboreality and habitat breadth) or activity cycle variables are significantly associated with cranial shape variation in the principal phylogeny-associated analysis nor in any of the alternative topological and branch length configurations in the sensitivity analyses.

Nonfeeding ecology (life history and environmental variables)

Significant correlation between cranial shape and age at sexual maturity is present in Carnivora (P = 0.02), mainly driven by caniforms (dogs, bears, weasels, etc.; P = 0.03; Table 1 and fig. S12A). This correlation, statistically significant only in uniform branch length configurations, is associated with a dorsally raised rostrum, narrower postorbital constriction, widening of the cranial vault, and expansion of the occipital crests in species with relatively delayed age at sexual maturity (Fig. 1 and Table 2). Sexual maturity age itself is significantly correlated to maximum longevity, although longevity alone does not covary significantly with cranial shape (Tables 1 and 4). Last, mean monthly precipitation directly correlates significantly with cranial shape across all Carnivora in three of the additional topology/branch length configurations tested in the sensitivity analyses (P = 0.02 to 0.05; table S4 and fig. S12B), mainly driven by significant patterns within feliforms (P = 0.01 to 0.05; fig. S12B), with elongation of the rostrum, narrowing of the postorbital constriction, narrowing of the zygomatic arch, and narrowing of the cranial vault in species living in higher precipitation regions (Fig. 1 and Table 2). Precipitation is significantly correlated with another environmental variable, mean temperature, and with terrestriality/arboreality (locomotor mode), even though neither temperature nor terrestriality/arboreality independently showed significant covariation with cranial shape (Tables 1 and 3). Otherwise, in testing for autocorrelation among the shape-covarying ecological variables themselves, statistically significant pairwise correlations are lacking for all, except for between trophic level and mean precipitation, in which herbivores tend to live in regions of higher precipitation than carnivores, but neither of those trophic distributions is statistically different from that of omnivores.

Table 4 Engineering-based performance measure shifts associated with ecological traits.

Category abbreviations as in Table 1.

View this table:

Engineering-based performance measures

Theoretical cranial models representing the extremes of shape covariation, relative to the four significant ecological traits plus CS, were tested in three different sets of FE simulations. The first simulations, “nonmasticatory” dorsal bending tests, showed no significant differences in cranial stiffness (strain energy density, in Joules) or safety factor [maximum von Mises stress, in megapascals (MPa)] measures for any of the ecological traits (table S6). There is much higher uncertainty [wider confidence intervals (CIs)] in nonmasticatory simulations compared to masticatory simulations.

The second simulations, of temporalis muscle–driven bilateral canine bites, showed no difference in cranial stiffness in skull models representing shape changes associated with size allometry (CS increase). Increased overall stiffness is observed in four other ecological traits: with narrower diet breadth, earlier sexual maturity, lower precipitation, or shift toward herbivory, respectively (Table 4 and table S7). Cranial safety factor (maximum von Mises stress) also increased for each of those four ecological traits, whereas mechanical efficiency (ratio between output and input forces) increased with narrower diet and lower precipitation and decreased in the two other ecological shifts (earlier sexual maturity and herbivory).

In the third set of simulations, masseter-driven bilateral canine bites, no difference in stiffness is observed in CS-associated shape changes. Narrower diet and earlier sexual maturity are related to lowered stiffness, whereas lower precipitation and shifts toward herbivory are associated with higher stiffness. Safety factor exhibited no significant change with diet breadth or age at sexual maturity, but increased in the other two ecological traits plus CS. Finally, mechanical efficiency showed no significant difference with changes in diet breadth or age at sexual maturity; CS increase, lower precipitation, and relative shift from omnivory toward herbivory were associated with higher efficiency, and relative shift from carnivory toward herbivory resulted in lower efficiency (Table 4 and table S8).


Hypothesis 1

Our analyses demonstrate that cranial shapes in Carnivora not only are linked to feeding ecology patterns, as might be expected, but also correlate to life history and environmental variables, thus indicating that cranial adaptations for feeding ecology are closely associated with and may be constrained by interactions with life history and environment [cf., Eisenberg (21)]. Furthermore, cranial shape of ecological specialists (narrow diet breadth) can be distinguished from that of generalists (broad diet breadth) at the all-Carnivora level and in feliform species, and different trophic levels can be distinguished in the all-Carnivora and feliform-only partitions as well. These findings, part of the first comprehensive morphofunctional analysis of 3D cranial shape data, corroborate previous conclusions from 2D analyses with only a much more limited landmark coverage of carnivoran skulls (22, 23). Even with our higher data coverage of the cranial regions, using nearly 200 anatomical fixed landmarks and semilandmarks in 3D shape space, feeding ecology variables still constitute only part of an even larger suite of variables that influence the shape of carnivoran skulls; other ecological or phylogenetic factors also can explain some observed variation in carnivoran skull shapes (24).

Results demonstrate that differences in some life history traits, specifically age at sexual maturity, have a significant association with adult cranial shapes across the Carnivora (Table 1). Furthermore, mean monthly precipitation also is significantly associated with cranial shape, independent of feeding ecology (see below). However, lack of similar cranial shapes among taxa with the same feeding ecologies is not by itself evidence for lack of performance-based convergences in feeding adaptation. Instead, the presence of many-to-one, one-to-many, and many-to-many mapping of structure-function correlations across vertebrate structures suggests that convergence in dietary specialization and function (biomechanical performance properties) still could occur independently of any direct associations between skull shape and function (25, 26). As a case in point, although no significant cranial shape differences are observed between carnivore and omnivore trophic levels, significant differences do occur in both temporalis- and masseter-driven biomechanical performance (Fig. 2), suggesting that similar cranial shapes can nevertheless differ significantly in biomechanical traits. Such complex relationships may partly explain weak or absent correlations between skull shape and dietary factors in Carnivora (24). Therefore, our results support the recommendation that estimates of function or performance be included in tests of hypotheses that convergent skull shape suggests convergent feeding adaptation.

Fig. 2 Engineering-based performance measures analyzed in cranial models.

(A) Stiffness (strain energy density) in bending tests. (B) Safety factor [von Mises (VM) stress] in bending tests. (C) Stiffness in temporalis-driven tests. (D) VM stress in temporalis-driven tests. (E) Mechanical efficiency in temporalis-driven tests. (F) Stiffness in masseter-driven tests. (G) VM stress in masseter-driven tests. (H) Mechanical efficiency in masseter-driven tests. Error bars indicate 95% CI; if no error bars are visible, error is encompassed within the size of the mean value dots. Mpa, megapascals.

Hypothesis 2

Tests of cranial nonmasticatory and masticatory biomechanics using engineering-based estimates of performance indicate that all of the significant ecological covariates of cranial shape (CS, trophic level, age at sexual maturity, and mean precipitation) also show significant variation relative to biomechanical capability for mastication. The functional significance of both the temporalis muscle– and masseter muscle–driven bite simulations is evidenced by the fact that there are no significant differences when a nonmasticatory bending load is applied to the same cranial models (Fig. 2). This suggests that the overall mechanical strength of the carnivoran cranium has “built-in” levels of stiffness and safety factors that provide a basic level of mechanical protection regardless of feeding biomechanics, and this level of protection does not vary significantly even with substantial shape variation across the clade.

On the other hand, masticatory simulations show that performance gradients are present in cranial models associated with feeding-related ecological differences among species. Feliform ecological specialists with narrow diet breadth exhibit cranial shapes that are stiffer, more efficient, with a higher safety factor, and with emphasis on temporalis-driven bites. These biomechanical traits occur at the expense of a less stiff cranium for masseter-driven bites. In contrast, herbivory in Carnivora involves cranial shape changes that generate higher stiffness and safety factor for both temporalis- and masseter-driven bites, while deemphasizing mechanical efficiency. Carnivorous species show intermediate values of efficiency, stiffness, and safety factor between omnivory and herbivory in terms of temporalis-driven bites but have cranial shapes that show similar stiffness to that of omnivores and slightly higher efficiency in masseter-driven bites compared to both herbivores and omnivores (Fig. 2). This suggests that carnivory requires only moderately higher levels of biomechanical specialization relative to the requirements of herbivory. This observation makes sense when considering the material properties of the foods involved in each dietary category: Plants tend to require more mechanical force to masticate than do vertebrate soft tissues such as muscle (27).

Sensitivity of results to variations in phylogenetic tree topology and branch lengths

The simplest model we used, a phylogeny with uniform branch lengths, gives robust results relative to a varied array of sensitivity analyses incorporating branch lengths, particularly in significance of correlations of cranial shape with GMM/CS and feeding ecology parameters (see Supplementary Materials). Analyses incorporating branch lengths are sensitive to even small permutations in data inputs or assumptions (including fossil versus molecular clock age calibrations) or collapsing or expanding equal age nodes (using an equal-partition age assignment for expanded branches), and even resolving a phylogenetic polytomy in only one clade (family) can yield differences (losing significant correlations with feeding ecology for Carnivora and Feliformia, but gaining significance in correlations with precipitation for Carnivora and Feliformia). These analyses indicate that different ecological correlates of skull shape may be masked or highlighted depending on the particular phylogenetic model (tree topology and/or branch lengths) used in phylogenetic comparative analysis. Therefore, although the sensitivity analyses are largely congruent with the results from our main analysis and the major conclusions are similar among all analyses, the differences indicate that multiple analyses are useful in any study because model choice could have an influence on phylogenetic comparative analysis outcomes and inference if only a single phylogenetic model is applied.

Cranial performance and feeding ecology

The observed skull efficiency and stiffness trends in temporalis- and masseter-driven simulations reflect general feeding and behavioral patterns in carnivoran species. Feliform species categorized as dietary specialists (28) are either active hunters of vertebrates (Acinonyx jubatus, C. crocuta, Felis silvestris, Genetta piscivora, Neofelis nebulosa, and Panthera pardus) that could experience unpredictable loads during prey capture, or insectivores (Mungos mungo and Suricata suricatta) that have to pierce the hard and brittle exoskeletons of their prey. Both of those specialized ecomorphologies have biomechanical reasons to benefit from higher mechanical efficiency to deliver powerful killing or piercing bites and from stiffer crania with higher safety factors to withstand increased stresses from muscle and bite reaction forces compared to less specialized species.

The most specialized herbivores known among living carnivorans are the red panda (Ailurus fulgens) and giant panda (Ailuropoda melanoleuca); their convergent adaptations for consuming bamboo plants are well known, and functional convergence is supported by both genomic and biomechanical analyses (27, 29, 30). Although both hypercarnivores and herbivores have been shown to have higher bite forces (31), our results indicate that the underlying evolutionary transformations were achieved differently, with herbivores emphasizing masseter efficiency and carnivores emphasizing temporalis efficiency (Table 2). The presence of this biomechanical distinction in jaw muscle emphasis, together with the remarkable extent of morphological convergence between hypercarnivorous bone crackers and herbivores (27) and their association with stringent biomechanical requirements for mechanically demanding feeding tasks, suggests a prime example of complex mapping of structure to function in widely divergent feeding styles. Cranial morphological convergence in this case is associated with differential biomechanical emphasis depending on whether the feeding adaptation is to hard vertebrate tissues (including bone) or to tough plant fibers.

Cranial performance and nonfeeding ecological variables

Shape change associated with allometric size increase has the same covariation with cranial mechanical efficiency as a shift toward carnivory or herbivory, suggesting that selective pressure on overall size increase alone could drive trophic specialization. This is consistent with recent morphometric analyses that suggest size as a key factor in ecomorphological convergence in other mammal groups (26), and a similar pattern has been found in predatory birds (32). Similar allometry-to-diet linkages across distantly related vertebrate groups suggest that this may be a more general evolutionary phenomenon. Our new biomechanical findings also are consistent with the trend of body-size increase concomitant with morphological evolution toward hypercarnivory during Cenozoic carnivoran evolution (3), but the results here provide a new view grounded in 3D biomechanical data showing that functional benefits of size-associated shape changes, such as higher bite force, are accompanied by biomechanical changes that also enhance efficiency and cranial safety factor in similar ways as trophic level shifts. It is noteworthy that among the ecological variables, sexual maturity age is significantly correlated to size, suggesting that interactions between size and sexual maturity age covariation may amplify some biomechanical features (temporalis-driven safety factor and efficiency) and cancel out others during evolution (Table 3).

Temporalis muscle–driven performance gradients in cranial models of species from high versus low monthly precipitation habitats are exactly the same as the trend from broader to narrower diets in showing higher stiffness, safety factor, and efficiency. Similarly, the cranial model performance gradient from late to early sexual maturity age shares the same trend as specializations toward herbivory, with higher stiffness and safety factor but lower efficiency. Specific associations among variables also are observed in masseter muscle–driven simulations: Narrower diet and earlier sexual maturity share the same performance changes, whereas omnivory shifting toward herbivory shares the same trend as living in habitats with lower precipitation. Although unique performance profiles emerge when temporalis muscle– and masseter muscle–driven data are considered simultaneously, shared outcomes of functional changes in muscle-specific data partitions with nonfeeding ecological traits suggest that covariation that might have arisen from selection for those nonfeeding ecological traits could result in cranial shapes that are highly similar biomechanically to the observed outcomes due to feeding specializations. From a methodological perspective, the differences in simulation outcomes between bending, temporalis-driven, and masseter-driven loading scenarios show that output values are highly sensitive to loading condition. These results highlight the importance of estimating physiologically relevant and question-specific loading conditions when using FE simulations to test hypotheses of structure-performance relationships.

Previous research shows that caniform species exhibit more heterochronic shifts in cranial suture fusion than feliform species (33); this could provide one explanation for our finding that significant performance differences are correlated with cranial shape differences in adult specimens with varying age at sexual maturity (Table 4). Cranial suture closure in mammals marks cessation of bone growth around the suture, and heterogeneous distributions of suture closure timing across craniofacial sutures might affect overall cranial shape in adults, which then could directly influence biomechanical variables related to mastication. Sexual maturity age is correlated with maximum longevity, consistent with prolonged heterochrony as a source of cranial shape variation (Table 3). The cranial regions with the largest shape deviations in the gradient of younger to older ages at sexual maturity across species occur at the dorsal rostrum and occipital areas, which contain several cranial sutures at the extremes of fusion timing documented in carnivorans (33).

Genes such as Runx2 (runt-related transcription factor 2) that have previously been shown to be correlated to carnivoran facial length via heterochrony of ossification (34) represent potential genetic controls of overall cranial shape change via developmental heterochrony. As a test of this genetic control interpretation, we conducted phylogenetic generalized least squares (PGLS) regression analyses of Runx2 glutamine-alanine tandem repeat ratio against cranial shape for 16 taxa that overlapped with our data set. The analyses returned no statistical significant covariation (P = 0.49 to 0.80), but the broad association between shape and ontogeny that we document in this analysis warrants additional study to determine the relationships among cranial suture closure timing, adult cranial shape, genetic control of bone ossification and growth, and biomechanical variables relevant to ecological specialization. It is noteworthy that the statistically significant correlations for age at sexual maturity appearing in only a few of the topological and branch length permutations of the phylogenetically constrained analyses (Table 1 and table S4) indicate that age at sexual maturity is more weakly associated with cranial shape than is either cranial CS or feeding ecology in Carnivora. As such, it may be more difficult to identify specific developmental (heterochronic) mechanisms for explaining cranial structure-function transformations than for feeding ecology variables.

An association between precipitation and feeding ecomorphology has been demonstrated in diverse vertebrate groups [(3537) and references therein]. Field studies of primates also show significant influences of precipitation levels on food quality, which in turn can significantly affect fitness (measured by infant mortality rate) via variations in masticatory efficiency between individuals at different tooth wear stages (38, 39). Although these associations alone do not constitute direct causal factors for how precipitation influences cranial shape, we did find significant correlations between precipitation and both temperature and terrestriality/arboreality (Table 3). Species classified as “aboveground” dwelling tend to occur in high-precipitation, high-temperature biomes such as tropical forests and savanna woodlands (fig. S14). The correlation between aboveground- or arboreal-dwelling species and forested environments with high precipitation appears to hold at the faunal level (40), suggesting that this general relationship could explain our findings. Most research on morphological proxies of arboreality in mammals have focused on postcranial anatomy [for example, the studies by Hunt (41) and Dublin (42)], and studies of orbital orientation found no significant ordinal-level trends of cranial shape covariation with arboreality or other ecological variables (43, 44). Because we found that significant correlations between precipitation and terrestriality/arboreality encompass species with different diet breadths and trophic levels, the causal factors behind cranial shape covariation with precipitation-related aboveground dwelling or arboreality are likely common biological factors that vary across, rather than within, those feeding categories.

In the absence of comparative feeding performance data (for example, bite force) to directly test the structure-function relationship associated with precipitation levels, we identified the most biomechanically demanding food items or specific prey hunting behavior known in the species studied, using natural history accounts (45, 46) and relative biological material properties (47), to then test for correlations between the identified structure-function relationship and mechanical demand, as a proxy for feeding performance. In a post hoc test, we assessed the potential association between precipitation and mechanically demanding food items or hunting behavior (table S9) using PGLS regression and found no significant correlation between them in either the uniform branch length main analysis or any of the alternative sensitivity analyses (P = 0.08 to 0.68). In fact, none of the other continuous ecological (sexual maturity and longevity) or environmental (mean temperature) variables (sexual maturity age, P = 0.08 to 0.40; maximum longevity, P = 0.12 to 0.81; and temperature, P = 0.31 to 0.91) are significantly correlated with the presence of mechanically demanding food items or with specific hunting behavior, and cranial shape is very weakly correlated with mechanical demand (P = 0.07 to 0.93 across 10 topologies and branch lengths, except in sensitivity analysis 8, with a significance of P = 0.03) across the data set as a whole. The biomechanically significant cranial shape covariation with precipitation levels therefore may be related to nonfeeding biomechanical demands along the arboreality-precipitation spectrum, such as reducing overall skeletal weight of the head in arboreal settings relative to terrestrial species, or possibly from nonmechanical factors altogether, such as skull changes related to accommodation of modified sensory systems (visual, balance, olfaction, etc.) in arboreal and forested environments versus open habitats. Another potential link worthy of further investigation is the association between biodiversity and precipitation (or precipitation and temperature combinatorially). Two-thirds of the world’s biodiversity hotspots are located in high-precipitation forests (48), most of which are also warm, and constraints on cranial shape across a community of species (feeding-related or not) may be important in understanding how precipitation-related trait variation in ecological and life history variables causally affects cranial shape via species interactions in high-diversity environments.

Considered in sum, these results demonstrate the possibility of multiple evolutionary pathways that can result in similar changes in mastication-related performance traits, independent of feeding ecology: (i) Shape changes associated with allometric size increase mimic shifts toward adaptation for herbivory and carnivory; (ii) shifts from higher to lower annual precipitation environments show temporalis-driven biomechanical changes that mimic narrowing diet breadth and masseter-driven improvements that mimic shifts between herbivory and omnivory; and (iii) a shift toward earlier age at sexual maturity mimics biomechanical changes associated with herbivores relative to omnivores and carnivores. These associations suggest that direct selective pressure for those masticatory performance variables over evolutionary time across carnivoran lineages is neither exclusive nor necessary in achieving cranial morphological shifts toward biomechanical specialization (or generalization, depending on the direction of selection). At its hypothetical extreme, masticatory performance–changing cranial shape evolution could simply be a by-product of an overall covariation of species cranial shape with environmental or life history changes (Fig. 3). This covariation in turn may be related to a diverse array of nonfeeding selective pressures such as sexual selection for heterochrony in age at maturity, selection for sensory adaptations in arboreal environments, or other factors currently lacking sufficient data to test. The resulting biomechanical specialization then could lead to ecological specialization by providing the necessary morphological variation to be selected upon later, in association with ecological influences (4). Sexual selection, specifically male-male combat, has been identified as a key factor explaining the biomechanical properties of pinniped mandibles (49); the results of our analyses suggest that similar influences from nonfeeding ecological and environmental variables on carnivoran cranial structure-function are likely to have been much more widespread than previously recognized. The fact that we found little statistical support for further correlations between shape-correlated feeding and nonfeeding ecology variables suggests that there can be dietary convergence and possibly biomechanical performance convergences without the typically assumed, diet-driven skull shape convergences.

Fig. 3 Summary figure with directions for future research.

A flowchart represents the relationship between the environmental, life history, ecological, shape, and feeding performance variables analyzed in this study. Asterisks denote areas that might especially benefit from additional research to understand how life history and environmental variation directly influence cranial shape variation. Shaded arrows indicate likely feedback pathways that were not specifically analyzed in this study but, nevertheless, are important for future research into the relationships among feeding ecology, biomechanical performance, and cranial shape.

Evolutionary potential of cranial shape covariation

The significant association of sexual maturity and precipitation-related arboreality with masticatory performance suggests an additional explanatory mechanism for the macroevolutionary trend of iterative hypercarnivory in carnivoran evolution. Precipitation has been shown to significantly influence Cenozoic fossil mammal community composition through water-mediated environmental differences (50), but few mechanisms have been proposed or demonstrated that could directly link precipitation with changes in masticatory biomechanics [but see the studies by King et al. (38, 39)]. Our results show, contrary to previous hypotheses of climate-driven evolutionary specialization of Cenozoic carnivorans (5, 7, 23), that direct selection for size and feeding biomechanical specialization during the overall cooling and drying climatic trend of the Cenozoic is not an exclusive mechanism for the observed pattern of biomechanically modified cranial shapes in hypercarnivorous lineages. A decrease in precipitation (and corresponding decrease in aboveground dwelling or arboreal) correlates with the same changes in cranial biomechanics as would a shift from broader to narrower diets or from omnivory toward herbivory or carnivory (Fig. 2 and Table 2). Lower precipitation correlates with cranial shape changes that provide across-the-board increases in stiffness, safety factor, and mechanical efficiency of both temporalis- and masseter-driven mastication. Therefore, compared to the functional changes observed in across-feeding ecological variables (diet breadth), cranial shapes of species at the lower end of the precipitation gradient experience more biomechanically advantageous changes relative to high-precipitation species than observed between broad and narrow diet breadth taxa.

These results lead us to propose that both developmental (size and age at sexual maturity) and environmental (precipitation-related ecological variables) shifts played a more important role in generating cranial morphofunctional variation critical in the evolution of carnivoran feeding specialization than previously envisioned for carnivoran evolutionary history. Sexual selection and environmental changes could lead to evolutionary shifts in feeding-relevant cranial morphology that interacts with, but does not directly require, evolving predator-prey relationships or shifts in other feeding ecology parameters.

Implications for evolution of dietary specialization

The evolution of hypercarnivory in Carnivora is associated with decreased morphological disparity (51). Decreases in morphological disparity have been related to the increases in morphological trait integration, given the premise that selective pressure on specialized structural systems channels an increase in integration and a reduction of trait variability to maintain functionality during morphostructural changes in evolution (52). Evidence for the relationship between trait disparity and integration subsequently was considered weak in Carnivora, in part because of other potential sources (for example, environment) affecting disparity (53). However, our findings support the assertion that factors affecting trait disparity, such as in feeding ecology, need not come from a direct association between the function in question (for example, mastication) and the morphological structures of interest (for example, cranial shape). Instead, covariation between morphological shape changes (and its potential effects on morphological disparity) and environmental variables during periods of strong and/or persistent directional environmental change may provide an indirect pressure driving evolution of masticatory performance via trait covariation, thereby indirectly influencing trait disparity.

Changing prey bases and hunting strategy responses may still represent important sources of direct selective pressure for ecological and biomechanical specialization, but, as our results demonstrate, cranial shape covariation with nonfeeding ecological variables alone can significantly affect the same cranial biomechanical performance measures that distinguish diet breadth and trophic levels in Carnivora as a whole. The significant correlations identified between some ecological and life history variables and skull shape suggest that those shape changes correlated with nonfeeding ecological variables could be involved in the evolution of feeding specialization in several ways (Fig. 3). Shape correlation arising from nonfeeding ecology traits could provide the variation upon which selection can act to alter cranial morphology feeding performance relationships, or as sources of constraint on those variations, reducing the number of possible morphofunctional outcomes. The direct causal mechanisms of cranial morphology variation resulting from life history and environmental variation remain to be identified and constitute an important future research direction to clarify the evolutionary mechanisms underlying observed structure-performance trends in Carnivora (Fig. 3).


We discovered significant influences of nonfeeding ecological variables in masticatory biomechanical performance, and our findings demonstrate that direct associations of putative feeding ecology drivers for cranial morphology alone are not sufficient to infer coupling of cranial structure and function—other, nonfeeding drivers can change morphology in ways that result in biomechanical modifications relevant to feeding adaptations. Continued research on both the causal functional factors underlying cranial shape covariation with nonfeeding ecological variables, and the potential feedback on cranial shape–ecology covariation imposed by biomechanical and other functional trade-offs (54), would further clarify the relative importance of life history and nonfeeding and feeding ecology factors in the evolution of carnivoran feeding specializations (Fig. 3). From a biomechanical and functional perspective, our new findings suggest that mechanisms mediating the iterative evolution of carnivoran feeding ecomorphologies during Cenozoic environmental changes may be more complex than permitted by existing, widely accepted, size and feeding ecology–dominated narratives. Furthermore, because the performance variables we analyzed represent attributes shown to be important in the evolution of musculoskeletal biomechanics in general, similar processes would be expected across a wide spectrum of vertebrate groups.


Our data set consisted of 53 species of Carnivora; we sampled across all living families with the goal of representing the broad taxonomic and morphological diversity in the order (table S1). Adult, wild-caught skulls were used preferentially for analysis whenever available. Skulls were high-resolution CT–scanned using the GE v|tome|x s scanner at the American Museum of Natural History (AMNH) for model construction; where available, published models from previous studies were used instead of rescanning those specimens (15, 30, 5558). 3D models of crania were constructed from the CT scans using Mimics Research (Materialise) and exported as *.stl files. One hundred ninety fixed landmarks and semilandmarks were placed on each cranium using Landmark Editor and exported as text files for analysis using the R programming environment (packages Geomorph and Morpho) (59). The landmarks covered specific homologous anatomical regions (as fixed points) and analogous regions (as a grid patch, each containing 15 landmarks covering the muscle attachment sites for the temporalis and masseter muscles, frontal region, maxillary regions of the rostrum, and hard palate) (table S2) (60). Because semilandmarks could be either treated in the same way as fixed anatomical landmarks or adjusted (or slid) along designated curves and surfaces, separate analyses were conducted treating either all as fixed or each as containing six sliding points and nine anchoring points on each grid patch (fig. S1 and table S2). Fixed semilandmarks were used “as is” with the fixed anatomical landmarks in subsequent superimposition, whereas the sliding semilandmark data sets were further analyzed in two ways: using the bending energy or Procrustes criteria. Results were compared among the three methods of treating semilandmarks. No significant differences were observed in the resulting morphospace created by the superimposition using generalized Procrustes analysis followed by PC analysis. Therefore, the fixed landmark data set, which makes the simplest set of assumption that all landmarks are fixed, was used for all subsequent analyses.

Phylogenetic context

A phylogenetic tree of all 53 species examined was generated on the basis of the most recent comprehensive molecular phylogeny of Carnivora (61), with all species interrelationships resolved (other than a polytomy of Eupleridae species) and used as the “base phylogeny” for the main analysis subsequent results and discussion (fig. S2). Reliable branch length and divergence age data are not available for all of the species or nodes/clades examined, so this base phylogeny (with a conservative assumption of all branch lengths being equal/uniform) was used in subsequent comparative analyses. However, recognizing that branch lengths can influence ancestral state reconstructions and inferences of rates of change and magnitudes of evolutionary changes (62), we implemented a sensitivity analysis to assess the potential impacts of alternative topologies and incorporation of varying branch length estimates on the statistical outcomes obtained in subsequent analyses, using nine different alternative topology/branch length configurations (figs. S3 to S11 and Supplementary Materials text). These sensitivity analyses assessed the effects of assigning branch length estimates using a variety of methods, data sets, and assumptions, including resolving a polytomy for Eupleridae species, incorporating alternative suites of fossil calibrations versus molecular clock divergence age estimates, and collapsing branches with equal divergence age estimates or expanding the resultant age polytomies by partitioning the age difference equally along each branch between the earliest and latest diverging nodes. Phylogenetic topologies and branch lengths were constructed in Mesquite and exported to R in the simple Newick format. Nine ecological traits of carnivoran species were extracted from the PanTHERIA database of ecological traits for mammals (28). The traits analyzed were chosen on the basis of completeness of coverage across the 53 species examined in this study; only traits that are recorded for more than 80% of the species (or n ≥ 42 species) were included in the analysis. The ecological traits are activity cycle, terrestriality/arboreality, habitat breadth, diet breadth, trophic level, longevity, age at sexual maturity, mean precipitation, and mean monthly temperature.

A series of PGLS analyses were conducted to determine the relationships between cranial shape and ecological traits, while simultaneously taking into account the influence of phylogenetic nonindependence on the covariation among residuals. The superimposed shape data set was first regressed against cranial CS to check for size allometry. The residuals from the shape-CS PGLS analysis were extracted and used as size-free shape data in all subsequent PGLS analyses between shape and ecological variables. PGLS regression analyses were conducted for all nine ecological traits for the Carnivora and on separate caniform and feliform partitions of the carnivoran taxa data set, respectively. Additional PGLS analyses were conducted between ecological variables found to be significantly (P ≤ 0.05) associated with cranial shape.

Theoretical shape morphing

The cranial shape variations represented by significant ecological correlations (group means for categorical data; extreme values for continuous data) were exported as landmark deviations from the consensus cranial shape of the entire Carnivora data set for functional simulation purposes. For continuous variables (cranial CS, age at sexual maturity, and mean precipitation), the two species with the maximum and minimum values were identified, respectively, and the landmark deviations corresponding to the PGLS-fitted values of those two species were used to “morph” the consensus shape model into theoretical cranial models. For discrete variables (dietary breadth and trophic level), the landmark deviations of the categorical consensus shape from the global (whole data set) consensus shape model were used for morphing the consensus shape model into theoretical models. A PC analysis of the superimposed shape data set was conducted, identifying the Madagascar fossa Cryptoprocta ferox as the species closest to the consensus shape for all Carnivora. Therefore, all theoretical shape morphing operations and analyses were conducted using Cryptoprocta cranial morphology as the average shape. Shape deviations associated with the ecological traits were mapped onto the consensus cranium mesh (3D model constructed from primary scan data) represented by Cryptoprocta by using the 190 landmarks to deform the consensus shape model according to a thin-plate spline criterion implemented in Landmark Editor (13).

The anatomical landmarks used in this study were distributed only over the surface of the skull (that is, there were no internal landmarks); thus, the morphed models only captured variation in surface shape, not variations in internal thickness of bones. This landmark arrangement produced relative overthickening of small-CS shape models (and correspondingly overestimating performance) and overthinning of large-CS shape models (underestimating performance). This approach has the effect of reducing the performance differences that would otherwise be observed between them, thereby making comparisons of models on the CS spectrum more conservative/homogenized in terms of CS-related bone thickness variation. All deformed meshes then were saved and imported into Geomagic Wrap (3D Systems), where they were remeshed/improved to remove triangular plate elements that were highly skewed by the deformation process. Metric deviations between meshes representing categorical means, or extremes of continuous trait variables, were visualized using the “deviations” function in Wrap.

FE simulations

The improved cranium meshes then were saved as stereolithography (*.stl) files and imported into Strand7 (Strand7 Pty Ltd.) FE analysis software for solid meshing. Using the protocol of Tseng and Flynn (63) to account for systematic variation in FE-measured performance outcomes, low-, medium-, and high-resolution, four-noded tetrahedral meshes were generated for each of 11 theoretical cranium meshes (representing theoretical morphotypes spanning the extremes for ecological variables assessed: large size, small size, broad diet, narrow diet, early maturity, late maturity, high precipitation, low precipitation, trophic carnivore, trophic herbivore, and trophic omnivore) (table S5). The performance variables measured (mechanical efficiency, safety factor, and strain energy) from each theoretical mesh were averaged over the different resolution models, and 95% CIs were calculated using the t distribution (because of small sample size) to assess statistically significant differences between species models by determining a lack of overlap between the 95% CI ranges in a given pairwise comparison. Three sets of FE simulations were conducted: First, a nonmasticatory dorsal bending test was implemented by fixing one node each at the left (zero degrees of translational freedom) and right (a single degree of translational freedom, in the long axis connecting the two constraint nodes) occipital condyles and by placing dorsally directed nodal forces at both canines, with each nodal force (in Newtons) equivalent to 1% of the numerical value of the total surface area (in square millimeters) of the models (for example, a model with a surface area of 100 mm2 would have a nodal force equal to 1 N) as a biologically reasonable estimate. Nodal force magnitude was arbitrarily chosen because all comparisons are relative, not absolute; the set ratio of nodal force as 1% of total surface area of the models satisfies the isometric requirement for comparing stress values between FE models with different absolute areas (64), which requires the input force–to–area ratio to be constant across models being compared. This bending simulation provides a nonphysiological, engineering-based strength test to assess the general performance of the structure.

The second set of simulations represented a temporalis muscle–driven bite. Regions of the temporalis muscle attachment were identified on all models based on comparative anatomy, and arbitrarily assigned, biologically reasonable forces equivalent to 1% of total model surface area were then applied equally to both left and right sides using Boneload (65). Muscle forces were oriented toward centroid sites on the insertion sites of the mandible; centroids of the mandibular muscle insertion sites were calculated using the “center of gravity” function in Wrap. Cryptoprocta was used as the reference mandible for all theoretical models; that reference mandible was scaled linearly to fit each theoretical model by adjusting its length to create occlusion between the upper and lower incisors and to close articulation of the temporomandibular joints (TMJs), and by adjusting width so that the first lower molar and last upper premolar (the carnassial pair) were articulated. The mandible then was propped open at a 30° angle to represent a typical gape for Carnivora. Single nodal constraints were placed at the TMJs (that is, one side with no translational freedom, the other side with a single degree of translational freedom in the long axis of the TMJ), as in the study by Dumont et al. (66), and at both tips of the canine to simulate a bilateral canine bite. The canines were the only homologous teeth present across all 53 species in this analysis and, thus, were the only teeth that were associated with homologous fixed landmarks across taxa. Because no homologous landmarks could be placed at the cheek teeth, we interpreted engineering-based performance measures using only bilateral canine simulations and did not test premolar or molar bites.

The third set of simulations was essentially identical to the second set, except that the masseter muscles were activated instead of the temporalis. These two major jaw-closing muscles (masseter and temporalis) were analyzed separately because data were not available on relative physiologic cross-sectional area across the entire species data set. All FE models were theoretical deformations that do not correspond to actual species. By analyzing temporalis and masseter bites separately, we were better able to identify cranial shape changes that differentially influenced the biomechanical effectiveness of each of those muscle groups. All models were assumed to be isotropic as a simplification (11, 67), with a Young’s modulus of 20 GPa and a Poisson’s ratio of 0.3; homogeneous models provide simplifications of bone properties for a complex and dynamic tissue (67), but simplified cranial models, nevertheless, have been demonstrated to replicate general stress and strain distribution patterns in ex vivo validation experiments of mammal skulls, as well as closely estimating empirically measured bite forces (11, 12, 68). Three engineering-based performance metrics previously shown to be correlated to ecological traits (10) were used: (i) total strain energy density (as a measure of stiffness; lower strain energy means higher stiffness of structure), adjusted using the equation of Dumont et al. (64) to account for model differences in surface area and volume; (ii) maximum von Mises stress (as a measure of structural safety factor; the lower the maximum stress, the higher the safety factor), with the top 5% of values excluded to account for the effect of point loads creating artificially high stress values in close proximity to nodal constraints (68); and (iii) mechanical efficiency (defined as the ratio of output bite force measured at the canine tips to the total input force at the muscle attachment regions). Statistical significance was assessed by the lack of overlap between the 95% CI ranges of the models (no overlap = statistically significant difference at the P ≤ 0.05 level). Statistical analyses were conducted in R, Microsoft Excel, or Minitab.


Supplementary material for this article is available at

Additional Methodological Details and Results

table S1. List of specimens used in the study.

table S2. List of anatomical landmarks and semilandmarks used in the GMM analyses.

table S3. Ecological attributes of species analyzed.

table S4. Summary of sensitivity analyses of PGLS regression using different phylogenetic tree topology and branch length configurations.

table S5. Area and volume attributes of the theoretical morphed models analyzed.

table S6. Output values of FE simulations in the “bending test” analysis.

table S7. Output values of FE simulations in the “temporalis-driven” analysis.

table S8. Output values of FE simulations in the “masseter-driven” analysis.

table S9. List of the highest mechanical demands from food items experienced by each carnivoran species, ranked on the basis of stiffness and density.

fig. S1. Anatomical landmarks used in the GMM analyses.

fig. S2. Main analysis.

fig. S3. Sensitivity analysis 1.

fig. S4. Sensitivity analysis 2.

fig. S5. Sensitivity analysis 3.

fig. S6. Sensitivity analysis 4.

fig. S7. Sensitivity analysis 5.

fig. S8. Sensitivity analysis 6.

fig. S9. Sensitivity analysis 7.

fig. S10. Sensitivity analysis 8.

fig. S11. Sensitivity analysis 9.

fig. S12. Contour plot of the first two PC axes against life history and environmental variables.

fig. S13. 3D plots of PC scores of cranial shape data.

fig. S14. Plot of mean annual precipitation versus mean annual temperature of species studied.

References (6972)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: M. Hill, H. Towbin, and C. Grohé assisted with CT scanning at the AMNH and NSF-funded AMNH Microscopy and Imaging Facility. E. Westwig and N. Duncan assisted with mammalogy loans. D. Su, B. Van Valkenburgh, T. Rowe, and the Digimorph library provided additional CT scans. J. Finarelli provided some fossil calibration dates for the sensitivity analysis. S. Naeem, an anonymous reviewer, and D. Polly provided highly encouraging and constructive comments that improved the articulation of the overall framework and implications of the study. Funding: This research was funded by NSF DEB (Division of Environmental Biology) 1257572 and the Frick Fund (Division of Paleontology, AMNH). Author contributions: Z.J.T. and J.J.F. conceived the study, collected the data, conducted the analyses, and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. All FE models (in Strand7 format) generated from this study and all alternative phylogenetic configurations tested in the sensitivity analyses are stored at the Dryad repository (doi:10.5061/dryad.4t5q5). All numerical data used in making the figures and tables are included in the Supplementary Materials file. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article