Ecological correlates of cranial evolution in the megaradiation of dipsadine snakes
BMC Ecology and Evolution volume 23, Article number: 48 (2023)
Dipsadine snakes represent one of the most spectacular vertebrate radiations that have occurred in any continental setting, with over 800 species in South and Central America. Their species richness is paralleled by stunning ecological diversity, ranging from arboreal snail-eating and aquatic eel-eating specialists to terrestrial generalists. Despite the ecological importance of this clade, little is known about the extent to which ecological specialization shapes broader patterns of phenotypic diversity within the group. Here, we test how habitat use and diet have influenced morphological diversification in skull shape across 160 dipsadine species using micro-CT and 3-D geometric morphometrics, and we use a phylogenetic comparative approach to test the contributions of habitat use and diet composition to variation in skull shape among species.
We demonstrate that while both habitat use and diet are significant predictors of shape in many regions of the skull, habitat use significantly predicts shape in a greater number of skull regions when compared to diet. We also find that across ecological groupings, fossorial and aquatic behaviors result in the strongest deviations in morphospace for several skull regions. We use simulations to address the robustness of our results and describe statistical anomalies that can arise from the application of phylogenetic generalized least squares to complex shape data.
Both habitat and dietary ecology are significantly correlated with skull shape in dipsadines; the strongest relationships involved skull shape in snakes with aquatic and fossorial lifestyles. This association between skull morphology and multiple ecological axes is consistent with a classic model of adaptive radiation and suggests that ecological factors were an important component in driving morphological diversification in the dipsadine megaradiation.
Understanding the dynamics and causes of ecological and morphological diversification during major radiations is a key challenge in evolutionary biology. The relationship between ecology and morphology is central to our ability to test the role of ecological opportunity and other factors in mediating lineage and phenotype diversification during such radiations [1,2,3]. However, the ecological basis of phenotypic variation – while readily identifiable with relatively simple phenotypic traits  – presents an acute challenge when considering complex, highly-integrated and/or modular morphological structures [5, 6] with multiple functional modalities, such as the vertebrate skull. Previous work on vertebrate skulls has revealed an unexpected complexity to the form-function-ecology relationship [7,8,9]. Obtaining a clear picture of how this highly intricate structure diversifies in parallel with ecological factors has the potential to illuminate how and why some vertebrate clades have become so much more diverse than others, even when those other lineages appear to have had similar ecological opportunity and biogeographic context.
In this article, we describe the evolutionary dynamics of the skull during one of the most diverse continental vertebrate radiations: the dipsadine snakes, which have undergone an extraordinary evolutionary explosion in the neotropics , diversifying into at least 806 species [11, 12] and a wide variety of ecological roles. In many cases, dipsadines account for over 50% of the species richness within neotropical snake communities [10, 13,14,15]. The (new world or western hemisphere) dipsadine megaradiation has received relatively little attention from a macroevolutionary perspective, despite its extreme ecological and morphological diversity (Figs. 1 and 2).
Dipsadines exhibit a vast range of dietary and habitat use profiles: diet can include snails , slugs , arthropods , annelids , bird eggs , mammals , birds , snakes and lizards , and fish  among other items, and individual species are characterized by varying degrees of specialization on those resources. Habitat use is also extremely variable in this clade, ranging from subterranean  to surface , arboreal , aquatic  and intermediate microhabitat use [29, 30]. Dipsadines encompass much of the ecological diversity seen across all snakes, and they have independently evolved many specialized ecologies that occur in distantly-related and biogeographically-separate snake lineages. Even within dipsadines, several specialized ecologies appear to have evolved multiple times (e.g. aquatic habitat use, molluscivory; Fig. 2). The recurrence of similar ecologies across the dipsadines suggests that they are a good system in which to test whether parallel phenotypic evolution has occurred in response to common selection pressures, as predicted by the ecological theory of adaptive radiation .
We focus on skull evolution because skull shape has important functional consequences for numerous ecological axes in vertebrates, including prey acquisition, handling, and ingestion [31, 32] habitat use , locomotion , sexual competition [35, 36], mate choice , and defense . Numerous studies have investigated the dynamics of skull evolution with respect to ecology in various taxa, including mammals [7, 9, 35, 39, 40], birds [8, 41], ray-finned fishes [42, 43], and lizards [44, 45]. In general, these studies have shown a significant link between ecology and shape, although the importance of ecology in explaining morphological variation appears to vary across specific factors (e.g. diet, habitat, etc.) as well as across clades.
Numerous functional innovations distinguish the snake skull from that of their limbed (and limbless) “lizard” relatives, as well as from the skulls of most other vertebrates [46,47,48,49,50]. Perhaps most significantly, snake skulls are characterized by highly mobile elements (kinesis), partly accounting for their ability to consume large prey relative to their body size (Fig. 3) [38, 49, 51]. Snakes are thus able to exploit unique dietary niches that more gape-limited taxa are unable to capitalize on [38, 46]. These innovations may have provided snakes with a source of intrinsic ecological opportunity that allowed them to successfully radiate in novel environments [1, 3].
Many intriguing examples of convergent morphologies have been described within particular snake subclades , suggesting that skull morphology in snakes is, at least in some contexts, both evolutionarily labile and responsive to ecological pressures. One of the most iconic examples is found in snail eating specialists, namely, the new world dipsadine genera of Dipsas and Sibon and the old world genus Pareas. These geographically and taxonomically disparate groups are so similar in morphology that they were once thought to make up a single subfamily, although molecular and morphological data now show that they are separated by nearly 45 million years of independent evolution [12, 52]. Subsequent studies have corroborated the impact of diet on the evolution of snake head  and skull morphology [54,55,56,57,58,59,60] as well as the impact of habitat use [53, 59, 61,62,63,64,65,66], foraging mode and locomotion  and prey subjugation mode  in driving morphological diversification in the snake skull and head.
In our analyses of dipsadine skull evolution, we assess the breadth of morphological variation across the radiation in a phylogenetic comparative framework and we determine whether diet and habitat use are significant ecological predictors of this variation. A significant relationship between ecological traits and morphology would suggest that the extraordinary species richness of dipsadines has resulted in part from diversification along these ecological axes, thus suggesting a prominent role for adaptive radiation  in the assembly of neotropical snake faunas.
In the non-trophic module analyzed for all 160 species, the first two PC axes accounted for 46.3% of shape variation, with PC1 accounting for 27.2% and PC2 accounting for 19.1% (Fig. 4, Fig. S2). Low PC1 values were associated with very broad, short skulls while high PC1 values were associated with very narrow, elongate skulls. PC2 tracked the extent to which the anterior portion of the skull was expanded and the extent to which the “snout” was elongated. Low PC2 values were associated with narrow anterior portions of the skull and elongate snouts while high PC2 values were associated with a severely expanded anterior region and short snout (Fig. 4, Fig. S2). Descriptions of the PC axes for the remaining skull modules can be found in Additional file 4: Appendix S4 of the supporting material.
Description of skull morphospace
Within the non-trophic module, semi-fossorial snakes tend to cluster towards the right on PC1, indicating generally narrower skulls (Fig. 4A), with extreme examples of this morphology apparent in the genera Apostolepis, Phalotris, and Elapomorphus. Heterodon and Xenodon dorbignyi are the exceptions to this general pattern (labeled in Fig. 4); in contrast to all other semi-fossorial species, they exhibit extremely wide skulls with a short snout. Aquatic snakes, representing species from three phylogenetically disparate clades, cluster very closely together in morphospace, indicating highly similar morphologies. No clear patterns for other habitat use groups emerge, although there is some overlap between cryptozoic and semi-fossorial taxa. Patterns in the non-trophic module morphospace with respect to diet are less clear; species specializing on mollusks and annelids seem to deviate in morphology from other species, however, the morphological correlates involved are not clear. Fish eating species also cluster to some extent (Fig. 4B, Fig. S2). A detailed account of the results we obtained from morphospace analyses for each of the remaining modules can be found in the supporting material (Additional file 4: Appendix S4).
When considering all modules, aquatic habitat use most often resulted in tight morphological clustering (as well as in deviation from other groups in the non-trophic and maxillary modules) and semi-fossorial habitat use most often resulted in deviations from other habitat use groups. In terms of diet, mollusk and annelid specialization most often resulted in a deviation from other morphologies.
For the majority of skull modules, PC1 and PC2 cumulatively explained more than 70% of shape variation (see discussion in Additional file 4: Appendix S4, Table S1), although this is at least partly a function of the number of landmarks considered per module (Table S1). For modules where this was not the case, such as the non-trophic module, variation in additional PC axes was qualitatively determined to either lack any clear trends across ecological groups or if trends were visible, they were redundant with trends exhibited by the first two PC axes. We therefore chose not to address variation in axes beyond PC2 in this manuscript.
PGLS and PPLS analyses
Our PGLS analyses corroborated the trends visually identified in morphospace (Table 1). When considering multivariate shape data with procD.pgls, habitat use was a statistically significant predictor of shape in all 8 modules (Table 1). Primary diet was a statistically significant predictor of multivariate shape in the maxilla, mandible, and palatine modules, although in these modules it had a smaller effect size than habitat use (Table 1). Size significantly predicted shape in all modules except the ectopterygoid, but generally had a smaller effect size than habitat use (Table 1).
Analyses of univariate shape data with the gls implementation of PGLS generated similar results when compared with multivariate shape analyses. Habitat use significantly predicted shape in all but the quadrate and palatine module when considering PC1 shape data and in all modules when specifying PC2 as the response variable. Primary diet significantly predicted shape (PC1) in the maxilla and pterygoid. For PC2, primary diet was a significant predictor of shape in the non-trophic module, maxilla, ectopterygoid, and palatine. Size significantly predicted shape in the non-trophic module, quadrate, and mandible for PC1 and the quadrate and mandible in PC2 shape data.
The phylogenetic partial least squares analysis using diet proportion data was largely consistent with the results obtained in the PGLS analyses. For the non-trophic module (r-PLS=0.6, P=0.01), maxilla (r-PLS=0.52, P=0.02), quadrate (r-PLS=0.45, P=0.048), mandible (r-PLS=0.69, P=0.001), pterygoid (r-PLS=0.49, P=0.01), and palatine (r-PLS=0.49, P=0.01), shape significantly covaried with diet. Consistent with most PGLS analyses, covariance of shape with diet in the ectopteryoid (r-PLS=0.37, P=0.14) and supratemporal module (r-PLS=0.35, P=0.25), were non-significant.
Robustness of PGLS results
We assessed the robustness of the results obtained from the PGLS analyses through the application of a combined approach utilizing both real and simulated data. We compared PGLS analyses using the true grouping configuration (i.e. shape data associated with true corresponding ecological data) to those from datasets with categorical groupings randomly permuted across the tips. Note that this procedure both (1) removes any associations between shape and ecology present in the true data, but also (2) removes any phylogenetic signal that may be present in the ecological data. We also tested the effect of sample size (number of species) on the robustness of the results by repeating this test using varying sized subsets of the original dataset. 1000 repetitions of this simple permutation test were run for each of 25, 50, 75, 100, and 125 species subsets of the non-trophic module shape data using both procD.pgls and gls functions. We incorporated the same model as in our main PGLS results for the dataset, including size as a covariate.
For randomly permuted datasets, we find that p-values are skewed towards zero for both the habitat use and diet variable (Fig. 5 – center, Fig. S5 – center), thus indicating an increased Type I error rate for our analyses. This trend seems to become exacerbated as the data subset size increases. Note, however, that for the habitat variable, which significantly explained shape in the original analysis of the non-trophic module with procD.pgls (P=0.001), p-values across all real dataset subsets (Fig. 5 – top left, Fig. S5 – top left) are substantially lower than those with randomly permuted ecological groupings (Fig. 5 – top center, Fig. S5 – top center).
To determine whether the elevated Type I error rates observed under permutations were a property of the underlying PGLS framework or of our data, we conducted a third test using simulated data. First, we estimated the covariance matrix for the real non-trophic module shape data using the ratematrix function in the R package “geiger v2.0.10” [67, 68]. We then used the sim.char function in the same R package to create a simulated dataset of n=148 “species”, containing a multivariate response variable with n=116 variables (equal to the number of principal components of the non-trophic module). We simulated this multivariate response variable under Brownian motion on the same phylogeny used throughout this study, utilizing the covariance structure of the true shape data as estimated previously. We paired this simulated shape dataset with the true ecological data, size data, and phylogeny and repeated the random permutation robustness test as before.
Our results show that when retaining all factors constant (ecological data, size data, phylogeny, shape covariance structure) but simulating the response variable under multivariate BM, PGLS analyses on the datasets behave precisely as expected, with a median p-value of ~0.5 recovered with both functions and across all data subset sizes for both the habitat use and diet variable (Fig. 5 – right, Fig. S5 – right). These results suggest that PGLS may be sensitive to non-Brownian evolution of the underlying morphological traits. It is possible that many complex morphological traits that have been assessed in a PGLS framework show similar departures from underlying assumptions – namely, Brownian motion – such that the resulting distribution of p-values is unreliable. Robustness tests similar to those performed here are rarely, if ever, included in PGLS analyses to assess relationships between ecological predictors and multivariate morphological response variables. We thus recommend that researchers utilizing PGLS implementations for complex shape data consider the potential sensitivity of these methods to deviations from null assumptions.
Recent phylogenetic evidence [12, 69] suggests that dipsadines originated in the Old World before dispersing to the New World relatively early in their evolutionary history, subsequently diversifying into the tremendous diversity of species seen today. The exceptional diversity of trophic ecologies and foraging modalities within the dipsadines suggests that the diversification of this group may have been facilitated by ecological opportunity [1, 3, 70]. It is possible that South and Central America were, at the time of dipsadine colonization, relatively depauperate in highly-derived caenophidian snakes with evolutionarily versatile prey subjugation and processing phenotypes [10, 60], thus setting the stage for the rapid evolution of novel trophic ecologies observed in the group . The apparently rapid lineage diversification in the dipsadines is also coupled with fast rates of trophic diversification and expansion of trophic space , further supporting the idea that the overall tempo and mode of this radiation reflects a response to ecological opportunity.
Nonetheless, it is overly simplistic to assume that dipsadines radiated in an ecological vacuum resulting from the simple absence of other snake lineages. Several major clades of "advanced" snakes were likely contemporaneous with dipsadines as they radiated , including colubrines , elapids, and viperids. Although the elapids in South America show substantially less diversity in morphology, diet, and habitat use than the dipsadines, related lineages of elapid snakes have undergone dramatic radiations in body form and ecology in other regions, most notably Australia [61, 73]. Hence, a key outstanding question in the assembly of New World reptile communities is whether dipsadine snakes radiated in response to extrinsic ecological opportunity, or whether they evolved novel phenotypes that facilitated subsequent ecological and lineage diversification even in the presence of potential competitor lineages (e.g., colubrine snakes; ). Several dipsadine lineages have acquired striking morphological  and histochemical [75, 76] adaptations that enable them to feed on unusual prey that are rarely consumed by other snake lineages, consistent with the hypothesis that the dipsadine bauplann is characterized by greater evolutionary versatility  than many other snake lineages.
We found that habitat and diet were significant predictors of shape in many dipsadine skull modules, across both multivariate and univariate shape data. Despite the statistical anomalies discussed in the above section, the p-values associated with ecological variables that significantly explained shape data in our analyses were substantially lower than with randomized data for the same modules/variables (Fig. 5, Fig. S5). Thus, our results suggest that ecological factors contribute to overall variation in skull shape across the dipsadine megaradiation. This is not to say that ecology is necessarily the most important factor in explaining variation in skull shape in dipsadines. Size was a significant predictor of shape across many modules, although it generally had a weaker signal than habitat use (Table 1). This is not an unexpected result, as many studies have shown that allometry often plays a role in explaining morphological disparity in skulls among vertebrates [8, 58, 78,79,80,81]. Although an in-depth discussion of allometric trends is outside the scope of this study, particularly given the dominance of ecological variables in explaining the data, we provide allometric regression plots in the supplemental results (Figs. S6, S7). We also cannot exclude the possibility that other untested factors such as pleiotropy, integration with other traits, or developmental constraints may play a more dominant role in explaining skull morphology than the variables tested here, and some studies have indeed shown that these factors can have a strong influence on skull evolution [6,7,8]. For instance, Bright and colleagues  found that in bird skulls, classically thought of as model examples of adaptation, the beak and braincase are highly integrated and developmentally constrained structures correlated with size, rather than independently evolving elements controlled by diet. Nonetheless, our results provide quantitative support for a correlation between ecology and morphology, long considered a central feature of adaptive radiation .
We showed that habitat use is a significant predictor of shape in a greater number of skull regions in comparison to diet and size, generally with higher effect sizes, indicating that habitat use may impose stronger selective pressures on the skull as a whole. Deepak et al.  recovered a similarly strong influence of habitat over diet in studies of natricine snake head morphology. In our study, the most distinct ecological groupings in morphospace included aquatic and semi-fossorial species, a finding that is consistent with Watanabe et al.  who recovered rampant convergence in skull shape within these ecological groups across squamates as a whole.
Savitzky  observed widespread convergence in aquatic/piscivorous snakes, noting that these species exhibit long quadrates in particular, which he hypothesized increased gape size for more rapid and effective consumption of difficult to handle fish prey. One study , corroborated this qualitative assessment with a morphometric study showing that among natricines, fish-eating species had relatively longer quadrates that appeared to significantly reduce prey-handling time. This finding was mirrored by Silva and colleagues  that found increased quadrate length in an aquatic coral snake relative to two terrestrial species. In our study, aquatic species also showed evidence of longer quadrates and tight clustering in morphological space with respect to quadrate shape (Fig. 6D1). We also found that aquatic snakes tightly clustered and tended toward an elongate skull in the analysis of the non-trophic module (Fig. 4A), a finding corroborated by several studies that speculated this to either function in increasing hydrodynamic efficiency or handling of fish prey [53, 55, 63, 83]. In cetaceans and crocodilians, aquatic clades of vertebrates that have both evolved varying degrees of elongation in the skull, McCurry and colleagues  found that extreme examples of elongation (i.e. gharials and river dolphins) were convergently evolved as a result of diet rather than specific aquatic habitat. As piscivory and adaptations for aquatic proclivities seem to be tied in vertebrates, disentangling these factors in dipsadine snakes will take further work and finer-scale ecological data.
The semi-fossorial group, with the exception of Heterodon and Xenodon dorbignyi, clustered fairly neatly away from other groups and towards a narrower skull in the non-trophic module analysis (Fig. 4A). This is consistent with Savitzky’s  finding that burrowers tend to have narrower, more rigid skulls, hypothesized to allow fossorial snakes to more effectively penetrate dense soils (see also  for detailed discussion on fossorial snake skull anatomy). Most semi-fossorial snakes also exhibit strong differentiation from other habitat use groups in their quadrate morphology – a finding corroborated by Palci et al.  in their study of squamate quadrate morphology, which again found that semi-fossorial snakes and lizards had significantly different quadrate morphology from terrestrial and aquatic species. With the exception of Heterodon spp., we found that semi-fossorial snakes tended to have short, squat quadrates (the opposite of aquatic snakes), which may be tied to functional constraints associated with evolving elongate, but more compact skulls appropriate for soil penetration. Mediolaterally narrower skull shapes for burrowing versus non-burrowing species have also been recorded in fishes  and skinks , and dorsoventrally flattened skulls in salamanders . Individuals with narrower skulls were shown to be faster burrowers in the caecilian Schistometopum thomense  and in amphisbaenians , although wider skull shapes are associated with higher absolute burrowing forces . A more angular, wedge-shaped rostrum, exhibited by the semi-fossorial Atractus, Heterodon, and Xenodon (E, G, and H in Fig. 2, respectively), has been found in other fossorial snakes , fossorial lizards , burrowing fish , and forward-burrowing frogs . In general, fossoriality seems to be a strong constraint on vertebrate skulls [96,97,98,99,100] often driving rampant homoplasy [45, 98, 101].
Our findings with respect to the diet variable corroborated some of the classic exemplars of trophic specialization in dipsadines. Mollusk and in particular snail specialists possess highly derived and unique traits that aid them in extracting and consuming their prey [52, 74, 76]. Our findings are consistent with this assessment, as in PC analyses of the ectopterygoid, maxilla, and palatine – all trophic components of the skull – mollusk specialists cluster in morphospace to some extent (Figs. 6, 7, S2, S4). Morphology associated with annelid prey specialization in snakes has not been discussed extensively to our knowledge. However, we suspect that the tendency we see for annelid specialists to deviate in morphospace in analyses of mandible and pterygoid shape (Figs. 7, S3, S4) relate to trophic constraints imposed by slippery, difficult to handle prey.
Our results suggest several instances of convergent evolution, evidenced both quantitatively in morphospace (Fig. 4) and upon inspection of the skulls themselves (Fig. 2). One of the most conspicuous examples of this phenomenon involves the semi-fossorial outliers: Heterodon and Xenodon dorbignyi. These two taxa represent phylogenetically and geographically disparate lineages that visibly converge on the same area of morphospace, both clustering, however, at the extreme opposite end of morphospace when compared to all other semi-fossorial species (Fig. 4A). Although both species are burrowers [102, 103], they probably spend most of their time above ground  and may use their expanded and prominent rostrum to unearth toad prey during the daytime . This is a very different life history than most truly fossorial or semi-fossorial snakes and may account for their morphological distinctiveness relative to other purportedly fossorial/semi-fossorial snakes. Strikingly similar skulls with expanded and upturned rostra have also evolved in another diurnal toad forager, Leioheterodon [105, 106], a phylogenetically distant lineage from Madagascar. Having independently evolved in several phylogenetically and geographically distinct lineages, it is probable that this highly specialized morphotype has strong ecological correlates, despite the fact that our data were unable to detect them.
Unfortunately, we lack even basic natural history data for many snakes [107, 108] and it is probable that the coarse ecological state codings used here have masked complexity in the relationship between habitat, diet, and skull morphology. As several authors have noted [108,109,110], better-quality and finer-scale ecological data are badly needed to understand snake diversification more generally . The paucity of ecological data for snakes results in part from their inherently cryptic nature and the difficulty in obtaining large sample sizes . Their subsequently poor reputation as ecological study models has resulted in comparatively few snake-focused ecological publications [107, 108], but see . A concerted effort to gather more and better ecological data for snakes is more than warranted, particularly when considering this taxon’s potential as a system for understanding large-scale biodiversity phenomena. Snakes comprise many overlooked but spectacular examples of adaptive radiation [61, 71], key innovation [48, 49, 60], evolutionary arms races [113,114,115], and Mullerian-Batesian mimicry rings [72, 116, 117]. As our understanding of snake natural history increases, we believe that they will come to represent a rich source of insight into macroevolutionary and ecological dynamics.
With dipsadine snakes, habitat use and diet are strong predictors of skull shape. This finding is consistent with the idea that dipsadines have radiated in response to ecological opportunity, an idea that is further supported by the tempo and mode of trophic niche evolution within the group . More evidence is required to definitively support the idea that the spectacular ecological and phenotypic diversity of the dipsadines is the result of adaptive radiation, however, the findings of this work serve to satisfy one of the requirements, namely, a correlation between ecology and morphology. We find that habitat ecology may play a more prominent overall role in constraining skull shape in dipsadines when compared to diet, although both significantly correlated with the shape of many regions of the skull. Fossorial and aquatic dipsadines separated from other habitat use groups most strongly; we argue that this is likely associated with the strong constraints imposed by the difficult mediums (soil and water, respectively) that these snakes tend to move through. Lastly, we identify an instance where our coarse ecological data were unable to illuminate the circumstances of a probable instance of convergence. As one of the largest continental vertebrate radiations, there is an acute need for comprehensive ecological, phenotypic, and phylogenetic information from dipsadines such that we might better understand both the evolutionary causes and macroecological consequences of their diversification.
We generated cranial images for 160 species of dipsadines using x-ray micro-computed tomography (hereafter, CT) scanning technology. We digitally constructed 3-dimensional surface models for each skull and quantified their shape using geometric morphometrics. We compiled ecological data from the primary literature and through a recently published predator-prey interaction database for snakes . We used principal components analysis to quantify and visualize the variation in skull shape across the dataset, after removing the effects of rotation and minimizing size differences in the shape data. We then used phylogenetic generalized least-squares (PGLS) and phylogenetic partial least-squares (PPLS) to explore the contribution of ecological variables to the observed variation in skull shape in a phylogenetic framework, while accounting for potential allometric effects. Finally, we applied a set of permutation and simulation based analyses to explore the robustness of the relationships between ecology and cranial morphology.
We collected 3D morphological images for 160 species of dipsadine snakes using micro-CT. We selected one representative specimen, free of skeletal damage, for each species (specimens examined in Additional file 1: Appendix S1). Specimens utilized were obtained primarily from the publicly available holdings of the University of Michigan Museum of Zoology – a full list of specimens used, with corresponding museum institution codes and catalog numbers, is available in the supplementary material (Additional file 1: Appendix S1). We selected only adult specimens, assessed by comparing material for the species, as juvenile snake skull morphology is significantly different from adult skull morphology in many cases [118,119,120,121] and specifically for dipsadines [122, 123], which could mislead an analysis of skull morphology across species. Many juvenile snakes exhibit classically paedomorphic features such as proportionally enlarged eyes and braincase and a larger head relative to their body size [118,119,120,121,122,123], and these traits often easily distinguish juveniles from adults when sufficient comparative material is available. We did not include species where adequate comparative material was not available and all specimens appeared to exhibit paedomorphic traits.
We obtained micro-CT scans primarily using a Nikon XT H225ST μCT machine, though several were imaged on a Scanco Medical µCT 100. Our scanning parameters varied depending on the specimen and facility, but ranged from 10-30 μm voxel size, 85-100 kV, 80-200 mA, 500-1600 projections, and 1-2 frame averaging. We reconstructed tomograms from raw radiograph projections on the Nikon XT H225ST machines using proprietary Nikon CT pro 3D software. We used the program Avizo to create 3-dimensional surface models on which landmarks could be placed.
We chose landmark locations (n=73) (Fig. 8, Additional file 2: Appendix S2) based on easily identifiable and homologous locations on skull elements that were present and unambiguous across all taxa. We attempted to capture the maximal amount of observed morphological variation across species with our landmarking scheme. We placed landmarks on both sides of the cranium, as there is evidence that landmarking only single halves of bilaterally symmetrical structures can influence the accuracy of geometric morphometrics analyses . Disjoint trophic structures that we analyzed separately were landmarked right side only. We landmarked surface models using Stratovan Checkpoint.
We obtained habitat use data through general literature surveys as well as from personal experience in observing species in the field (Additional file 3: Appendix S3). We qualitatively assessed each of the 160 dipsadine species as being one of the following: semi-fossorial (partially subterranean dwelling), cryptozoic (highly cryptic, found below leaf litter, logs etc.), terrestrial, semi-arboreal, arboreal, semi-aquatic, and aquatic.
We compiled dietary data in two separate ways. For a smaller species subset (n=68), we assessed the data quantitatively in the form of proportion of each prey category consumed. For the entire dataset (n=160), we assessed diet categorically in the form of primary prey category, utilizing a combined quantitative and qualitative approach. We extracted quantitative dietary data from SquamataBase, a database consisting of predator-prey interactions for over 1,200 species of snakes . Each of the records in this database represents a single observation of a prey item consumed by a snake. These records were either obtained from published literature or through dissections of museum specimens. These quantitative dietary data were available for 129 of the species in the dataset. Of these 129 species with quantitative dietary data, 68 species were represented by greater than 10 dietary records. We compiled the dietary data for this 68 species subset in the form of proportion prey consumed by each species in each of 10 categories: amphibians, reptiles, birds, mammals, fish, reptile eggs, bird eggs, annelids, and mollusks. Our diet-coding scheme closely follows prey taxonomy. For the entire 160 species dataset, we then assessed diet categorically. For species with quantitative diet data available, we assigned primary diet category based on which prey type accounted for over 50% of a given species diet. For species in which no one category made up more than 50% of the total prey composition, we assigned the category of “generalist”. For the remaining 31 species for which no quantitative diet data was available, we assessed primary diet category qualitatively based on the available literature. All ecological coding designations and corresponding references are available in Additional file 3: Appendix S3.
Morphology across species will to some extent be correlated as a consequence of shared evolutionary history; species thus do not represent independent data points that can be directly compared [125,126,127,128]. We accounted for phylogenetic structure in the data through the use of comparative methods incorporating phylogenetic relationships, extended for use with geometric morphometric data [129,130,131,132]. The phylogeny used in all analyses was a recent time-calibrated molecular phylogeny for 1263 species of advanced caenophidian snakes, produced utilizing a supermatrix approach in a maximum likelihood framework . From the full 160 species dataset, 148 species were present in the Zaher tree; consequently, we restricted statistical analyses to this 148 species subset but visualized morphospace for the full 160 taxon dataset.
We assessed variation in skull shape across the dataset after removing the effects of rotation and minimizing size differences in the landmark data. This was done by visualizing species positions in the morphospace defined by the first two principal components of the Procrustes-transformed landmark data. We then tested whether habitat use and diet composition were significant predictors of skull shape. We included size as an additional explanatory variable in our analyses, to account for any possible effects of allometry.
For analysis, we subset landmark coordinates for a total of 8 skull modules that we then analyzed in isolation: non-trophic elements (braincase, postorbitals, prefrontals, nasals and premaxilla), maxilla, ectopterygoid, supratemporal, quadrate, mandible, pterygoid, and palatine (Fig. 3). Because there are many kinetic elements in the snake skull that are capable of rotation and translation, a global analysis of overall skull shape would be confounded by the arbitrary position that these elements were preserved in. Various methods are available for removing the effects of arbitrary rotation and translation in articulated 3D structures [133,134,135]. However, the large number of extremely loose articulations in the snake skull, combined with their mobility on multiple axes, makes implementing some of these methods (i.e. ) inhibitive and ineffective for the 3-dimensional snake skull . Although other methods (e.g. [134, 135]) can accommodate complex, linking chains like those that occur in the snake cranium, we nonetheless chose to analyze the most mobile elements individually. This both removed the effect of arbitrary element rotation as a result of kinesis and allowed us to more easily detect fine-scale variation in each element independently.
We conducted analyses in R 4.2.2, utilizing the packages “geomorph v.4.0.5” [137, 138], “RRPP v1.3.1” [139, 140], “ape v.5.7” , “nlme v3.1-162” [142, 143], “phytools v1.5-1” , and “car v3.1.1” . We also created certain visualizations using functions from “diversitree v.0.9-16” . We provide all relevant data (landmarks, ecological data etc.) along with an annotated R script containing code used for all procedures described in this paper in the Dryad Digital Repository, under the https://doi.org/10.5061/dryad.6q573n63n.
We aligned landmarks and minimized size differences for each skull module across all 160 species using Generalized Procrustes Analysis [147, 148]. We conducted principal components analyses for all aligned modules separately, plotting PC scores for principal components 1 and 2 and differentially coloring points based on habitat use and primary diet to visualize variance in skull shape with respect to these ecological variables. We obtained the mean shape for the non-trophic skull module across all species and warped it using the thin-plate spline approach  to represent the shape extremes along the plotted PC axes. This essentially allows for the interpretation of shape variation tracked by each principal component axis. We also conducted principal components analyses for the subset of species with proportional diet data (n=68) across all aligned skull modules. We then plotted PC scores for PC 1 and 2 as pie chart points in morphospace representing each species’ diet composition (figures available in Additional file 4: Appendix S4).
Ecological predictors of cranial evolution
We tested the effects of habitat use and primary diet on skull shape using phylogenetic generalized least squares (PGLS) with permutation. As size (allometry) has been shown to be a significant contributing factor to shape in many groups [58, 78,79,80,81], we designated size as an additional explanatory variable. We used log of centroid size (the square root of the sum of squared distances of the landmarks to their centroid, obtained from general Procrustes analyses), as a proxy for skull size. We implemented this PGLS analysis using the function “procD.pgls” in geomorph [129, 131, 132]. We used a single model, with three independent variables (habitat use, primary diet, and size) and one dependent variable (skull shape). Modeling independent variables together allowed us to compare the predictive strength of these variables on the shape of each skull module. Because R2 values are difficult to interpret for multivariate data [Adams, pers. comm.], we computed and compared the effect size as a metric for comparing the strength of signals. As procD.pgls can accommodate 3-dimensional, aligned landmark coordinates, we used these to represent skull shape. We calculated covariance across species as a result of phylogenetic relationships under a Brownian motion model of evolution.
We applied a second PGLS implementation using “gls” in the nlme package, “Anova” in the car package, and simplified shape data, in order to provide an additional layer of robustness to the results. Whereas with procD.pgls we utilized the unaltered high-dimensional shape dataset, here we used PC scores for principal components 1 and 2 of skull shape as the shape variables in two separate analyses. In order to account for phylogenetic structure, we specified an evolutionary covariance matrix under Brownian motion, using the function “corBrownian” from the ape package. We utilized the same model structure and variables as with procD.pgls.
We assessed the correlation between diet in the form of proportional data and skull shape using two-block partial least squares in a phylogenetic context (PPLS). This method assesses the degree to which two sets of variables covary in a phylogenetic context, assessing significance with a permutation test [150, 151]. Thus, it is not a test that assesses the effect of independent variables on a dependent variable as in PGLS, but rather the degree to which the two sets of variables are correlated, with no assumption of directionality . We implemented PPLS using the function “phylo.integration” in geomorph [130, 132].
Availability of data and materials
We have included a zipped archive of all data files and scripts in our submission that we would like to make available during the review process. All data and associated scripts for analysis is also archived in the Dryad Data Repository, under the following https://doi.org/10.5061/dryad.6q573n63n.
Phylogenetic generalized least-squares
Phylogenetic partial least-squares
Schluter D. The ecology of adaptive radiation. Oxford (UK): Oxford University Press; 2000. p. 296.
Benkman CW. Divergent Selection Drives the Adaptive Radiation of Crossbills. Evolution. 2003;57(5):1176–81.
Stroud JT, Losos JB. Ecological Opportunity and Adaptive Radiation. Annu Rev Ecol Evol Syst. 2016;47:507–32.
Grossnickle DM, Newham E. Therian mammals experience an ecomorphological radiation during the Late Cretaceous and selective extinction at the K-Pg boundary. Proc R Soc B Biol Sci. 1832;2016(283):20160256.
Andjelković M, Tomović L, Ivanović A. Morphological integration of the kinetic skull in Natrix snakes. J Zool. 2017;303:188–98.
Felice RN, Watanabe A, Cuff AR, Noirault E, Polk D, Witmer LM, Norell MA, O’Connor PM, Goswami A. Evolutionary Integration and Modularity in the Archosaur Cranium. Integr Comp Biol. 2019;59:371–82.
Marroig G, Cheverud JM. A Comparison Of Phenotypic Variation And Covariation Patterns And The Role Of Phylogeny, Ecology, And Ontogeny During Cranial Evolution Of New World Monkeys. Evolution. 2001;55:2576–600.
Bright JA, Marugán-lobón J, Cobb SN, Rayfield EJ. The shapes of bird beaks are highly controlled by nondietary factors. P Natl Acad Sci U S A. 2016;113:5352–7.
Arbour JH, Curtis AA, Santana SE. Signatures of echolocation and dietary ecology in the adaptive evolution of skull shape in bats. Nat Commun. 2019;10:2036.
Cadle JE, Greene HW. Phylogenetic Patterns, Biogeography, and the Ecological Structure of Neotropical Snake Assemblages. In: Ricklefs RE, Schluter D, editors. Species diversity in ecological communities: Historical and geographical perspectives. Chicago (IL): University of Chicago Press; 1993. p. 281–93.
Uetz, P, Freed P, Aguilar R, Reyes F, Hošek J, editors. The Reptile Database. 2022. http://www.reptile-database.org. Accessed 30 Dec 2022.
Zaher H, Murphy RW, Arredondo JC, Graboski R, Machado-Filho PR, Mahlow K, et al. Large-scale molecular phylogeny, morphology, divergence-time estimation, and the fossil record of advanced caenophidian snakes (Squamata: Serpentes). PLoS One. 2019;14(5):e0216148.
Da Silva NJ, Sites JW. Patterns of diversity of neotropical squamate reptile species with emphasis on the brazilian amazon and the conservation potential of indigenous reserves. Conserv Biol. 1995;9:873–901.
Guedes TB, Nogueira C, Marques OA. Diversity, natural history, and geographic distribution of snakes in the Caatinga. Northeastern Brazil Zootaxa. 2014;3863(1):1–93.
Rabosky DL, von May R, Grundler MC, Davis Rabosky AR. The Western Amazonian richness gradient for squamate reptiles: Are There Really Fewer Snakes and Lizards in Southwestern Amazonian Lowlands? Diversity. 2019;11(10):199.
Kofron CP. Systematics of Neotropical gastropod-eating snakes: the dimidiata group of the genus Sibon, with comments on hte nebulata group. Amphib Reptil. 1990;11:207–23.
Angarita-Sierra T. Ninia atrata. Catálogo anfibios y Reptil Colomb. 2017;3(2):30–7.
Stafford PJ, Henderson RW. Ecological Traits of the Colubrid Snake Conophis Lineatus Concolor (Guarda Camino) in the Yucatán Peninsula. South Am J Herpetol. 2006;1(3):210–7.
Balestrin RL, Di-Bernardo M, Moreno AG. Feeding ecology of the neotropical worm snake Atractus reticulatus in southern Brazil. Herpetol J. 2007;17:62–4.
Marques OA, Oliveira JL. Predacao de ovos por Rachidelus brazili (Serpentes, Colubridae, Pseudoboini). Presented at: XXIII Congresso Brasileiro de Zoologia. Cuiabá: Universidade Federal de Mato Grosso; 2004.
Gaiarsa MP, de Alencar LR, Martins M. Natural History Of Pseudoboine Snakes. Pap Avulsos Zool. 2013;53:261–83.
Greene HW, Jaksic FM. The feeding behavior and natural history of two Chilean snakes, Philodryas chamissonís and Tachymenis chilensis (Colubridae). Rev Chil Hist Nat. 1992;65:485–93.
Rodriguez-Robles JA, Mulcahy DG, Greene HW. Feeding Ecology of the Desert Nightsnake, Hypsiglena torquata (Colubridae). Copeia. 1999;1999(1):93–100.
Obst FJ, Richter RK, Jacobs U. The completely illustrated atlas of reptiles and amphibians for the terrarium. Neptune (NJ): TFH Publications; 1988. p. 830.
Franca FG, Mesquita DO, Nogueira CC, Araújo AF. Phylogeny and Ecology Determine Morphological Structure in a Snake Assemblage in the Central Brazilian Cerrado. Copeia. 2008;2008(1):23–38.
Schwartz A, Henderson RW. Amphibians and Reptiles of the West Indies: Descriptions, Distributions, and Natural History. Gainesville (FL): University of Florida Press; 1991. p. 720.
Silveira AL, Ribeiro LS, Dornas TT, Fernandes TN. New records of dispas albifrons (Serpentes, Dipsadidae) in the atlantic forest of Minas Gerais, Brazil, with morphological data. Herpetol Notes. 2018;11:809–15.
Mount RH. The reptiles and amphibians of Alabama. Auburn (AL): Auburn University; 1975. p. 347.
Marques OA, Eterovic A, Sazima I. Snakes of the Brazilian Atlantic Forest: An Illustrated Field Guide for the Serra do Mar Range. Ribeirão Preto (BR): Holos Editora; 2004. p. 205.
Sampaio IL, Santos CP, França RC, Pedrosa IM, Solé M, França FGR. Ecological diversity of a snake assemblage from the atlantic forest at the south coast of paraíba, northeast Brazil. Zookeys. 2018;787:107–25.
Herrel A, Schaerlaeken V, Meyers JJ, Metzger KA, Ross CF. The evolution of cranial design and performance in squamates: Consequences of skull-bone reduction on feeding behavior. Integr Comp Biol. 2007;47:107–17.
Pierce SE, Angielczyk KD, Rayfield EJ. Patterns of Morphospace Occupation and Mechanical Performance in Extant Crocodilian Skulls : A Combined Geometric Morphometric and Finite Element Modeling Approach. J Morphol. 2008;269:840–64.
Herrel A, Meyers JJ, Vanhooydonck B. Relations between microhabitat use and limb shape in phrynosomatid lizards. Biol J Linn Soc. 2002;77:149–63.
Kohlsdorf T, Grizante MB, Navas CA, Herrel A. Head shape evolution in Tropidurinae lizards : does locomotion constrain diet ? J Evol Biol. 2008;21:781–90.
Preston BT, Stevenson IR, Pemberton JM, Coltman DW, Wilson K. Overt and covert competition in a promiscuous mammal : the importance of weaponry and testes size to male reproductive success. Proc R Soc B. 2003;270:633–40.
Lappin AK, Husak JF. Weapon Performance, Not Size, Determines Mating Success and Potential Reproductive Output in the Collared Lizard ( Crotaphytus collaris ). Am Nat. 2005;166:426–36.
Cooper CE, Vitt LJ. Female mate choice of large male broad-headed skinks. Anim Behav. 1993;45:683–93.
Greene HW. Snakes: The Evolution of Mystery in Nature. London (UK): University of California Press, Ltd; 1997. p. 366.
Dumont ER, Herrel A, Medellín RA, Vargas-Contreras JA, Santana SE. Built to bite: Cranial design and function in the wrinkle-faced bat. J Zool. 2009;279:329–37.
Santana SE, Dumont ER, Davis JL. Mechanics of bite force production and its relationship to diet in bats. Funct Ecol. 2010;24:776–84.
Kulemeyer C, Asbahr K, Gunz P, Frahnert S, Bairlein F. Functional morphology and integration of corvid skulls – a 3D geometric morphometric approach. Front Zool. 2009;6(1):1–4.
Buser TJ, Burns MD, López JA. Littorally adaptive ? Testing the link between habitat, morphology, and reproduction in the intertidal sculpin subfamily Oligocottinae (Pisces : Cottoidea). PeerJ. 2017;5:e3634.
Buser TJ, Sidlauskas BL, Summers AP. 2D or Not 2D? Testing the Utility of 2D Vs. 3D Landmark Data in Geometric Morphometrics of the Sculpin Subfamily Oligocottinae (Pisces; Cottoidea). Anat Rec. 2017;301:806–18.
Stayton CT. Morphological Evolution of the Lizard Skull: A Geometric Morphometrics Survey. J Morphol. 2005;263:47–59.
Barros FC, Herrel A, Kohlsdorf T. Head shape evolution in Gymnophthalmidae: Does habitat use constrain the evolution of cranial design in fossorial lizards? J Evolution Biol. 2011;24:2423–33.
Gans C. The Feeding Mechanism of Snakes and Its Possible Evolution. Am Zool. 1961;1:217–27.
Kardong KV. “Protovipers” and the Evolution of Snake Fangs. Evolution. 1979;33(1):433–43.
Savitzky AH. The Role of Venom Delivery Strategies in Snake Evolution. Evolution. 1980;34:1194–204.
Cundall D, Greene HW. Feeding in snakes. In: Schwenk K, editor. Feeding: Form Function, and Evolution in Tetrapod Vertebrates. San Diego (CA): Academic Press; 2000. p. 293–333.
Lillywhite HB. How Snakes Work: Structure, Function and Behavior of the World’s Snakes. Oxford (UK) and New York (NY): Oxford University Press; 2014. p. 256.
King RB. Predicted and observed maximum prey size - Snake size allometry. Funct Ecol. 2002;16:766–72.
Savitzky AH. Coadapted Character Complexes Among Snakes. Am Zool. 1983;23:397–409.
Fabre AC, Bickford D, Segall M, Herrel A. The impact of diet, habitat use, and behaviour on head shape evolution in homalopsid snakes. Biol J Linn Soc. 2016;118:634–47.
Vincent SE, Brandley MC, Herrel A, Alfaro ME. Convergence in trophic morphology and feeding performance among piscivorous natricine snakes. J Evol Biol. 2009;22:1203–11.
Hampton PM. Comparison of Cranial Form and Function in Association with Diet in Natricine Snakes. J Morphol. 2011;272:1435–43.
Andjelković M, Tomović L, Ivanović A. Variation in skull size and shape of two snake species (Natrix natrix and Natrix tessellata). Zoomorphology. 2016;135:243–53.
Klaczko J, Sherratt E, Setz EZ. Are diet preferences associated to skulls shape diversification in xenodontine snakes? PLoS One. 2016;11(2):e0148375.
Sherratt E, Sanders KL, Watson A, Hutchinson MN, Lee MS, Palci A. Heterochronic Shifts Mediate Ecomorphological Convergence in Skull Shape of Microcephalic Sea Snakes. Integr Comp Biol. 2019;59(3):616–24.
Watanabe A, Fabre AC, Felice RN, Maisano JA, Müller J, Herrel A, et al. Ecomorphological diversification in squamates from conserved pattern of cranial integration. Proc Natl Acad Sci U S A. 2019;116(29):14688–97.
Westeen EP, Durso AM, Grundler MC, Rabosky DL, Davis Rabosky AR. What makes a fang? Phylogenetic and ecological controls on tooth evolution in rear-fanged snakes. BMC Evol Biol. 2020;20(1):1–5.
Grundler MC, Rabosky DL. Trophic divergence despite morphological convergence in a continental radiation of snakes. Proc R Soc B Biol Sci. 2014;281(1787):20140413.
Esquerre D, Keogh JS. Parallel selective pressures drive convergent diversification of phenotypes in pythons and boas. Ecol Lett. 2016;19:800–9.
Segall M, Cornette R, Fabre AC, Godoy-Diana R, Herrel A. Does aquatic foraging impact head shape evolution in snakes? Proc R Soc B. 2016;283:20161645.
Silva FM, da Costa Prudente AL, Machado FA, Santos MM, Zaher H, Hingst-Zaher E. Aquatic adaptations in a Neotropical coral snake- A study of morphological convergence. J Zool Syst Evol Res. 2018;56:382–94.
Da Silva FO, Fabre AC, Savriama Y, Ollonen J, Mahlow K, Herrel A, et al. The ecological origins of snakes as revealed by skull evolution. Nat Commun. 2018;9:376.
Palci A, Caldwell MW, Hutchinson MN, Konishi T, Lee MSY. The morphological diversity of the quadrate bone in squamate reptiles as revealed by high-resolution computed tomography and geometric morphometrics. J Anat. 2020;236:210–27.
Harmon LJ, Weir JT, Brock CD, Glor RE, Challenger W. GEIGER: investigating evolutionary radiations. Bioinformatics. 2008;24:129–31.
Pennell M, Eastman J, Slater G, Brown J, Uyeda J, Fitzjohn R, et al. Geiger v2.0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. Bioinformatics. 2014;30:2216–8.
Figueroa A, Mckelvy AD, Grismer LL, Bell CD. A Species-Level Phylogeny of Extant Snakes with Description of a New Colubrid Subfamily and Genus. PLoS One. 2016;11(9):e0161070.
Yoder JB, Clancey E, Des Roches S, Eastman JM, Gentry L, Godsoe W, et al. Ecological opportunity and the origin of adaptive radiations. J Evol Biol. 2010;23:1581–96.
Grundler MC, Rabosky DL. Rapid increase in snake dietary diversity and complexity following the end-Cretaceous mass extinction. PLoS Biol. 2021;19(10):e3001414.
Davis Rabosky AR, Cox CL, Rabosky DL, Title PO, Holmes IA, Feldman A, et al. Coral snakes predict the evolution of mimicry across New World snakes. Nat Commun. 2016;7(1):11484.
Greer AE. The Biology and Evolution of Australian Snakes. Chipping Norton (AU): Surrey Beatty & Sons Pty Ltd; 1997. p. 358.
dos Santos MM, da Silva FM, Hingst-Zaher E, Machado FA, Zaher HE, da Costa Prudente AL. Cranial adaptations for feeding on snails in species of Sibynomorphus ( Dipsadidae : Dipsadinae ). Zoology. 2017;120:24–30.
de Oliveira L, da Costa Prudente AL, Zaher H. Unusual Labial Glands in Snakes of the Genus Geophis Wagler, 1830 (Serpentes- Dipsadinae). J Morphol. 2014;275:87–99.
Zaher H, de Oliveira L, Grazziotin FG, Campagner M, Jared C, Antoniazzi MM, et al. Consuming viscous prey: A novel protein-secreting delivery system in neotropical snail-eating snakes. BMC Evol Biol. 2014;14(58):1–28.
Pigliucci M. Is evolvability evolvable? Nat Rev Genet. 2008;9(1):75–82.
Chamoli U, Wroe S. Allometry in the distribution of material properties and geometry of the felid skull: Why larger species may need to change and how they may achieve it. J Theor Biol [Internet]. 2011;283:217–26.
Paluh DJ, Bauer AM. Phylogenetic history, allometry and disparate functional pressures influence the morphological diversification of the gekkotan quadrate, a keystone cranial element. Biol J Linn Soc. 2018;125:693–708.
Gray JA, Sherratt E, Hutchinson MN, Jones MEH. Changes in ontogenetic patterns facilitate diversification in skull shape of Australian agamid lizards. BMC Evol Biol. 2019;19:7.
Chatterji RM, Hipsley CA, Sherratt E, Hutchinson MN, Jones MEH. Ontogenetic allometry underlies trophic diversity in sea turtles (Chelonioidea). Evol Ecol. 2022;36:511–40 (Springer International Publishing).
Deepak V, Gower DJ, Cooper N. Diet and habit explain head-shape convergences in natricine snakes. J Evol Biol. 2022;00:1–13.
Herrel A, Vincent SE, Alfaro ME, Van Wassenbergh S, Vanhooydonck B, Irschick DJ. Morphological convergence as a consequence of extreme functional demands: Examples from the feeding system of natricine snakes. J Evol Biol. 2008;21:1438–48.
McCurry MR, Evans AR, Fitzgerald EMG, Adams JW, Clausen PD, McHenry CR. The remarkable convergence of skull shape in crocodilians and toothed whales. Proc R Soc B. 2017;284:20162348.
Strong CRC, Palci A, Caldwell MW. Insights into skull evolution in fossorial snakes, as revealed by the cranial morphology of Atractaspis irregularis (Serpentes- Colubroidea). J Anat. 2021;238:146–72.
Evans KM, Larouche O, West JL, Gartner SM, Westneat MW. Burrowing constrains patterns of skull shape evolution in wrasses. Evol Dev. 2023;25:73–84.
Stepanova N, Bauer AM. Phylogenetic history influences convergence for a specialized ecology: comparative skull morphology of African burrowing skinks (Squamata; Scincidae). BMC Ecol Evol. 2021;21:86.
Ehmcke J, Clemen G. The structure and development of the skull of Costa Rican plethodontid salamanders (Amphibia: Urodela). Ann Anat. 2000;182:537–47.
Teodecki EE, Brodie ED Jr, Formanowicz DR Jr, Nussbaum RA. Head dimorphism and burrowing speed in the African caecilian Schistometopum thomense (Amphibia: Gymnophiona). Herpetologica. 1998;54:154–60.
dos Hohl L SL, de Loguercio MF C, Sicuro FL, de Barros-Filho JD, Rocha-Barbosa O. Body and skull morphometric variations between two shovel-headed species of Amphisbaenia (Reptilia: Squamata) with morphofunctional inferences on burrowing. PeerJ. 2017;5:e3581.
Lowie A, De Kegel B, Wilkinson M, Measey J, O’Reilly JC, Kley NJ, et al. Under pressure: The relationship between cranial shape and burrowing force in caecilians (gymnophiona). J Exp Biol. 2021;224:jeb242964.
Herrel A, Lowie A, Miralles A, Gaucher P, Kley NJ, Measey J, et al. Burrowing in blindsnakes: A preliminary analysis of burrowing forces and consequences for the evolution of morphology. Anat Rec. 2021;304:2292–302.
dos Hohl L SL, de Loguercio MF C, Buendia RA, Almeida-Santos M, Viana LA, Barros-Filho JD, et al. Fossorial gait patterns and performance of a shovel-headed amphisbaenian. J Zool. 2014;294:234–40.
Berra TM, Allen GR. Burrowing, Emergence, Behavior, and Functional Morphology of the Australian Salamanderfish Lepidogalaxias salamandroides. Fisheries. 1989;14(5):2–10.
Keeffe R, Blackburn DC. Comparative morphology of the humerus in forward-burrowing frogs. Biol J Linn Soc. 2020;131:291–303.
Navas CA, Antoniazzi MM, Carvalho JE, Chaui-Berlink JG, James RS, Jared C, et al. Morphological and physiological specialization for digging in amphisbaenians, an ancient lineage of fossorial vertebrates. J Exp Biol. 2004;207:2433–41.
Vanhooydonck B, Boistel R, Fernandez V, Herrel A. Push and bite: Trade-Offs between burrowing and biting in a burrowing skink (Acontias percivali). Biol J Linn Soc. 2011;102:91–9.
Sherratt E, Gower DJ, Klingenberg CP, Wilkinson M. Evolution of Cranial Shape in Caecilians (Amphibia: Gymnophiona). Evol Biol. 2014;41:528–45.
Kazi S, Hipsley CA. Conserved evolution of skull shape in Caribbean headfirst burrowing worm lizards (Squamata: Amphisbaenia). Biol J Linn Soc. 2018;125:14–29.
Gomes Rodrigues H, Šumbera R, Hautier L. Life in Burrows Channelled the Morphological Evolution of the Skull in Rodents: the Case of African Mole-Rats (Bathyergidae, Rodentia). J Mamm Evol. 2016;23:175–89.
Marcy AE, Hadly EA, Sherratt E, Garland K, Weisbecker V. Getting a head in hard soils: Convergent skull evolution and divergent allometric patterns explain shape variation in a highly diverse genus of pocket gophers (Thomomys). BMC Evol Biol. 2016;16:207.
Edgren RA. The Natural History of the Hog-Nosed Snakes Genus Heterodon : A Review. Herpetologica. 1955;11:105–17.
Tozetti AM, de Oliveira RB, Pontes GM. Defensive repertoire of Xenodon dorbignyi (Serpentes, Dipsadidae). Biota Neotrop. 2010;9:157–63.
Plummer MV, Mills NE. Society for the Study of Amphibians and Reptiles Spatial Ecology and Survivorship of Resident and Translocated Hognose Snakes ( Heterodon platirhinos ). J Herpetol. 2000;34:565–75.
Glaw F, Vences M. A Field Guide to the Amphibians and Reptiles of Madagascar. 3rd ed. Köln (DE): Vences & Glaw Verlag; 2007. p. 496.
Spain M, Fuller G, Allard S. Effects of Habitat Modifications on Behavioral Indicators of Welfare for Madagascar Giant Hognose Snakes (Leioheterodon madagascariensis). Anim Behav Cogn. 2020;7(1):70–81.
Mullin SJ, Siegel RA, editors. Snakes: Ecology and Conservation. Ithaca and London: Cornell University Press; 2009. p. 384.
Alencar LR, Gaiarsa MP, Martins M. The Evolution of Diet and Microhabitat use in Pseudoboine Snakes. South Am J Herpetol. 2013;8(1):60–6.
Greene HW. Organisms in nature as a central focus for biology. Trends Ecol Evol. 2005;20(1):23–7.
Vitt LJ, Pianka ER. Deep history impacts present-day ecology and biodiversity. Proc Natl Acad Sci U S A. 2005;102(22):7877–81.
Durso AM, Willson JD, Winne CT. Needles in haystacks: Estimating detection probability and occupancy of rare and cryptic snakes. Biol Conserv. 2011;144:1508–15.
Grundler MC. SquamataBase: a natural history database and R package for comparative biology of snake feeding habits. Biodivers Data J. 2020;8:e49943.
Brodie ED, Brodie ED. Tetrodotoxin resistance in garter snakes: an evolutionary response of predators to dangerous prey. Evolution. 1990;44(3):651–9.
Brodie ED, Feldman CR, Hanifin CT, Motychak JE, Mulcahy DG, Williams BL, et al. Parallel arms races between garter snakes and newts involving tetrodotoxin as the phenotypic interface of coevolution. J Chem Ecol. 2005;31(2):343–56.
Balchan NR. Resistance to rattlesnake venoms in an eastern Colorado rodent community [master’s thesis]. Greeley: University of Northern Colorado; 2021. p. 143.
Sanders KL, Malhotra A, Thorpe RS. Evidence for a Müllerian mimetic radiation in Asian pitvipers. Proc R Soc B Biol Sci. 2006;273:1135–41.
Davis Rabosky AR, Cox CL, Rabosky DL. Unlinked Mendelian inheritance of red and black pigmentation in snakes: Implications for Batesian mimicry. Evolution. 2016;70(4):944–53.
Vincent SE, Moon BR, Herrel A, Kley NJ. Are ontogenetic shifts in diet linked to shifts in feeding mechanics? Scaling of the feeding apparatus in the banded watersnake Nerodia fasciata. J Exp Biol. 2007;210:2057–69.
Palci A, Lee MSY, Hutchinson MN. Patterns of postnatal ontogeny of the skull and lower jaw of snakes as revealed by micro-CT scan data and three-dimensional geometric morphometrics. J Anat. 2016;229:723–54.
Esquerre D, Sherratt E, Keogh JS. Evolution of extreme ontogenetic allometric diversity and heterochrony in pythons, a clade of giant and dwarf snakes. Evolution. 2017;71(12):2829–44.
Patterson M, Wolfe AK, Fleming PA, Bateman PW, Martin ML, Sherratt E, et al. Ontogenetic shift in diet of a large elapid snake is facilitated by allometric change in skull morphology. Evol Ecol. 2022;36:489–509.
Murta-Fonseca RA, Fernandes DS. The skull of Hydrodynastes gigas (Duméril, Bibron & Duméril, 1854) (Serpentes: Dipsadidae) as a model of snake ontogenetic allometry inferred by geometric morphometrics. Zoomorphology. 2016;135:233–41.
Abegg AD, Passos P, Mario-da-Rosa C, dos Azevedo W S, Malta-Borges L, de Moura Bubadué J. Sexual dimorphism, ontogeny and static allometry of a semi-fossorial snake (genus Atractus). Zool Anz. 2020;287:95–104.
Cardini A. Lost in the Other Half: Improving Accuracy in Geometric Morphometric Analyses of One Side of Bilaterally Symmetric Structures. Syst Biol. 2016;65:1096–106.
Ridley M. The explanation of organic diversity. Oxford: Oxford University Press; 1983. p. 272.
Felsenstein J. Phylogenies and the comparative method. Am Nat. 1985;125:1–15.
Huey RB. Phylogeny, history and the comparative method. In: Feder ME, Bennett AF, Burggren WW, Huey RB, editors. New directions in ecological physiology. Cambridge: Cambridge University Press; 1987. p. 76–101.
Harvey PH, Pagel MD. The comparative method in evolutionary biology. Oxford: Oxford University Press; 1991. p. 248.
Adams DC. A method for assessing phylogenetic least squares models for shape and other high-dimensional multivariate data. Evolution. 2014;68:2675–88.
Adams DC, Felice R. Assessing phylogenetic morphological integration and trait co-variation in morphometric data using evolutionary covariance matrices. PLoS One. 2014;9(4):e94335.
Adams DC, Collyer ML. Permutation tests for phylogenetic comparative analyses of high-dimensional shape data: what you shuffle matters. Evolution. 2015;69:823–9.
Adams DC, Collyer ML. Multivariate comparative methods: evaluations, comparisons, and recommendations. Syst Biol. 2018;67:14–31.
Vidal-García M, Bandara L, Keogh JS. ShapeRotator: An R tool for standardized rigid rotations of articulated three-dimensional structures with application for geometric morphometrics. Ecol Evol. 2018;8:4669–75.
Collyer ML, Davis MA, Adams DC. Making heads or tails of combined landmark configurations in geometric morphometric data. Evol Biol. 2020;47:193–205.
Rhoda D, Polly PD, Raxworthy C, Segall M. Morphological integration and modularity in the hyperkinetic feeding system of aquatic-foraging snakes. Evolution. 2021;75:56–72.
Rhoda D, Segall M, Larouche O, Evans K, Angielczyk KD. Local Superimpositions Facilitate Morphometric Analysis of Complex Articulating Structures. Integr Comp Biol. 2021;61(5):1892–904.
Baken E, Collyer M, Kaliontzopoulou A, Adams D. geomorph v4.0 and gmShiny: enhanced analytics and a new graphical interface for a comprehensive morphometric experience. Methods Ecol Evol. 2021;12:2355–63.
Adams DC, Collyer ML, Kaliontzopoulou A, Baken E. Geomorph: Software for geometric morphometric analyses, R package version 4.0.4. 2022. Available from: https://cran.r-project.org/package=geomorph
Collyer ML, Adams DC. RRPP: An r package for fitting linear models to high-dimensional data using residual randomization. Methods Ecol Evol. 2018;9:1772–9.
Collyer ML, Adams DC. RRPP: Linear Model Evaluation with Randomized Residuals in a Permutation Procedure, R package version 1.1.2. 2021. Available from: https://cran.r-project.org/package=RRPP
Paradis E, Schliep K. Ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics. 2019;35:526–8.
Pinheiro J, Bates D, R Core Team. nlme: Linear and Nonlinear Mixed Effects Models, R package version 3.1–162. 2023. Available from: https://CRAN.R-project.org/package=nlme
Pinheiro JC, Bates DM. Mixed-Effects Models in S and S-PLUS. New York (NY): Springer; 2000. p. 528.
Revell L. phytools: An R package for phylogenetic comparative biology (and other things). Methods Ecol Evol. 2012;3:217–23.
Fox J, Weisberg S. An R Companion to Applied Regression. 3rd ed. Thousand Oaks (CA): Sage; 2019. p. 608.
FitzJohn RG. Diversitree: Comparative phylogenetic analyses of diversification in R. Methods Ecol Evol. 2012;3:1084–92.
Gower JC. Generalized Procrustes analysis. Psychometrika. 1975;40:33–51.
Rohlf FJ, Slice DE. Extensions of the Procrustes method for the optimal superimposition of landmarks. Syst Zool. 1990;39:40–59.
Bookstein FL. Principal Warps: Thin-Plate Splines and the Decomposition of Deformations. IEEE Trans Pattern Anal Mach Intell. 1989;11(6):567–85.
Sampson PD, Streissguth AP, Barr HM, Bookstein FL. Neurobehavioral effects of prenatal alcohol: part II partial least squares analysis. Neurotoxicol Teratol. 1989;11:477–91.
Streissguth AP, Bookstein FL, Sampson PD, Barr HM. The enduring effects of pre- natal alcohol exposure on child development: birth through seven years, a partial least squares solution. Ann Arbor (MI): University of Michigan Press; 1993. p. 336.
Rohlf JF, Corti M. Use of Two-Block Partial Least-Squares to Study Covariation in Shape. Syst Biol. 2000;49:740–53.
We thank Greg Schneider for facilitating the use of specimens from the University of Michigan Museum of Zoology, which were crucial in completing this study. Walter Schargel, Rudolph von May, and Dimitrios Pandelis assisted in diagnosing some of the statistical anomalies encountered. We especially thank Dean Adams for detailed comments and discussion of the statistical analyses that were important in developing our approach to statistical robustness. Ray Chatterji provided invaluable discussion of methods in allometry that improved our modeling approach. Alison Davis Rabosky and Cody Thompson provided comments on preliminary versions of this manuscript, which improved the research presented here.
Data collection and analysis was supported by a David and Lucile Packard Foundation Fellowship to Daniel Rabosky and the oVert TCN project (oVert TCN; NSF DBI-1701714). Fieldwork was also supported by a Schweitzer Honor Travel Grant and EEB Summer Research and Travel Award to Gregory Pandelis.
Ethics approval and consent to participate
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Specimens used in this study, with corresponding museum voucher numbers.
A list and description of all landmark positions utilized in this study.
Codings for habitat and major diet ecology for all species analyzed in this study, including literature references used to infer these data.
Variance explained by the first ten PC axes across skull modules, for the full 160 species dataset. Note that variance explained for most skull modules, even those with high numbers of landmarks, drops significantly after PC4. Fig. S1. A snapshot of the ecological diversity of the Dipsadine mega-radiation with taxon labels. Fig. S2. Principal components of shape (PC1, PC2) in three modules: non-trophic (A), maxilla (B), ectopterygoid (C), with respect to proportional diet composition, illustrated as pie-chart points. Fig. S3. Principal components of shape (PC1, PC2) in three modules: supratemporal (A), quadrate (B), mandible (C), with respect to proportional diet composition, illustrated as pie-chart points. Fig. S4. Principal components of shape (PC1, PC2) in two modules: pterygoid (top) and palatine (bottom) with respect to proportional diet composition, illustrated as pie-chart points. Fig. S5. P-values from gls (nlme) analyses of real datasets with true corresponding ecological data (left), those from analyses of real datasets paired with randomly permuted ecological data (middle), and those from analyses of BM-simulated shape data paired with randomly permuted ecological data (right). P-values are lower than expected for randomly permuted real datasets, an issue exacerbated with increased sample size. Retaining everything constant but simulating the shape data under BM produces expected p-values. Fig. S6. Allometry plots (regression of size against shape) in four modules: non-trophic (A), maxilla (B), ectopterygoid (C), and supratemporal (D) with respect to habitat use (1) and primary diet (2). Allometry appears to be the strongest in the non-trophic module, with looser allometric relationships visible in the maxilla and supratemporal. A rough pattern of isometry is visible in the ectopterygoid. Fig. S7 . Allometry plots (regression of size against shape) in four modules: quadrate (A), mandible (B), pterygoid (C), and palatine (D) with respect to habitat use (1) and primary diet (2). Allometry appears to be the strongest in the quadrate and mandible, with a looser pattern between size and shape in the pterygoid and palatine.
About this article
Cite this article
Pandelis, G.G., Grundler, M.C. & Rabosky, D.L. Ecological correlates of cranial evolution in the megaradiation of dipsadine snakes. BMC Ecol Evo 23, 48 (2023). https://doi.org/10.1186/s12862-023-02157-3