Complex adaptive landscape for a "Simple" structure: The role of trade-offs in the evolutionary dynamics of mandibular shape in ground squirrels.

Trade-offs are inherent features of many biomechanical systems and are often seen as evolutionary constraints. Structural decoupling may provide a way to escape those limits in some systems but not for structures that transmit large forces, such as mammalian mandibles. For such structures to evolve in multiple directions on a complex adaptive landscape, different regions must change shape while maintaining structural integrity. We evaluated the complexity of the adaptive landscape for mandibular shape in Marmotini, a lineage of ground squirrels that varies in the proportions of seeds and foliage in their diets, by comparing the fit of models based on traits that predict changes in mandibular loading. The adaptive landscape was more complex than predicted by a two-peak model with a single dietary shift. The large number of adaptive peaks reflects a high diversity of directions of shape evolution. The number of adaptive peaks also reflects a multiplicity of functional trade-offs posed by the conflicting demands of processing foods with various combinations of material properties. The ability to balance trade-offs for diets with different proportions of the same foods may account for diversification and disparity of lineages in heterogeneous environments. Rather than constraints, trade-offs may be the impetus of evolutionary change.

Trade‐offs are present in any system where the same morphological change can have positive and negative functional consequences. Trade‐offs have been the focus of many studies in evolutionary morphology, from analyses of simple lever arms (Herring and Herring 1974; Bramble 1978; Greaves 1978; Emerson and Radinsky 1980) to multielement linkage systems (Liem 1973; Westneat 1994; Alfaro et al. 2004) and other complex structures (Vermeij 1973; Santana et al. 2011; Stayton 2011). The simplest trade‐offs are fundamental properties of the mechanical system, such as the inverse relationship between mechanical advantage and velocity ratio of a lever. One such system is the mammalian mandible, which is often treated as a class 3 lever with the muscle resultant between the bite point and the joint (Herring and Herring 1974; Hylander 1975; Dumont 1997; Taylor 2002; Young and Badyaev 2010). All other things being equal, the optimal mandibular length might be determined by the balance of selection pressures on bite force (mechanical advantage) and closing speed (velocity ratio). In more complex systems, there may be multiple trade‐offs associated with the interacting components, but similar principles can be applied to infer optimal morphologies as weighted averages of theoretical ideal shapes (Shoval et al. 2012).

Trade‐offs have been seen as constraints on evolution because any potential change to a system with trade‐offs will have both positive and negative functional consequences. Only extreme specialists restricted to performing a function that requires only one of the conflicting attributes will be able to evolve unimpeded toward the optimal morphology for that function. Functional decoupling, allocating functions with conflicting demands to different parts, has long been seen as a path to circumvent trade‐offs (Liem 1973; Lauder 1981; Alfaro et al. 2004). Each part is still subject to its own trade‐offs, but the combination may be able to undergo complimentary changes that minimize the most severe negative consequence, exchanging that cost for another that is less onerous. Decoupling the variation of functionally related parts also has the benefit of expanding the potential morphospace of the lineage because it increases the number of possible combinations of variations of the independent parts (Vermeij 1973).

Decoupling is not without its own trade‐offs. One of those trade‐offs is between the ability to generate independent variation permitting novel combinations and the increased potential to generate variations that are not adaptive and that may impede evolution in the direction of selection (Bürger 1986; Wagner 1988). This cost may be reduced in systems in which patterns of variation are themselves variable or flexible (Hansen and Houle 2008; Marroig et al. 2009). Another trade‐off of decoupling is the loss of structural integrity due to the elimination of physical connections. When one of the primary functions of a structure is to sustain or transmit large forces (crushing jaws, armor), an increase in mechanical independence may not be functionally viable. The advantages of structural continuity for resisting mechanical deformation may explain the progressive simplification of the lower jaw in the amniote lineage leading to mammals (Kemp 1972; Rubidge and Sidor 2001).

The simplicity of the mandibular structure in mammals has not prevented the diversification of its shape in multiple lineages (Monteiro and Nogueira 2009; De Esteban‐Trivigno 2011; Prevosti et al. 2012; Morales‐García et al. 2021; Wang et al. 2021). Several factors may have contributed to this evolutionary flexibility, including the development of the mandible from multiple primordia (Bhaskar 1953; Atchley and Hall 1991; Tomo et al. 1997; Ramaesh and Bard 2003), a plethora of genetic factors influencing development at multiple spatial scales (Gaunt 1964; Herring and Lakars 1981; Zhao et al. 1994; Depew et al. 1999; Ruest et al. 2003; Anthwal et al. 2008), and the multiple localized zones of deposition and resorption that mold the adult shape (Robinson and Sarnat 1955; de Buffrénil and Pascal 1984; He and Kiliaridis 2003; Mavropoulos et al. 2005; Sun and Tee 2011). Thus, there is a suite of localized signaling factors that have the greatest effect within narrow anatomical domains and that could provide the decoupling necessary for mandible shapes to diverge and diversify as positions and relative sizes of muscles and teeth change.

Regionalized gene expression, along with variation in positions of muscles and teeth, may have provided the potential to evolve novel shapes in response to changes in the balance of conflicting functional demands, enabling the evolutionary diversification of mammalian mandibles. The main lineages of tree squirrels may not have exploited that potential, instead taking advantage of the functional versatility of their morphology (Zelditch et al. 2020), but ground squirrels in the tribe Marmotini have undergone more substantial shifts in diet and habitat (Thorington et al. 2012; Zelditch et al. 2015), which may have provided the impetus for divergence to new adaptive peaks. Marmotini encompasses the majority of extant ground squirrel species, including all ground squirrels in temperate or colder grasslands of North America, northern Asia, and Europe (Thorington et al. 2012; Zelditch et al. 2015). Phylogenetic relationships indicate that the ancestral diet of Marmotini was primarily granivorous and that two branches shifted to diets that include large quantities of leaves from grasses and forbs (Zelditch et al. 2017). An analysis of this group demonstrated that granivores and folivores occupy different adaptive peaks (McLean et al. 2018). In this study, we tested additional hypotheses to determine whether those groups might occupy multiple adaptive peaks, reflecting different resolutions of the trade‐offs posed by the changes in diet and habitat. We used a likelihood‐based approach to infer the number of optima for mandibular shape (Clavel et al. 2015), with parametric bootstrapping to assess differences in model fit (Boettiger et al. 2012). The models in this study included the simple models commonly used in studies of evolutionary dynamics (single rate Brownian motion – BM‐1, early burst – EB, and single peak Ornstein‐Uhlenbeck process – OU‐1) and more complex models based on four traits expected to predict adaptive divergence in mandible shape (Diet, Habitat, Size and Tooth morphology). From the best fitting models and the phylogeny, we inferred the directions of shape changes that occurred in the transitions between optima. Then, we examined these patterns of shape change to identify differences in the combinations of shape changes in anatomically and functionally distinctive regions of the mandible.

Methods

These analyses use a subset of the data published by Zelditch et al. 2015; https://doi.org/10.5061/dryad.kq1g6) consisting of coordinates of 14 landmarks and 84 semilandmarks collected from mandibles photographed in lateral view (Fig. 1). We also used the phylogeny reported in that study, which was based on five mitochondrial genes (16S, 12S, COII, COIII, and Cyt‐b) and three nuclear genes (C‐myc, IRBP, and RAG1), and covered 66% of extant sciurid species and several outgroups. In brief, the tree topology and divergence times were estimated simultaneously in BEAST using substitution models obtained from PartitionFinder and calibration points as in Mercer and Roth (2003); following Markov chain Monte Carlo analyses, a maximum credibility consensus tree was generated from the sampled trees. Further details may be found in Zelditch et al. 2015). The tree was pruned to include only the Marmotini for which morphometric data were available (Fig. 2). The data set for those species includes 750 adults representing all 13 recognized extant genera (all monophyletic) and 75 of the 95 included species (Table 1). Most genera are represented by at least 67% of their species; only Spermophilus is represented by less than half of its species (7 of 15).

Figure 1

Landmarks (large blue circles) and semilandmarks (small green circles) on a representative specimen of Callospermophilus lateralis. Bar = 1 cm.

Figure 2

Phylogenetic relationships of included taxa, following Zelditch et al. 2015.

Table 1

Taxonomic sampling. Numbers of species included in the study compared to the Total number of recognized extant species (Thorington et al. 2012) and as a proportion of that total for the whole clade and for all currently recognized genera

Taxon	Number	Total	Proportion
Marmotini	75	95	0.80
Sciurotamias	1	2	0.50
Tamias	25	25	1.00
Notocitellus	2	2	1.00
Ammospermophilus	3	5	0.60
Callospermophilus	3	3	1.00
Otospermophilus	3	3	1.00
Marmota	10	15	0.67
Spermophilus	7	15	0.47
Urocitellus	10	12	0.83
Poliocitellus	1	1	1.00
Ictidomys	2	3	0.67
Xerospermophilus	3	4	0.75
Cynomys	5	5	1.00

All analyses described in this study were performed using packages in R (R_Core_Team 2018–2021). Coordinates of landmarks and semilandmarks were superimposed by Generalized Procrustes analysis, sliding semilandmarks to minimize bending energy (Green 1996; Bookstein 1997; Zelditch et al. 2012), using geomorph (Adams et al. 2021). Size was measured as the centroid size of the jaw, which is highly correlated with body size (Zelditch et al. 2015). Following superimposition, the mean shape and size were computed for each species, and mean sizes were log‐transformed. Shape disparity was computed as the average squared Procrustes distance of each species’ mean shape to the mean shape for all species, equivalent to the sum of variances over all superimposed coordinates (Zelditch et al. 2003; Zelditch et al. 2012).

MODEL SELECTION

We used a likelihood‐based approach to model selection, implemented in mvMORPH (version 1.1.4, Clavel et al. 2015). This approach finds the optimal parameter values for a model by maximizing the likelihood of observing the data under the model (Hansen 1997; Martins and Hansen 1997; Butler and King 2004). The Akaike information criterion (AIC) penalizes the likelihood by the number of parameters so that an improvement in fit is judged to be meaningful only if the increase in likelihood is large relative to the increase in the number of parameters (Butler and King 2004). The AIC corrected for small sample sizes (AICc) provides additional protection against overfitting when the sample size is small relative to the number of dimensions in the data (Hansen et al. 2008). Even with those protections against overfitting, likelihood‐based methods tend to favor complex models over simpler ones (Boettiger et al. 2012; Ho and Ane 2014; Cooper et al. 2016). This bias tends to increase with the dimensionality of the data (Adams and Collyer 2018).

Shape data are multivariate, and complex models can have many more parameters than the number of species in the study; however, using PCA to reduce dimensionality can reduce the bias favoring overly complex models. To avoid excessive reduction of dimensionality, which can introduce other biases (Uyeda et al. 2015), we reduced the data set to the seven PCs needed to describe at least 90% of the shape variance among species means. Each subsequent PC described less than 1.7% of the variation. In addition, the number of estimated parameters was reduced further by constraining the PCs to be adaptively independent (using decomp = “diagonal”), following a previous study of similar data (Zelditch et al. 2020), which found that the optimizer usually converged for constrained models but rarely did for unconstrained models. Each model was fit to the seven PCs as a single multivariate set.

To determine whether the difference in fit between the two models was meaningful, we implemented a parametric bootstrap approach outlined by Boettiger et al. (2012). The difference in fit of two competing models to the data was quantified as δ = −2(logL 0 − logL 1), where L 0 is the likelihood of the data under the simpler model and L 1 is the likelihood under the more complex model. Data were simulated 1000 times under each model using the rate parameters estimated by that model, and both models were fit to each simulation. This produced two distributions of δ, each representing expected differences in the fit of those models to the data if the data were generated by one of the models. The smaller the overlap between the two distributions, the greater the power of the analysis to distinguish between data sets generated under those models. Comparing the observed δ (the fit of the models to the observed data) to the distribution obtained for one set of simulations gives the likelihood of a difference in fit that large if the data had evolved under that model. If the observed δ fell within the distribution of values expected under the simpler model, that model could not be rejected in favor of the more complex one. Conversely, if the observed δ was higher than expected under the simple model, the more complex model was judged to fit the data significantly better. When two models had the same number of parameters, the same procedure was used to test whether the one with the higher likelihood was a significantly better fit. These analyses were performed using model simulation and fitting functions in mvMorph.

The procedures outlined above identify which model is the best fit to the data, but that model is not necessarily a good fit to the data. As a further check on how well the model fit the data, we evaluated how closely the disparity predicted by the model approximated the observed disparity, and we determined the probability of the model generating the observed disparity from the distribution of values across simulations under the model. For Ornstein‐Uhlenbeck (OU) models, we also calculated the rate coefficient, α, which can be interpreted as the strength of the pull to the optimum or translated into the time needed to move half the distance to the optimum, the phylogenetic half‐life, t ½ = ln(2)/α (Hansen 1997; Hansen et al. 2008). A small α means a weak pull to the optimum and a long time to reach it. When t ½ is large relative to the age of the lineage, there is a low probability of reaching the optimum, and therefore, the model is not a good explanation for the data. mvMORPH does not report the value of α for a multidimensional data set; however, the values for each individual dimension can be used to calculate the value for multidimensional data from the sums of the diagonals of the σ2 and stationary variance (vy) matrices: vy = σ2/2α (Zelditch et al. 2020).

MODELS

Three of the models evaluated in this study are the simple models commonly included in studies of the evolutionary dynamics of morphology: a single rate Brownian motion model (BM‐1), a single peak OU model (OU‐1), and the early burst model (EB). We also fit more complex models, including multiple‐peak OU models and corresponding multiple‐rate BM models, i.e., with rate changes rather than peak shifts mapped to the branches. Most of the complex models were based on four traits expected to predict adaptive divergence in mandible shape: diet, habitat, size, and tooth morphology. For each trait, we determined the number of distinct states predicting differences in how the mandibles are used to process foods. We hypothesize that these functional differences predict different optima for mandible shape due to their implications for the loading regime and other demands imposed on the bone. The criteria used to determine which species shared an adaptive regime are outlined below.

The Diet model (Fig. 3A) is based on the expectation that the material properties of foods influence feeding mechanics and thus the optimal shape of the mandible (Hylander 1979; Rosenberger 1992; Anapol and Lee 1994; Ross et al. 2012). The available data on ground squirrel diets (summarized in (Thorington et al. 2012)) suggest that they can be classified into two groups, granivores and folivores (McLean et al. 2018). Granivores have a diet composed largely of seeds, including nuts, along with a wide variety of other plant parts and some small animals; folivores have a diet that includes a large proportion of leaves from grasses and forbs and a smaller proportion of stiff, bulky items such as seeds, twigs, and underground storage organs. Although we have designated species with these diets as granivores and folivores to reflect the dominant foods in those diets, we acknowledge that ungulates with similar diets are usually characterized as intermediate or mixed feeders because they are between grazers that eat mainly grasses and browsers that eat both woody and herbaceous dicots (Hofmann and Stewart 1972; Cantalapiedra et al. 2014; Toljagic et al. 2018). Because leaves have several properties that make simple orthal crushing less effective than a chewing cycle that includes some horizontal shearing (Lucas and Luke 1984; Thiery et al. 2017), the lineages of folivorous ground squirrels are expected to shift to a novel mandibular shape. The phylogenetic distribution of ground squirrel diets suggests that shifts to folivory occurred in two lineages: one including only Marmota and the other including all descendants of the most recent common ancestor of Spermophilus and Cynomys (hereafter, the spermophiline clade). In the spermophiline clade, the branches leading to Ictidomys and Xerospermophilus are inferred to have diverged toward secondarily less folivorous diets; thus, there are four peak shifts.

Figure 3

Models based on a single ecological or morphological trait. In all cases, the phylogeny is shown with the same branching order as in Figure 2, and branches are color‐coded to indicate which lineages are on the same peak. (A) Diet: black – granivore, green – folivore; (B) Habitat: black – woodlands, green – grasslands, magenta – arid scrub and desert; (C) Size: lavender – small, black – medium, red – large, purple – extra large; and (D) Tooth morphology: black – Sciurus‐like, purple – Marmota, blue – Spermophilus‐like; red – Cynomys. Branches predicted to have peak shifts are labeled A, Ammospermophilus; C, Cynomys; I, Ictidomys; L, Callospermophilus; M, Marmota; O, Otospermophilus; S, spermophilines; T, Tamias; U, Urocitellus townsendii species group; and X, Xerospermophilus.

The Habitat model (Fig. 3B) is based on the expectation that changes in habitat will predict changes in diet and therefore changes in mandibular shape. The geographic distributions of ground squirrels suggest strong influences of cover and precipitation that tend to segregate them into three main habitat types: (a) woodlands, (b) grasslands, and (c) arid scrub and desert (Thorington et al. 2012). These habitats reflect a gradient of increasing aridity that predicts the relative abundances of fruits, seeds, and other foods; they also predict changes in the properties of any particular type of food (Jarman 1974; Hartley and DeGabriel 2016). These changes in food properties may not affect the magnitude of the peak force that must be applied, but instead increase the need to sustain the effective forces over a longer period or over a larger range of gape angles (Hylander 1979; Offermans and De Vree 1990; Williams et al. 2007). In our habitat model, there were three adaptive peaks corresponding to the three habitats and six peak shifts inferred from the phylogenetic distribution of habitat occupation. As in the diet model, shifts to the grassland habitat occurred in the Marmota and spermophiline lineages. Ictidomys and Xerospermophilus shifted from grasslands to more arid habitats, and two lineages (Ammospermophilus and Otospermophilus) appear to have made the transition to arid habitats directly from woodlands.

The Size model (Fig. 3C) reflects the premise that the optimal shape is a function of size. Under isometric scaling, lever arm ratios would be constant, but larger animals would have relatively stronger muscle forces because muscle force scales with cross‐sectional area. Therefore, animals evolving to larger sizes need to deviate from isometry by changing mandible shape to bear larger loads or by reducing muscle sizes and the corresponding muscle attachment areas to reduce the load. In contrast, animals evolving to smaller sizes would need to increase relative muscle size or compensate for relatively small muscles by changing muscle position to increase mechanical advantage or changing tooth shape to concentrate bite force in a smaller area. Changes in size also have implications for the gape angle required to eat a particular food item, which imposes additional changes in loading regime and mandible shape. Based on a prior optimization of size on the phylogeny (Zelditch et al. 2017), we identified four size peaks and eight peak shifts: two shifts to larger sizes and six shifts to smaller sizes. The size increases occurred on the branches leading to Marmota and Cynomys, producing species in two different ranges of sizes, occupying two different size peaks. One of the size decreases occurred on the branch leading to Tamias. Several other branches have substantial size reductions to size ranges that overlap the larger Tamias species; accordingly, all of these smaller species are inferred to occupy the same peak.

The Tooth model (Fig. 3D) reflects the premise that differences in tooth morphology are related to the material properties of the foods eaten (Lucas and Luke 1984; Strait 1993; Gailer et al. 2016) and the patterns of jaw movement in the chewing cycle (Janis 1979; Rensberger et al. 1984). Consequently, differences in tooth morphology predict differences in mandibular loading and therefore differences in mandibular shape. Based on previous descriptions of molar and premolar shapes (Bryant 1945; Goodwin 2009), we identified four distinct dentitions and three peak shifts. The primitive dentition present in the most basal extant lineages is very similar to that of tree squirrels (Sciurus‐like), which have relatively square teeth with low blunt cusps. In the spermophiline lineage, molars have much higher relief and are somewhat compressed antero‐posteriorly. These traits are developed to an even greater extent in Cynomys, which also has enlarged teeth at each end of the tooth row and other dental traits that separate it from other spermophilines. Marmota diverges from the Sciurus‐like dentition in a different direction, with greater mediolateral widening and a smaller increase in the relief of the crown morphology.

The factors underlying the four models in Figure 3 may interact, and the adaptive regimes may be predicted by a combination of traits better than by any single trait. For example, the loading regime imposed on a mandible while feeding on a particular diet may depend on the form of the teeth. Consequently, we formulated a model based on the intersection of the sets of taxa defined by the Tooth and Diet models, which produced a model (T+D) with five peaks: three for the three tooth morphologies in folivores (Marmota, Cynomys, and most spermophilines), and two for the two tooth morphologies found in granivores (one for all species with Sciurus‐like dentitions and the other for two genera with spermophiline dentitions, Ictidomys and Xerospermophilus). Similarly, we formulated compound models based on all other combinations of the four simple models (Table 2). Some possible combinations of models were rendered redundant by the nesting of Diet within Habitat, which occurred because all folivores live in grassland habitats; thus, D+H was equivalent to the Habitat model.

Table 2

List of models generated by combining single traits and their acronyms

Model	Acronym
Tooth+Diet	T+D
Tooth+Habitat	T+H
Tooth+Size	T+S
Diet+Habitat	D+H
Diet+Size	D+S
Habitat+Size	H+S
Tooth+Diet+Habitat	T+D+H
Tooth+Diet+Size	T+D+S
Tooth+Habitat+Size	T+H+S
Diet+Habitat+Size	D+H+S
Tooth+Diet+Habitat+Size	T+D+H+S

Because information about diet and ecology is incomplete and because theories of jaw function may not be the best model of jaw shape evolution, we used an expectation–maximization (EM) algorithm implemented in phylogeneticEM (version 1.4.0, (Bastide et al. 2017)) to identify potential shift locations that may indicate unanticipated functional transitions or combinations of them. The EM procedure alternates between estimating trait values at the ancestral nodes given a hypothesis of the model parameters and solving for the maximum likelihood estimate of the model parameters given the nodal values, iterating until convergence is achieved. Unlike mvMorph, phylogeneticEM treats the rate matrix as a scalar matrix, not as a more generalized diagonal matrix, and excludes the possibility of convergent evolution. The four best‐fitting models found by phylogeneticEM under its criteria (Fig. 4) were treated as additional hypotheses to be re‐evaluated in mvMorph using the same criteria as those used for the other models.

Figure 4

Models produced by phylogenetic expectation maximization, ordered by the number of peaks from four (A) to seven (D). The phylogeny is shown with the same branching order as in Figure 2, and branches are color‐coded to indicate which lineages are on the same peak. Branches predicted to have peak shifts are: Star, Marmotina; A, Ammospermophilus; C, Cynomys; M, Marmota; S, spermophilines; X, Xerospermophilus.

All four models generated by phylogeneticEM are very similar to the Tooth model; in fact, the five‐peak model (PEM5) differs from the Tooth model by the addition of a single peak (Fig. 5A and B). That addition is a novel peak for basal lineages of Marmotina, separating them from lineages that also have a predominantly granivorous diet and Sciurus‐like dentitions (i.e., Tamias and Sciurotamias). PEM5 is also similar to another 5‐peak model, T+D (Fig. 5C), which combines the Tooth and Diet models and adds a novel peak to the Tooth model for Ictidomys and Xerospermophilus, taxa with spermophiline teeth but diets that are more granivorous (less folivorous). In light of the similarity of these models, we also included an ad hoc 5‐peak model, PEM5+D (Fig. 5D), which also separated the relatively granivorous Ictidomys and Xerospermophilus from the more folivorous spermophilines but placed them on the same peak as the predominantly granivorous basal lineages of Marmotina.

Figure 5

Similarity of the PEM5 model to others in this study: (A) PEM5, (B) Tooth, (C) T+D, (D) PEM5+D. The phylogeny is shown with the same branching order as in Figure 2, and branches are color coded to indicate which lineages are on the same peak. All models share peak shifts at M – Marmota, S – spermophilines, and C – Cynomys. The star indicates a shift in basal Marmotina (PEM5 and PEM5+D). Ictidomys (I) and Xerospermophilus (X) occupy a unique peak in T+D; they share a peak with basal Marmotina in PEM5+D.

We also evaluated two models with only a single regime shift: one with the shift on the branch to Marmota and the other with the shift at the base of the subtribe Marmotina. These models are motivated by the large number of previous studies attributing large shifts in ecology and adaptive regimes to these branches (Bryant 1945; Black 1963; Casanovas‐Vilar and van Dam 2013; Ge et al. 2014; Zelditch et al. 2015; McLean et al. 2018). We added these models to evaluate whether either fit the distribution of mandible shape significantly better than a model of no divergence (OU‐1) and to test whether more complex models were significantly better than these models of a single regime shift.

Results

The OU model with the lowest AICc was PEM5+D (Table 3). Models with more peaks and parameters did not fit as well. The best fitting model was a considerable improvement over OU‐1, which fit much better than BM‐1 or EB. BM‐1 fit worse than the other simple models, and EB was a much smaller improvement over BM‐1 than was OU‐1. All of the multirate BM models evaluated for this study had higher likelihoods than BM‐1, but several had worse AICc values (Table 4). The multirate BM models all fit the data much more poorly than OU‐1; in fact, all but Size fit more poorly than EB.

Table 3

Comparison of OU models, BM‐1 and EB, ordered by number of parameters, then AICc. dAICc is the difference in AICc from the model with the lowest value, PEM5+D (bold). Shaded rows have more parameters and higher AICc than the model with the lowest AICc

Model	Peaks	Parameters	Loglikelihood	AICc	dAICc
T+D+H+S	9	98	1709.627	−3177.705	0.989
H+S	8	91	1680.999	−3141.329	37.365
PEM7	7	84	1674.466	−3148.477	30.217
T+D+S	7	84	1659.416	−3118.377	60.317
T+S	6	77	1661.818	−3142.763	35.931
PEM6	6	77	1658.497	−3136.122	42.572
T+H	6	77	1647.090	−3113.308	65.386
D+S	6	77	1641.276	−3101.680	77.014
PEM5+D	5	70	1670.294	−3178.694	0.000
PEM5	5	70	1664.867	−3167.840	10.854
T+D	5	70	1644.926	−3127.958	50.736
PEM4	4	63	1644.181	−3144.871	33.823
Tooth	4	63	1633.380	−3123.268	55.426
Size	4	63	1631.939	−3120.385	58.309
Habitat	3	56	1615.577	−3105.513	73.181
Diet	2	49	1604.107	−3099.898	78.796
Marmota	2	49	1601.763	−3095.211	83.483
Marmotina	2	49	1596.791	−3085.266	93.428
OU‐1	1	42	1583.342	−3075.191	103.503
EB	–	36	1539.749	−3002.039	176.655
BM‐1	–	14	1507.953	−2986.362	192.332

Table 4

Comparison of multirate BM models, ordered by loglikelihood. Shaded rows have AICc > BM‐1. The best fitting model, in bold, fits only slightly better than EB

Model	No. rate regimes	No. parameters	Log‐likelihood	AICc
H+S	8	63	1563.814	−2984.136
PEM7	7	56	1559.837	−2994.034
T+D+H+S	9	70	1559.659	−2957.424
T+D+S	7	56	1557.755	−2989.868
D+S	6	49	1554.471	−3000.626
PEM6	6	49	1552.237	−2996.159
T+S	6	49	1549.164	−2990.012
PEM5	5	42	1546.365	−3001.236
T+H	6	49	1545.686	−2983.056
PEM5+D	5	42	1545.667	−2999.840
T+D	5	42	1540.026	−2988.557
Size	4	35	1539.633	−3004.112
Tooth	4	35	1534.964	−2994.775
PEM4	4	35	1533.118	−2991.083
Habitat	3	28	1521.799	−2984.324
Marmotina	2	21	1519.049	−2994.261
Diet	2	21	1515.549	−2987.261
Marmota	2	21	1512.760	−2981.683

In addition to PEM5+D, we evaluated two models with five peaks: PEM5 and T+D (Table 3). All three models had separate peaks for Marmota and Cynomys and had Tamias on a separate peak from spermophilines (Fig 5), features they shared with the Tooth model (Fig. 3D). One difference between the 5‐peak models is whether Ictidomys and Xerospermophilus are on the same peak as other spermophilines, reflecting similar tooth morphologies (PEM5), or on a different peak reflecting their more granivorous diet (PEM5+D and T+D). The other difference is whether basal lineages of predominantly granivorous Marmotina are on the same peak as other taxa with Sciurus‐like dentitions (T+D) or on a different peak (PEM5 and PEM5+D). Although our current understanding of the phylogenetic relationships of Marmotini requires different numbers of evolutionary changes between peaks in these models, the models are equivalent with respect to the complexity of the adaptive landscape and the number of parameters to estimate.

The difference in AICc values between the Tooth and T+D models suggests that the latter fits the data slightly better (−4.69), but parametric bootstrapping produced broadly overlapping confidence intervals for the difference in fit of those models to simulations under the simpler tooth model (Fig. 6A). In fact, 29% of simulations under the Tooth model and 48.6% of simulations under the T+D model produced a δ as high as that observed. The results indicate that assigning the relatively granivorous spermophilines Ictidomys and Xerospermophilus to a separate peak does not significantly improve the fit to the data. The PEM5 model, which did not separate Ictidomys and Xerospermophilus from other spermophilines but did assign the predominantly granivorous basal lineages of Marmotina to a separate peak from other taxa with similar diets and dentitions (Tamias and Sciurotamias), did significantly improve the fit (Fig. 6B). The PEM5+D model, which placed more granivorous spermophilines (Ictidomys and Xerospermophilus) with more basal granivorous Marmotina, showed an even larger improvement (Fig. 6C). PEM5 and PEM5+D improved on T+D by similar margins (Figs. 6D and 6E), and again, the margin was larger for PEM5+D. The difference in fit between PEM5 and PEM5+D is not significant (Fig. 6F), even though a difference in AICc of −10.85 usually can be taken as indicating a very substantial improvement in analyses of simpler data. However, the frequency of δ as high as observed was 10.3% for PEM5 and 78.9% for PEM5+D, suggesting that the observed disparity is much more likely to have been produced by the model that placed the more granivorous spermophilines with more basal granivorous Marmotina (PEM5+D) than by the model that placed them with the other, more folivorous spermophilines (PEM5).

Figure 6

Analyses of relative fit by parametric bootstrap: (A) Tooth vs. T+D, (B) Tooth vs. PEM5, (C) Tooth vs. PEM5+D, (D) T+D vs. PEM5, (E) T+D vs. PEM5+D, (F) PEM5 vs. PEM5+D. Each frame shows two distributions of δ, representing differences in the fit of the models to simulations produced using the parameters of one model. The simpler tooth model (fewer parameters) is on the left in comparisons to it (A–C); in comparisons of models of equal complexity (same numbers of parameters), the model with higher AICc is on the left. Bars below the curves show the overlap of the 95% confidence intervals (Table 5), which reflects the power of the analysis to discriminate between the models. The heavy dashed line marks the observed δ for the difference in fit to the data; an observed δ greater than the upper limit of the confidence interval for the model with the lower range indicates that the model with the higher range is a significantly better fit to the data.

Table 5

Evolutionary rate (α) and the relationship between predicted and observed disparity (0.00399) for the models judged to be most informative about the dynamics of mandibular shape evolution

					Confidence Interval
Model	α	t ½(Ma)	Expected disparity (mean of simulations)	Proportion of observed disparity	Lower	Upper
OU‐1	0.035	19.85	0.00220	55.0%	0.00157	0.00364
Marmota	0.054	12.94	0.00251	62.9%	0.00174	0.00404
Diet	0.100	6.92	0.00318	79.6%	0.00216	0.00467
Teeth	0.138	5.01	0.00301	75.4%	0.00202	0.00473
T+D	0.145	4.77	0.00304	76.2%	0.00212	0.00440
PEM5	0.482	1.44	0.00383	96.0%	0.00333	0.00437
PEM5+D	0.393	1.76	0.00388	97.1%	0.00341	0.00434

The two best‐fitting models (PEM5 and PEM5+D) are very similar to the Tooth model, one of the single trait models. Comparison of the Tooth model to the Diet model, the simplest of the single trait models and the one that is most similar to the Tooth model, indicates that the Tooth model is a significantly better fit (Fig. 7A). Because the Tooth model is unique among the single trait models in having divergent shifts on the branches to three groups of folivores (Marmota, the spermophilines, and Cynomys), this result demonstrates that the divergence of those three branches is an important component of the dynamics of mandible shape evolution in ground squirrels. Comparison of the Diet model to OU‐1 (Fig. 7B) demonstrates that divergence of folivores from granivores, even without divergence among folivores, greatly improves on a model of no divergence. It is less clear that divergence of Marmota, alone, fit much better than OU‐1 (Fig. 7C). Although the difference in fit of the models to the simulations had broadly overlapping distributions, the data are considerably more likely under the Marmota model.

Figure 7

Analyses of the relative fit of the Diet and OU‐1 models to selected others: A) Diet vs. Tooth, B) OU‐1 vs. Diet, C) OU‐1 vs. Marmota. Each panel is organized as in Figure 6.

The poor fit of the OU‐1 model to the mandible shape data is underscored by its low value for α, which represents a phylogenetic half‐life (t ½) that is approximately 75% of the age of the most recent common ancestor of Marmotini (Table 5). Not only is the inferred evolutionary rate implausibly slow, but the disparity of mandible shape predicted by the OU‐1 model is approximately half of the observed value and is not encompassed by the 95% confidence interval computed from the simulations. The Marmota model is a relatively small improvement that predicts disparity less than 2/3 of the observed value and a t ½ that is approximately 50% of the age of the common ancestor. The Diet, Tooth, and T+D models are more substantial improvements with little to distinguish one from the others. All three have α values three to four times that of OU‐1, which still represents a rather slow approach to the optimum (t½ = 4.8 – 6.9 Ma). The predicted disparities are improved as well but are still < 80% of the observed value. The two best‐fitting models (PEM5 and PEM5+D) predict disparities that are > 95% of the observed value and t½ < 2 Ma.

To further investigate the differences between the two best‐fitting models, we examined the values of θ estimated in the simulations of those models that were used in the parametric bootstrap analyses. Those values were reported as scores on the PCs of the original data, allowing comparison across models and interpretation of the scores as shape differences in the original space of the data. In all simulations of the better fitting model (PEM5+D), the peaks for the three groups of folivores were widely separated from each other as well as from the two groups of granivores, but in 33.5% of the simulations, the peaks for the folivores were at distinctly larger values on PC1 (Fig. 8A). For each group of folivores (Marmota, Cynomys, and the spermophilines), the difference in scores between the two sets of results was so great that the ranges of values did not overlap. The larger values on PC1 were associated with larger values on PC2 for Marmota (Fig. 8A) and on PC3 for Cynomys (Fig. 8B). The folivore peaks differed from the granivore peaks primarily in having a relatively expanded angular process and relatively short coronoid and condyloid processes (Fig. 8C), with Marmota and Cynomys diverging more than the folivorous spermophilines. In both sets of results, the position of Marmota on PCs 2 and 3 indicates its optimal shape differs from basal Marmotina in a different direction than that of the spermophilines, including a relatively shorter and broader condyloid process, a less reduced coronoid process and a deeper angular process. In contrast, the optimal shape for Cynomys has a more elongated angular process with a relatively less shortened condyloid and a more shortened and straightened coronoid.

Figure 8

Locations of model optima (θ) on PCs of the species means for 1000 simulations of the two best supported models. A) Locations of optima on PCs 1 and 2 for simulations of PEM5+D, the best supported model. The upper panel shows a subset of 665 results that have optima at relatively low values of PC1 for Marmota and Cynomys; the lower panel shows 335 simulations with values for Marmota and Cynomys outside the range observed in the upper panel. Ellipses encompass 95%, 75%, and 50% of the cases in each cluster. Blue, Tamias, and Sciurotamias (granivores not in Marmotina); red, basal Marmotina; green, Marmota; purple, spermophilines (including granivores Ictidomys and Xerospermophilus); and orange, Cynomys. (B) Locations of optima for PEM5+D on PCs 1 and 3, with upper and lower panels as in (A). (C) Shapes representing the largest positive scores on each axis are shown as deformations of the mean shape. (D) Locations of optima on PCs 1 and 2 for simulations of PEM5, the second‐best supported model. The upper panel shows a subset of 739 results that have optima at relatively low values of PC1 for Marmota and Cynomys; the lower panel shows 261 simulations with values for Marmota and Cynomys outside the range observed in the upper panel. Ellipses encompass 95%, 75%, and 50% of the cases in each cluster. Blue, Tamias, and Sciurotamias (granivores not in Marmotina); red, granivorous Marmotina, including Ictidomys and Xerospermophilus; green, Marmota; purple, folivorous spermophilines; and orange, Cynomys. (E) Locations of optima for PEM5 on PCs 1 and 3, with upper and lower panels as in (D).

In PEM5, the granivorous spermophilines, Ictidomys, and Xerospermophilus were assumed to be on the same peak as the folivorous spermophilines, not on the peak with basal Marmotina as in the other model. The simulations of this model also produced two distinct sets of results (Fig. 8D and E). For this model, the peaks for spermophilines, Marmota, and Cynomys were at high PC1 values in 26.1% of the results. The assignment of Ictidomys and Xerospermophilus to the same peak as other spermophilines rather than the peak with other granivorous Marmotina had little effect on the relative positions of those peaks. The similarity of the peak positions between this model and PEM5+D explains why the difference in fit to the data was not meaningful and was not associated with a substantial difference in expected disparity or evolutionary rate.

The differences between the optimal shapes and the functional implications of those differences may be more easily appreciated by comparing the evolutionary transformations that can be inferred from them (Fig. 9). Although Figure 8 demonstrates that both models would entail similar changes in the direction of shape evolution, we used the optima for PEM5+D because it is the better fitting model of two equally complex models. A somewhat more consequential issue is the inference of the ancestral peak. Because Tamias and Sciurotamias are on different branches from the root node, the peak occupied by them could be interpreted as the ancestral peak, in which case, the first peak shift in this lineage is from that peak to the one occupied by basal Marmotina (Fig. 9A). This shift would have involved a relative shortening of the condyloid process and a relative expansion of the upper portion of the angular process. An alternative that should be considered given uncertainties about the position of Sciurotamias within Marmotini (Fabre et al. 2012; Ge et al. 2014; Sinitsa et al. 2021) is that Sciurotamias and Tamias converged on the peak they occupy from the peak occupied by basal Marmotina, increasing the relative length of the coronoid and reducing the relative size of the upper angular process. The shape transformation from Tamias to basal Marmotina would have expanded the masseter attachment area and shortened output lever arms for a relatively stronger masseter‐driven bite, with the trade‐off of a smaller gape and reduced contribution to bite force from the temporalis muscle. The opposite transformation would have produced relatively elongated jaws with shorter input arms allowing a larger gape but producing a relatively weaker bite. This combination would increase cheek pouch capacity, compensating for a weaker incisor bite with the ability to transport more seeds to the safety of a burrow; it would also have the advantage of a larger gape and faster closing, which would be beneficial for more predaceous species.

Figure 9

Shape changes between optima. (A) Ancestral peak to basal Marmotina, (B) basal Marmotina to basal spermophilines, (C) basal Marmotina to Marmota, (D) basal spermophilines to Cynomys.

From basal Marmotina, there were two different transitions associated with increased folivory, one in the divergence of spermophilines and the other leading to Marmota. The shape changes in the transition to the spermophiline peak (Fig. 9B) involved a different combination of local changes than were involved in the divergence of the basal Marmotina (Fig. 9A). There were similar transformations of the angular process, but in spermophilines, it was associated with a relative shortening and posterior shift of the coronoid process rather than a relative shortening of the condyloid process. These changes would have reduced input lever arms without a coincident shortening of the output lever arms and reduction of gape, alleviating a trade‐off in the relationship between input and output lever arms. The shift to the Marmota peak was a much larger shape change, and every region of the mandible was transformed (Fig. 9C). The condyloid process became relatively shorter, and the angular process became relatively deeper as well as longer. In addition, the coronoid process was straightened, and the region below the molars became relatively deeper. These shape changes, which are associated with a large increase in body size, are consistent with increasing relative bite force and with increasing ability to resist deformation under loading.

The shift from the spermophiline peak to that occupied by Cynomys was also associated with a size increase, although a smaller one than in Marmota. However, the direction of shape change in Cynomys (Fig. 9D) was more similar to the divergence of the spermophilines. Unique to the divergence of Cynomys was rotation of the condyloid process that raised the condyle above the tooth row, increasing the input arm for the muscles attached to the larger angular process while reducing the input arm for the muscles attached to the smaller coronoid process. These changes in the posterior processes were associated with a relative shortening of the diastema as in Marmota and similarly augmented incisor bite force beyond that due to increased body size.

In summary, each of the transformations of mandibular shape involved a unique modification of the shape of the angular process, and the transformation of the angular process was combined with changes in different sets of other local regions. In all cases, at least one other posterior process was involved. In some cases, one or both components of the horizontal ramus (diastema and molar region) were involved. The ability to independently evolve changes in local regions of the mandible enabled the lineages within Marmotini to meet different combinations of potentially conflicting mechanical demands imposed by new foods and new environments.

Discussion

Our results demonstrate that the adaptive landscape for mandibular shape in ground squirrels is more complex than predicted by a simple model of a single dietary shift. A two‐peak model with folivores and granivores on different peaks does fit better than a single peak model, but a five‐peak model with both folivores and granivores occupying multiple peaks fits even better. The two granivore peaks were not predicted by any of our a priori models. On one of those peaks, the species are heterogeneous in all the traits used to formulate those models. On the other peak, almost all species are miniatures (Tamias), which may indicate that the evolution of body size was an important factor in the evolution of mandibular shape in that lineage, as has been demonstrated in other taxa (Anapol and Lee 1994; Taylor 2002; Monteiro and Nogueira 2009; Zelditch et al. 2017). Sharing the peak with Tamias is the much larger Sciurotamias, which may be eating foods that are large for its size, or it may be convergent with Tamias in other important aspects of jaw function. The presence of three optima for folivores is more readily explained, in large part that is because each peak is associated with a unique tooth morphology that may represent a different solution to the challenges of eating leaves as well as a difference in the types or proportions of leaves eaten. Our results for folivorous ground squirrels are consistent with previous analyses showing that divergence among folivores reflects not only the different properties of the leaves that are eaten but also the properties of other foods that are included in the diet (Taylor 2002; Cantalapiedra et al. 2014; Wang et al. 2021). In addition, two of these optima are associated with shifts in body size, providing further demonstrations that the optimal shape is partly a function of size. When all of the adaptive peaks for folivores and granivores are considered, it is clear that the diversification of ground squirrels was more complex than a single shift to a novel dietary resource.

The adaptive landscape may be even more complex than the best fitting model. Although the predicted disparity is close to the observed value, the value of α equates to a phylogenetic half‐life that is several times longer than is predicted from rates of adaptive evolution within populations (Lande 1985, 1986; Lynch 1990; Hansen 2012). This result could be due to low or sporadic rates of shape evolution early in the transition, reflecting the pace of environmental changes (Toljagic et al. 2018). Climate data (Zachos et al. 2001) suggest that slow rates of net environmental change could explain the low evolutionary rates of lineages that diverged early in the Miocene, as did some older branches of granivores in the basal Marmotina, but slow evolutionary rates of lineages diversifying in the later Miocene and Pliocene (both granivores and folivores) may reflect rapid fluctuations in climate during this period of more rapid net change. A previous analysis of evolutionary rates in all Sciuridae showed considerable differences in rate among monophyletic groups within Marmotini, including striking differences among lineages of similar age (Zelditch et al. 2015). It is possible, therefore, that the value of α in this study reflects the combination of lineages with several different rates, some that may be ascribed to patterns of climate or other environmental change.

Our results for ground squirrels are strikingly different from the results of a similar study on tree squirrels (Zelditch et al. 2020). The diets of tree squirrels in the core lineages differed mainly in the replacement of hard, bulky foods such as nuts with softer seeds and fruits, and that study found that the distribution of their mandibular shapes was best fit by a single peak model. Traits that enhance the ability of tree squirrels to bite through harder or tougher shells do not substantially reduce the ability to bite through less resistant husks or the ability to crush weaker endosperm or fruit pulp with the molars. Thus, tree squirrels can occupy a single peak because the mechanical trade‐offs of increasing incisor bite force do not have costs with respect to foraging success. In contrast, the more challenging foods eaten by ground squirrels are leaves of grasses and forbs, which are tougher, thinner, and more flexible than nuts and seeds, but not stronger. Furthermore, the traits that enhance the ability to eat leaves reduce the ability to eat other foods: broader and flatter incisors that crop more leaves or stems per bite also distribute muscle force over a larger area, reducing the ability to bite harder foods; muscle orientations that improve shearing between the molars mean that a smaller proportion of the muscle force contributes to compression. Because there are different ways to change the positions of muscles and the arrangements and shapes of teeth, there are many possible ways of balancing the trade‐offs imposed by the conflicting mechanical demands of different foods, which makes the adaptive landscape complex.

The complexity of the musculature and dentition may explain why there is more than one adaptive peak possible for ground squirrels, but it does not necessarily follow that all possible peaks will be occupied. It has been proposed that different ancestral traits might make different adaptive peaks accessible, leading to “imperfect convergence” (Collar et al. 2014). The greater the complexity of the system, the more ways there might be for different ancestors to evolve toward the same functional abilities, and if those alternatives are functionally equivalent with respect to a critical measure of the performance of that system, there is no reason for those ancestors to converge on the same peak. However, an alternative explanation for related taxa occupying adaptive peaks that are similar but not identical is that those taxa are not convergent but divergent. They may have undergone generally similar shifts, such as including tougher or harder foods, but those foods are also different with respect to other properties that affect how they are handled. Similar patterns of overlapping diets have been described in other taxa (Hofmann and Stewart 1972; Ferrarezzi and do Amaral Gimenez 1996; Norconk and Veres 2011). These differences in diet lead to different optimal morphologies because the trade‐offs are weighted differently. In addition, different species may handle the same food differently as a consequence of their unique morphological features, reinforcing their divergent requirements for the optimal mandibular morphology. Instead of a many‐to‐one relationship between form and function in which different morphologies are equally capable, the relationship is many‐to‐many because each morphology must process many foods and many morphologies are used to process each food.

A question that remains is why and how do taxa with broadly overlapping diets maintain those relationships rather than specializing more narrowly. The answer may lie in the relative stability afforded by having the ability to use a variety of resources. Modeling studies suggest that the ability to be flexible in resource use may be the key to persistence for many taxa, especially those classically recognized as omnivores (Ingram et al. 2009, Ingram et al. 2012). Optimal foraging theory provides an additional rationale for balancing the ability to use challenging resources with the ability to continue using less challenging ‘preferred’ foods when opportunity permits (Robinson and Wilson 1998). This versatile specialist strategy may explain the 30 Ma years of stasis in some tree squirrel lineages (Zelditch et al. 2020). The underlying principle may also apply to taxa that are less specialized. Studies on African bovids suggest that their overlapping diets permit coexistence in heterogeneous environments through a mix of patch and resource partitioning (Hofmann and Stewart 1972; Cantalapiedra et al. 2014), with corresponding differences in mandibular morphology (Wang et al. 2021). Similar dynamics may explain how it is possible for 8–10 species of ground squirrels to occupy the western Great Plains and as many as 5 species to live in close proximity within that region (Thorington et al. 2012).

Dietary divergence may provide the impetus for mandibular shape to diverge, but functional differentiation of the dentition and adductor muscles is the mechanical foundation for differences in loading regime and optimal shape. The ability of the bone to respond to changes in the loading regime depends on the ability of local anatomical regions within the bone to respond individually to local changes in the loading regime while maintaining the structural integrity of the whole bone. Different functional transitions will require different combinations of local changes, thus producing different patterns of evolutionary integration. The analyses presented in this study were not designed to test hypotheses of integration because they did not directly assess hypotheses concerning the coordination of changes in different regions (cf. Zelditch and Goswami 2021); however, the results shown here do suggest models that could be tested in future studies. The transitions inferred from our best‐fitting model do not appear to be consistent with the two‐part model inferred from quantitative genetic analyses of lab mice (Bailey 1985; Leamy 1993; Klingenberg et al. 2001). Our results suggest that the components of the posterior region (the posterior processes) will not be well integrated and that there may be more integration between sections of different regions (i.e., the diastema and the angular process) than within regions. These results suggest patterns of integration similar to models from developmental biology (Atchley and Hall 1991) and the patterns of intraspecific variation in wild‐caught deer mice and fox squirrels (Zelditch et al. 2008, Zelditch et al. 2009). A similar variety of combinations of local changes was also seen in mandibles of bats that diverged in diet (Monteiro and Nogueira 2009). Such a variety of evolutionary patterns appears to support the contention that for a lineage to diverge, it must have the flexibility to change the genetic integration of those regions to allow for their divergent evolutionary responses (Arnold et al. 2001; Monteiro et al. 2005; Marroig et al. 2009).

Decoupling makes a larger morphospace available, allowing access to a greater diversity of morphological responses to evolutionary pressures (Liem 1973; Vermeij 1973; Lauder 1982). However, the diversity of challenges that morphological evolution must meet is a product of mechanical and functional trade‐offs. The trade‐offs comprise the impetus to diverge because it is not possible to meet both demands simultaneously. They also comprise the impetus to reorganize an unsuitable pattern of integration into a more modular pattern that allows one trade‐off to be exchanged for another that is less costly. Thus, trade‐offs are not an obstacle to change but a driver of innovation.

AUTHOR CONTRIBUTIONS

D.L.S. designed the study and drafted the manuscript with input from M.L.Z. D.L.S. and M.L.Z. collected the data. D.L.S. performed the analyses.

DATA ARCHIVING

Data are available on Dryad https://doi.org/10.5061/dryad.kq1g6.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Associate Editor: J. Clavel

Handling Editor: T. Chapman