Exploring the functional meaning of head shape disparity in aquatic snakes

Abstract Phenotypic diversity, or disparity, can be explained by simple genetic drift or, if functional constraints are strong, by selection for ecologically relevant phenotypes. We here studied phenotypic disparity in head shape in aquatic snakes. We investigated whether conflicting selective pressures related to different functions have driven shape diversity and explore whether similar phenotypes may give rise to the same functional output (i.e., many‐to‐one mapping of form to function). We focused on the head shape of aquatically foraging snakes as they fulfill several fitness‐relevant functions and show a large amount of morphological variability. We used 3D surface scanning and 3D geometric morphometrics to compare the head shape of 62 species in a phylogenetic context. We first tested whether diet specialization and size are drivers of head shape diversification. Next, we tested for many‐to‐one mapping by comparing the hydrodynamic efficiency of head shape characteristic of the main axes of variation in the dataset. We 3D printed these shapes and measured the forces at play during a frontal strike. Our results show that diet and size explain only a small amount of shape variation. Shapes did not fully functionally converge as more specialized aquatic species evolved a more efficient head shape than others. The shape disparity observed could thus reflect a process of niche specialization.

origin of shape disparity in organisms that demonstrate parallel evolution, we need to investigate the interplay between ecological and functional constraints. Feeding under water is a particularly interesting case as strong functional constraints are imposed by the physical properties of water. The feeding apparatus of fully aquatic vertebrates, such as fish, either has morphologically or functionally converged (i.e., manyto-one-mapping) in response to the hydrodynamic constraints involved during prey capture (Collar & Wainwright, 2006;Cooper et al., 2010;Stuart et al., 2017;Thompson et al., 2017;Wainwright et al., 2005;Winemiller, Kelso-Winemiller, & Brenkert, 1995). In the present study, we investigate the interplay between different selective pressures that may have generated shape diversity (i.e., disparity) in a complex and integrated system in an ecologically diverse group: snakes.
While convergence was expected, the head shape of aquatic foragers has diverged from their fully terrestrial relatives, but instead of converging toward a unique shape, this group demonstrates an unexpectedly large head shape variability (Segall, Cornette, Fabre, Godoy-Diana, & Herrel, 2016), ranging from very slender (e.g., Thamnophis sp.) to very bulky heads (Laticauda sp., Aipysurus sp.). Aquatically foraging snakes are both species-rich and ecologically rich and have fast rates of evolution (Sanders, Lee, Mumpuni, Bertozzi, & Rasmussen, 2013;Watanabe et al., 2019). To understand the origin and drivers of the morphological diversity of the head of snakes, we explore two hypotheses: (1) the head shape of aquatically foraging snakes has diversified in response to functional constraints related to diet specialization, (2) this diversification has been facilitated by a many-to-one mapping of form to function allowing multiple head shapes to be equally efficient at reducing the hydrodynamic constraints related to a strike under water.
First, we focused on the impact of diet-related functional constraints (i.e., manipulation and swallowing) on the head shape of snakes. Morphological adaptation to diet-related constraints is widespread in vertebrates, from the spectacular adaptive radiation in the beak of Darwin's finches, the head of cichlid fish (Cooper et al., 2010), and the skull and mandible of mammals (Monteiro & Nogueira, 2009). Snakes are gape-limited predators that swallow prey whole (Gans, 1961), meaning that the size and shape of their head directly impacts the size and shape of prey they can eat. As snakes are vulnerable to both predator attack and injuries by their prey during prey manipulation and intraoral prey transport, they must reduce the time spent swallowing their prey. Previous studies have demonstrated a link between dietary preference and head shape in snakes (Camilleri & Shine, 1990;Fabre, Bickford, Segall, & Herrel, 2016;Forsman, 1991Forsman, , 1996Klaczko, Sherratt, & Setz, 2016;Queral-Regil & King, 1998;Sherratt, Rasmussen, & Sanders, 2018;Vincent, Moon, Herrel, & Kley, 2007). Although most of these studies used taxonomic groups to characterize snake diets (e.g., mammals, fish, anurans, crustaceans), this may be insufficient from a functional point of view (Vincent, Moon, Shine, & Herrel, 2006). Therefore, we here classified diet by characterizing the shape of the main prey eaten by each species: elongated or bulky. The ingestion of bulky prey, such as frogs, is more difficult for snakes (Vincent, Moon, et al., 2006), and the results from previous work on viperids and homalopsids suggest that bulkyprey eaters should benefit from wider and broader heads compared to elongated prey eaters (Brecko, Vervust, Herrel, & Van Damme, 2011;Fabre et al., 2016;Forsman, 1991;Vincent, Herrel, & Irschick, 2004).
In contrast, to reduce ingestion time, elongated prey eaters might benefit from elongated jaws which would reduce the number of jaw cycles required to swallow a long prey (Vincent, Moon, et al., 2006). As head size is expected to directly impact feeding efficiency in gape-limited predators like snakes (Esquerré & Keogh, 2016;Forsman, 1996;Grundler & Rabosky, 2014), we also quantified the evolutionary allometry in our dataset.
In the second part of this study, we explored the functional implications of the observed shape variation. All considered species successfully capture aquatic prey despite the hydrodynamic constraints they face Van Wassenbergh et al., 2010). As these constraints are related to head shape (Fish, 2004;Godoy-Diana & Thiria, 2018;Koehl, 1996;Polly et al., 2016), we expected the observed morphological disparity to have functionally converged (i.e., have the same hydrodynamic profile) which would indicate a many-to-one-mapping of form to function (Wainwright et al., 2005). We here defined the aquatic strike as our function of interest, and the performance indicators are the drag and added mass coefficient (i.e., the hydrodynamic profile) associated with the head shape of snakes. Drag is the force that resists the motion and is involved in all locomotor behavior, whereas added mass is involved only during acceleration. While drag has been extensively studied (Bale, Hao, Bhalla, & Patankar, 2014;Fish, 1993Fish, , 2000Godoy-Diana & Thiria, 2018;Stayton, 2011;Webb, 1988), added mass has been mostly ignored to date despite evidence of its major role in energy expenditure during locomotion (Vogel, 1994). For instance, 90% of the resistive force generated by the escape response of a crayfish is caused by its own mass and added mass, while drag represents the remaining 10% (Webb, 1979). Vogel (1994) suggested that propulsion-based organisms should be under a selective regime that favors a reduction in acceleration reaction by reducing mass and/or added mass.
Both drag and acceleration reaction are linked to the properties of the fluid, the kinematics of the motion, a scaling component, and a shape component. If snakes are under selection and yet display a large head shape disparity, then we can expect a many-to-one mapping of form to function, with several shapes resulting in a reduction in both drag and added mass. To test this hypothesis, we measured the shape-component of both hydrodynamic forces (i.e., drag and added mass coefficient) of shapes that are representative of the morphological disparity of our dataset.
We first compared the head shape of 62 species of snakes that capture elusive aquatic prey under water by scanning the surface of the head of more than 300 specimens from museum collections.
We then used high-density 3D geometric morphometric and phylogenetic comparative analyses to test the impact of diet and size on the head shape of snakes. We subsequently 3D printed five models of the head of snakes corresponding to the extremes of the main axes of variability (i.e., the two first principal components and the mean shape). We built an experiment that mimics a frontal strike, and we calculated the hydrodynamic efficiency of each shape to assess whether morphological disparity is associated with a functional convergence.

| Specimens
We scanned the head of 316 snakes belonging to 62 species of snakes that consume elusive aquatic prey (e.g., fish, amphibians, crus-taceans…) using a high-resolution surface scanner (Stereoscan3D Breuckmann white light fringe scanner with a camera resolution of 1.4 megapixels) at the morphometric platform of the Muséum National d'Histoire Naturelle, Paris (see Figure 1; Supplementary Material 1 for a list of specimens). Only specimens with a well-preserved head and closed mouth were scanned to allow shape comparisons. We chose the species to cover the diversity of aquatic snakes across the phylogeny (Pyron & Burbrink, 2014). The phylogenetic tree of Pyron and Burbrink (2014) was pruned in Mesquite 3.03 (Maddison & Maddison, 2015) (Figure 1). We described the diet of each species based on the available literature and attributed a main prey shape to each species depending on the length and shape of the maximal cross section of the prey. We defined two categories: Elongated prey are the items with a nearly circular cross section and a body length of more than twice the size of the longest dimension of the cross section (e.g., eels, gobiid fish, caecilians, tadpoles, snakes); bulky prey have either a noncircular cross section or a short, stout body (e.g., flattened fish, anurans) or represent a manipulation challenge for snakes (e.g., crustaceans) ( Figure 1, Supplementary Material 1 for references and details on the attribution of prey shape). If several prey types are present in the diet, the favorite items are indicated by + or ++, and their shape was used to define the "prey shape" for the species in our analysis. If no preference was noted, the shape of the prey item that requires the most extensive manipulation is considered (e.g., "fish, amphibians": amphibians more constrained because of gape limitations: bulky).

| Geometric morphometrics
We created a template consisting of a set of 921 landmarks with 10 anatomical landmarks, 74 semilandmarks on curves corresponding to anatomical features, and 837 surface semilandmarks ( Figure 2).
We manually placed all the landmarks on the template (anatomical, curve, and surface landmarks), and only the anatomical landmarks and curve semilandmarks on all specimens using the Landmark software (Wiley et al., 2005). We ensured the reliability and repeatability of the landmark positioning (see Supplementary Material 2).
Next, the surface semilandmarks were projected on each specimen, and both curve and surface semilandmarks were relaxed and slid by minimizing the bending energy (Gunz & Mitteroecker, 2013) using the "Morpho" package (Schlager, 2017). We then obtained a consensus shape for each species by performing a generalized Procrustes analysis (GPA) for symmetrical shapes on all the specimens of each species using the function "procSym" of the "Morpho" package (R script available in Supplementary Material 3). Finally, we performed another GPA on all the species consensus shapes using the "geomorph" package (Adams, Collyer, & Kaliontzopoulou, 2020) to ensure that all the consensus shapes are in the same morphological space. We used Procrustes coordinates as the shape variable to run the statistical analyses.

| Statistical analyses
We estimated the phylogenetic signal in the head shape of snakes by using the multivariate version of the Κ-statistic: Κmult (Adams, 2014a) using the "geomorph" package. The statistical significance of the K mult was obtained by running 1,000 simulations.
The K mult indicates how much closely related species resemble one another (Adams, 2014a;Blomberg, Garland, & Ives, 2003). To test the impact of diet and allometry on the head shape of snakes, we performed a phylogenetic MANCOVA using the function procD.pgls in "geomorph" (Adams, 2014b). We used the Procrustes coordinates as response variable, the prey shape as cofactor, and the log-transformed centroid size as a covariate. As the body length of the species (snout-vent length) was strongly correlated with centroid size (Pearson's correlation: df = 60, t = 9.03, p < 10 -12 , R = .75), we only used the centroid size to test for allometry. We tested for an interaction between size and diet by adding interactions to the model.
We assessed the statistical significance of the variables by using 10,000 simulated datasets obtained by permuting the phenotypic data across the tips of the phylogeny. We extracted the shapes associated with allometry (named "smaller" and "larger") by using the function shape.predictor in "geomorph" (Adams, 2014b). The shapes associated with the different dietary groups (named "bulky" and "elongated") were obtained by performing a GPA on the species belonging to each dietary group. We extracted the resulting consensus along with their centroid sizes. Then, we performed another GPA on the rescaled consensus of the groups to obtain the models in the same morphospace. We then generated meshes from the different configurations using MeshLab (Cignoni et al., 2008) and compared them using the function meshDist in "Morpho." Because the shape variability might be structured by other factors than diet and size, we used an unsupervised pattern recognition based on Gaussian Mixture Modelling (GMM) implemented in the F I G U R E 1 Phylogenetic relationship (pruned from Pyron & Burbrink, 2014

| 3D models
To test our hypothesis of many-to-one mapping of form to function, we characterized the hydrodynamic profile of five head models that describe the main axes of variability in our dataset.
Thus, we chose to work on the extreme shapes described by the first two PCs, as these components represent 65.1% of the overall head shape variability, and the mean shape (Figures 3 and 4). PC1 represents more than 54.6% of the variability and separates species with long and thin heads on its negative part from species with bulkier and shorter heads on its positive part (Figure 3). PC2 represents 10.5% of the variability and separates species having a horizontally flattened head from species with a more circular head ( Figure 3).
Aquatic snakes strike at their prey with the mouth open at various angles depending on the species, ranging from 40° (Alfaro, 2002) to 80° (Herrel et al., 2008;Vincent, Herrel, & Irschick, 2005) (Supplementary Material 4). The opening of the mouth starts at the initiation of the strike (Alfaro, 2002)

| Experimental setup
To characterize the hydrodynamic profile of the models, we measured the forces opposing the impulsive motion of a snake during a frontal strike maneuver (Figure 4b, Supplementary Video 6).
We used the same protocol as in Segall et al., 2019 to be able to compare our results with theirs. The snake models were attached to the mobile part of an air-bearing rail by a force sensor (FUTEK LSB210+/-2 Lb). Consequently, when the mobile part moves, the model pushes on the sensor, which records the axial force applied (Figure 4b,c). To mimic a strike, we positioned two springs on each side of the mobile part of the rail that were manually compressed against a vertical plate and then suddenly released, producing the impulsive acceleration. We applied different compressions to the spring to generate a range of strike velocities and accelerations. We set a position sensor (optoNCDT1420, Micro-Epsilon) at the end of the track to record the position of the cart, and calculated the kinematics (i.e., velocity U (t) and the acceleration a (t) ) of each strike by derivation of the position using Equations (1) and (2) (Figure 4b,c).
(1) We filtered x (t) and U (t) using a moving average filter over 50 datapoints.
Both force and position sensors were synchronized and recorded at a frequency of 1 kHz. We performed approximately 60 trials for each model.

| Drag and added mass coefficients
Any object accelerated in a fluid undergoes three forces that oppose the motion: the steady drag (F d ), the acceleration reaction (F a ), and the solid inertia of the body (Brennen, 1982). The force F measured on our model by the sensor is the resulting force of these three components and can be expressed as follows Vogel, 1994): dt ,

F I G U R E 3
Scatter plot of the principal components one and two (PC1 & PC2) representing, respectively, 54.6% and 10.5% of the head shape variance among the 62 aquatically foraging snake species. Dots are shaped according to the preferred prey shape (oval = elongated, square = bulky prey), and color corresponds to the centroid size of the species in mm (color scale up-left corner). Shape variation represented by each PC is shown by the two extreme shapes superimposed at the bottom (red: PC1min, blue: PC1max) and on the left (pink: PC2min, yellow: PC2max) of the figure. The phylogenetic link between species represented by the lines was generated using the function phylomorphospace in "phytools" (Revell, 2012  where ρ is the density of water, U (t) the velocity at the instant of interest, a (t) is the acceleration of the strike, and S, m, and V, are, respectively, the projected frontal surface area, the mass and the volume of the models (Table 1), and C d , C a are, respectively, the drag and added mass coefficients.
We calculated the drag coefficient C d of each model by solving Equation (4) when the acceleration is null and U = U max . When a = 0, the force measured by the sensor is only the steady drag; thus, F = F d .
The force reaches a plateau, but as the signal is oscillating, we took the average value of this plateau as a measure of the steady drag force F d (Figure 4c). Then, we calculated the drag coefficient (C d ): The term 2F d /ρS was plotted against U 2 , and the linear regression coefficient corresponds to the drag coefficient of the models (Supplementary Material 7). This representation allows to visualize the experimental data and to check the consistency of the measurement. The Reynolds number range of our experiments is 10 4 -10 5 which is consistent with previous observations (Webb, 1988).
The added mass coefficient of each model, C a , was calculated at instant t when a = a max as it corresponds to F (t) = F max : where F d(t) is the instantaneous drag. We named the numerator of Equation (7): F M , such that: . We plotted F M /ρV, against the acceleration a so the linear regression coefficient corresponds to the added mass coefficient of the models (Supplementary 8).
As these two hydrodynamic coefficients (i.e., added mass and drag) are independent of the size of the object, they are the hydrodynamic properties of the shapes, and thus, we used them as indicators of the hydrodynamic efficiency of each shape.

| Morphometric analyses
The head shape of snakes showed a significant phylogenetic signal (p = .001, K mult = 0.37). Prey shape, size, and the interaction between the two factors show a significant impact on head shape (D-PGLS: P prey = 0.008, P size = 0.002, P prey*size = 0.004). Allometry and diet, respectively, represent 7.6% and 5.6% of the overall variation in our dataset (D-PGLS R 2 ).
In snake species that prefer bulky prey, the region between the eyes and the snout is enlarged compared to elongated prey-eating snakes ( Figure 5). The upper jaw is slightly enlarged at its rear part for bulky-prey eaters, whereas the lower jaw appears more robust in elongated prey eaters. The rear part of the head is enlarged in elongated prey eaters, especially on the sides, resulting in a more tubular shape, while the bulky-prey eaters show a reduction in the head girth in this region. The eyes of elongated prey eaters are also smaller ( Figure 5).
The shape variation due to allometry is characterized by larger species having an elongated snout and a smaller head-neck transition area, which gives them an overall slenderer head compared to smaller species (Figure 6). The rear part of the head in smaller species is bulkier, whereas the front part is narrower, providing them with a head shape that is more triangular. The upper jaw is wider at its rear in larger species, whereas the mandible is bulkier and shorter in smaller species. The eyes of smaller species are also smaller ( Figure 6). The shape variation range explained by diet is smaller than the variation explained by the allometry (Figure 5a and Figure 6a, scale values).
Bulky-prey eaters have a wider range of head sizes than elongated prey eaters but overall, snake species that specialize in elongated prey have smaller heads (Figure 7). The interaction between size and dietary preference highlighted in the linear model suggests that elongated prey-eating species have smaller heads and a shape that is a combination between Figure 5c and Figure 6c. The Gaussian mixture model applied to 90% of the variability (i.e., 7 first principal components) returned a unique component suggesting little or no structure described by mixtures of datasets with normal distributions: The species in our dataset form a single cluster, which is visible in the morphospace (Figure 3). The variability in head shape is mostly carried by outlier species.

TA B L E 1 Characteristics of each model
The shapes representing the maxima of the two PCs (PC1max, PC2max) have a smaller added mass coefficient (C a ) and a smaller drag coefficient (C d ) than the shapes corresponding to the minima (PC1min, PC2min) (Figure 8). The hydrodynamic coefficients of PCmax are close to those of the typical aquatic snake profile (shape resulting from a linear discriminant analysis ). In contrast, the C a of PCmin are close to the ones of nonaquatically foraging snakes but their C d is smaller. The mean shape occupies a special position in Figure 8 by having a small C a , close to the one of PC1max and PC2max, but an intermediate C d .

| D ISCUSS I ON
In the present study, we investigated the structure and the functional implications of morphological variability in a group of species that face strong environmental constraints. First, we looked at the relationship between morphology and functionally relevant biological traits (i.e., diet and size) and we demonstrated that only a small part of the shape variability is explained by the considered factors.
Diet and size contribute to the morphological variation to a different extent, size having a somewhat larger impact on shape than the type of prey eaten. The impact of the interaction between size and diet on the head shape is not easy to interpret, but elongated prey eaters tends to have smaller heads, while bulky-prey eater shows a broader F I G U R E 5 Shape variation between bulky versus elongated prey-eating snakes. a. heatmap of the deformation from the bulky to the elongated preyeating snake in percentage of the head length (%HL) of the bulky model: Dark blue shows areas where the elongated eater is smaller than the bulky one, and red shows areas where the elongated model is larger, b. superimposition of the two shapes, c. elongated prey eater shape, and d. bulky-prey eaters shape prey shape deformation bulky prey eaters elongated prey eaters superimposition range of sizes. The deformation patterns associated with diet and allometry (Figure 5a,b. and Figure 6a Ng, 2002). We here highlight evolutionary allometry in the head of aquatic snakes showing that smaller species have a more triangular head. The same allometric pattern has been found by Vincent and colleagues in an ontogenetic context (Vincent, Vincent, Irschick, & Rossell, 2006). Thus, it is possible that small snakes benefit from having shorter, bulkier heads. However, Vincent and colleagues did not report any advantage of this shape in terms of performance. We argue that some aspect of the performance measurements or the experimental design might be improved in order to investigate into greater detail the potential advantage of this shape. For instance, a functional advantage might be revealed by including prey that differ in length. Indeed, head shape might impact the overall transport time. This bulkier shape could also increase the gape size of smallheaded snakes, yet this remains to be measured directly (King, 2002).
An increase in gape size may allow small-headed snakes to feed on a broader range of prey sizes. Overall, diet and allometry explain a relatively small amount of the shape disparity. Finally, our dietary categories (bulky vs elongated) should be taken with caution as we estimated the shape of the main item in the diet of our species based on the literature. Yet, whereas the literature is rich for some species, it is scarce for others. Moreover, the functional relevance of the categories we used remains to be tested using feeding performance measurements or models.
The phylogenetic signal in our dataset is <1, suggesting that there is more variation in the head shape of snakes than expected under a pure Brownian Motion model of evolution (Adams, 2014a).
The biological significance of K mult lower than 1 has been of interest to the scientific community. Early hypotheses suggested it could be the result of selection, adaptation, or measurement error (Blomberg et al., 2003). A recent review by Adams and Collyer (2019) suggests that K mult < 1 might be obtained when the phylogenetic signal is concentrated in a subset of traits. Yet, to our knowledge, no procedure to test for this hypothesis is currently available (Adams & Collyer, 2019). Revell and colleagues simulated several evolutionary scenarios using univariate traits and demonstrated that a low phylogenetic signal might be produced by very different processes such as divergent or stabilizing selection (Revell, Harmon, & Collar, 2008).

Considering that previous studies by Sanders and colleagues and
Watanabe and colleagues, respectively, demonstrated high rates of evolution in sea snakes, and snakes more generally (Sanders et al., 2013;Watanabe et al., 2019), one of the most likely scenarios involves strong functional constraints and high rates of evolution.
The head of snakes plays a central role in multiple fitness-related functions such as protection of the brain and sensory organs, food acquisition and manipulation, and defense against predators. It is possible that the different functions are associated with different optimal phenotypes (Shoval et al., 2012). In that case, the optimization of one or the other function would depend on the environmental pressures with which animals are confronted (e.g., more predator, lesser density preys…). Testing this hypothesis would require detailed quantitative ecological, behavioral, and functional data that are unfortunately not available to date. Finally, while the work by Revell and colleagues provides an excellent overview of how evolutionary processes and associated parameters (e.g., mutation rate) can impact the phylogenetic signal, these simulations are based on single-peak optimum, a condition that is violated by any many-toone mapping and multifunctional systems, that could potentially involve several optimum peaks (Shoval et al., 2012).
The head of snakes fulfills many functions, one of which is to catch prey. We measured the hydrodynamic constraints that resist the forward attack of a snake under water. The higher the constrains, the higher energetic cost for the animal to move its head forward (Vogel, 1994)  Our results are also consistent with the drag, and added mass coefficients found in the literature for prolate spheroids (Vogel, 1994).
Drag coefficient has been calculated for a variety of other aquatic animals such as invertebrates (Alexander, 1990;Chamberlain & Westermann, 1976), fish (Webb, 1975), amphibians (Gal & Blake, 1988), turtles (Stayton, 2019), birds (Nachtigall & Bilo, 1980), and mammals (Fish, 1993(Fish, , 2000); yet, they are not directly comparable, as different reference areas (i.e., frontal surface (S in Equation (4)), wetted area or volume) can be used. Each method is relevant depending on the system: While for a duck, one would preferably use the wetted area, for a penguin one would use either the volume or the wetted area. However, for a feeding whale the frontal area might be more relevant. Furthermore, as the drag coefficient depends on the Reynolds number (Vogel, 1994), to be comparable, the drag coefficient must be calculated in the same range of Reynolds numbers.
This makes the comparison between animals difficult as both reference areas and Reynolds numbers are strongly depend on the biological model (Gazzola, Argentina, & Mahadevan, 2014). The added mass coefficient has rarely been measured for complex shapes (Chan & Kang, 2011;Lin & Liao, 2011), but it is also known to be related to the shape of the object (Vogel, 1994). Exploring the contribution of added mass during locomotion in a wider range of organisms and how this may impact their locomotion and shape evolution could open novel avenues for more energy efficient transient-motion bioinspired objects.
Our results indicate that head shape strongly impacts the drag coefficient associated with a frontal strike maneuver in aquatically foraging snakes and, to a smaller extent, the added mass coefficient.
Bulkier heads appear to have a better hydrodynamic profile than the slenderer shapes, but even the least efficient of the aquatic foragers (PC1min, PC2min) are more hydrodynamically efficient than the snakes that never forage under water (see orange dot in Figure 8).
Thus, our results invalidate the hypothesis of many-to-one mapping of form to function as an explanation of the observed morphological disparity. Yet, we demonstrate a partial many-to-one mapping that holds for some of the shapes tested (PC1min, PC2min & PC1max, PC2max) that resulted in a similar hydrodynamic performance.
Species that drive the positive part of the morphospace (Figure 3) are highly aquatic species (Hydrophiids, Homalopsids). Their short and bulky head shape is well adapted for transient motion under water, and these shapes are associated with the best hydrodynamic profile (i.e., smallest drag and added mass coefficients). The other negative extremes of the morphospace are represented by occasional aquatically feeding species or semiterrestrial species with a long and thin head shape. These shapes are associated with the largest drag and added mass coefficients of all shapes. As these species need to be efficient on both land and under water, this less hydrodynamic profile is not surprising. The long and thin head with larger eyes might allow them to have a larger binocular field of vision and thus to be able to target their prey more accurately, whereas more aquatic snakes might not primarily rely on visual cues and thus show a reduced eye size (Hibbitts & Fitzgerald, 2005 question-oriented to fit this study, it provides a solid physical base for further model-based work. Our biological model, the head of snake, is promising from both an evolutionary and a functional perspective. The multifunctional nature of the head of snakes and their ecological diversity imposes many mechanical and ecological constraints that act together in shaping the head of these animals. Future work should include the development of feeding models in order to measure performance related to food acquisition and swallowing in snakes, as feeding is probably the more constrained and fitness-relevant activity of a snakes' head. Developing such a model, combined with Computational Fluid Dynamic models, could allow the use of performance surfaces (Stayton, 2019) and may thus offer a more thorough understanding of the phenotypic disparity of the head in snakes and its relationship with functional demands.
Ultimately, such an approach should help in untangling the interplay between different selective pressures and phenotypic responses and the mechanisms that are at the origin of evolutionary processes such as invasion of new media, adaptation to new niches.

ACK N OWLED G M ENTS
We

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.