How head shape and substrate particle size affect fossorial locomotion in lizards

Substrates composed of granular media dominate the Earth and are an important habitat for many animals (European Environmental Agency, 2012; Hosoi and Goldman, 2015; Kinlaw, 1999; Post and Zobler, 2000; Sharpe et al., 2015). These granular substrates are highly variable, differing in particle size from silt to sand to gravel, particle shape from round to angular, moisture from dry to saturated, degree of compaction, and density and frictional properties (Cho et al., 2006; Li et al., 2013; McInroe et al., 2016; Mehta and Barker, 1994). As such, although the behavior of controlled and homogeneous granular substrates is reasonably well understood (Zhang and Goldman, 2014), their variation and heterogeneity results in complex behavior that is not (Li et al., 2013). The diversity of animals that use these substrates includes mollusks, various worms, fishes, amphibians, reptiles and mammals (reviewed by Dorgan, 2015; Hosoi and Goldman, 2015; Kinlaw, 1999). Granular substrates provide animals with habitat for foraging, and a refuge from heat, aridity and predators, and these factors can affect how deep and quickly they bury themselves (Arnold, 1995; Benesch and Withers, 2002; Pough et al., 1997). When animals move through these substrates, both the animal and the substrate deform and move (Aguilar and Goldman, 2015). However, our knowledge of how factors such as particle size and shape affect locomotor performance in different animals remains rudimentary.

Animals move both on top of and through granular substrates and how the substrate affects their locomotion depends on the situation and species. For example, the locomotion of some lizards is inhibited by granular substrates, likely because the substrate shifts under load (Brandt et al., 2015; Renous et al., 2008). Particle size also affects locomotion, with intermediate particle sizes, consistent with sand, allowing some lizard species to run faster than they do on finer or coarser particle substrates (Bergmann et al., 2017; Li et al., 2011). However, a force exerted on a granular substrate can also compact it, providing a solid platform from which an animal can push off, and allowing locomotor performance comparable to that on a solid substrate (Aguilar and Goldman, 2015; Bergmann and Irschick, 2010; Li et al., 2012; Mazouchova et al., 2010).

Fossorial locomotion should be more sensitive to substrate properties than locomotion on the surface because it involves displacing more substrate (Dorgan, 2015). A number of behaviors comprise fossorial movement, including sand burial, sand diving, sand swimming and burrowing or tunneling (Arnold, 1995; Ducey et al., 1993; Gans, 1974; Maladen et al., 2016). The size and shape of the animal, and the speed with which it moves affect the fossorial strategy (Hosoi and Goldman, 2015). Fossorial movement is energetically expensive, with the cost of transport being as much as 350 times that of surface locomotion (Wu et al., 2015). Substrate moisture and compaction also increase the force needed to penetrate a substrate, affecting the behavior and locomotor performance of fossorial animals (Ducey et al., 1993; Sharpe et al., 2013, 2015).

Many lizards have evolved fossoriality convergently and use their heads to penetrate the substrate (Barros et al., 2011; Baumgartner et al., 2008; Morinaga and Bergmann, 2020). Lizard heads have evolved considerable variation in shape (Stayton, 2005; Watanabe et al., 2019). Understanding how this shape variation affects fossorial locomotion requires consideration of features that are functionally relevant (Foster et al., 2018). Lizard heads are generally conical, and considerable variation has evolved in head diameter, snout pointiness and tapering of the snout (Fig. 1; Barros et al., 2011), features that physics experiments suggest would affect how easily a granular substrate is penetrated (see below). Different head shapes may also be associated with different locomotor strategies. For example, amphisbaenians use their thick heads as a shovel to move substrate laterally or vertically, depending on which plane the head is flattened in (Gans, 1974; Hohl et al., 2014), skinks such as Chalcides ocellatus use intermittent lateral oscillations of a less specialized head to wriggle into the substrate (Sharpe et al., 2015), while skinks such as Scincus scincus use high amplitude and frequency undulation of the body and a pointy snout to fluidize the substrate around them (Baumgartner et al., 2008; Maladen et al., 2009). In turn, head shape and the ability to penetrate different substrates affects habitat selection. For example, gymnophthalmid lizards that burrow in tougher soils have shorter and blunter heads (Barros et al., 2011). The amphisbaenian Trogonophis wiegmanni selects coarse sand substrates over ones with finer particles (Martín et al., 2013). The skink Lerista labialis selects less packed sand found on the crests of sand dunes (Greenville and Dickman, 2009). Thus, fossorial lizards comprise an excellent system to study how evolved variation in head shape translates into substrate penetration performance.

Although we know little about the role of head shape in facilitating fossorial locomotion, physics experiments that drive penetrators of various shapes through granular substrates provide some insight. This work has demonstrated that granular materials yield depending on their internal frictional properties and jam when a force is exerted on them by the creation and collapse of force bridges between adjacent particles, leading to fluctuations in forces experienced by the penetrator (Albert et al., 2001b, 1997, 1999; Katsugari and Durian, 2013). Objects that are narrower in diameter and streamlined require less force to penetrate granular substrates, and, therefore, give clues as to what shapes might facilitate fossorial locomotion in animals (Albert et al., 2001a; Goldman and Umbanhowar, 2008). For example, conical shapes penetrate substrates with lower force than do spheres or disks (Brzinski et al., 2013; Clark and Behringer, 2013). Drag forces experienced by the penetrator increase linearly with depth of penetration, increase with the square of their diameter, but do not change with velocity of movement, at least across biologically relevant velocities (Albert et al., 2001a, 1999; Katsugari and Durian, 2013; Stone et al., 2004). However, most of this work has used only a single, controlled substrate and non-biological penetrator shapes. A systematic manipulation of substrate properties and biological penetrator shapes is needed to understand how these factors affect fossorial locomotion in a biologically relevant context.

Here, we took a systematic experimental approach using physical models informed by evolved morphology to understand how head shape may affect fossorial performance in a variety of substrates. First, we distilled evolved head shape diversity of 152 species of fossorial lizards into three fundamental characteristics: head diameter, and the slope and pointiness of the snout (Barros et al., 2011; Storr et al., 1999). We quantified the range of variation in these metrics, and tested whether they evolved in a correlated manner in the sampled species, and whether different clades are characterized by different head shapes. We expected that narrower heads would evolve to be tapering and pointy to minimize the force needed to penetrate a substrate (Albert et al., 2001a; Goldman and Umbanhowar, 2008). However, given the diversity of strategies for fossorial locomotion and clades represented in the data, there may instead be substantial clade differences in head shape and little correlated evolution among aspects of head shape. We then used the variation in evolved head shapes to inform the construction of a series of physical models that allowed us to systematically and independently manipulate head diameter, slope and pointiness. Finally, we measured the force required for these models to penetrate a series of substrates varying from silt to sand to gravel (ISO, 2002), with both spherical (glass beads) and irregular (natural rock) particle shape. We tested the hypotheses that particle size, particle shape and head shape characteristics affect penetration force. Low penetration forces suggest adaptive head shapes for a particular substrate. We predicted that penetration forces would be lowest in particles with spherical shapes and sizes consistent with coarse sand (Bergmann et al., 2017), and for heads that are narrow, gradually tapering and pointy (Brzinski et al., 2013).

We measured the heads of 508 individuals belonging to 152 species of head-first burrowing and sand swimming lizards that we studied from natural history collections listed in the acknowledgements and in our dataset, available from the Dryad digital repository (https://doi.org/10.5061/dryad.mkkwh710d). We photographed each specimen using a Nikon D90 12.3MP (Nikon, Tokyo, Japan) camera with a Vivitar 100 mm macro lens (Sakar International, Edison, NJ, USA) from the dorsal and lateral view, and included a ruler for scale. Then we used ImageJ (http://imagej.nih.gov/ij/) to take measurements from the images to 0.01 mm precision. From the dorsal view, we measured head width, head length and internasal distance. From the lateral view, we measured head height and nasal height. Head width was the distance at the widest point of the head. Head length was the distance from the tip of the snout to the occiput of the head. Internasal distance was the distance between the external nares. Head height was the distance from the ventral to dorsal surface of the head taken at the level of the middle of the eye. Finally, nasal height was the vertical distance from the external naris to the ventral edge of the lower jaw.

We used these measurements to calculate three variables from both dorsal and lateral views: head diameter, rostral slope and rostral pointiness (Fig. 1). Head diameter was simply head width or head height, depending on the view. Rostral slope was the internasal distance or nasal height subtracted from the head diameter, all divided by the head length. The lower the rostral slope, the more gently tapering the snout. Rostral pointiness was the internasal distance or nasal height divided by the head diameter. Here, lower values indicate a pointier snout, with our index being the inverse of rostral angulation used by Barros et al. (2011). As dorsal photographs give information about the lateral characteristics of the snout and lateral photographs give information about the dorsoventral characteristics of the snout, we refer to our measurements based on what information they provide to avoid confusion. Hence, lateral rostral pointiness is internasal distance divided by head width measured from the dorsal photograph, and dorsoventral rostral pointiness is nasal height divided by head height from the lateral photograph (Fig. 1). All specimen and head shape data are available from the Dryad digital repository (https://doi.org/10.5061/dryad.mkkwh710d).

To gain an understanding of how head diameter, rostral slope and pointiness were related, we ran an evolutionary correlation analysis among these traits using code published in Morinaga and Bergmann (2017). We ran analyses in R v3.5.2 (http://www.R-project.org/), using the ‘geiger’ package (Pennell et al., 2014). Briefly, we calculated a pairwise evolutionary variance–covariance matrix using our trait data and the ultrametric phylogeny of Zheng and Wiens (2016), pruned to include the species for which we had data, and converted these to correlation matrices. To test whether evolutionary correlations were significant, we simulated 10,000 null datasets under a Brownian motion model with rate parameters matching our empirical data, and calculated P-values as the number of null datasets that had correlations greater than or equal to those from the empirical data, divided by 10,000 (Morinaga and Bergmann, 2017). We also tested whether clades of lizards represented in our dataset differed significantly in head shape variables. We did not feel it appropriate to test this using a phylogenetic ANOVA (Garland et al., 1993) because our groupings were clades, based on the phylogeny itself. In this situation, one would expect that group differences would be undetectable because a Brownian motion model of evolution results in more closely related species being more similar (Felsenstein, 1985). Therefore, we used standard ANOVA, followed by pairwise t-tests with a Bonferroni correction for multiple comparisons to identify clade differences.

We constructed 28 physical models that covered the range of evolved fossorial lizard head shapes, allowing us to systematically and independently manipulate the diameter, slope and pointiness to test the effects of these factors on penetration force in granular media. We constructed these models out of aluminium rods using a lathe to make symmetrical conical ends. These models incorporated all combinations of the maximum, minimum and an intermediate value of evolved head diameter, rostral slope and pointiness (Figs 1 and 2). One of these models also approximated the average evolved morphology (5 mm diameter, 0.55 slope, 0.33 pointiness). We tapped all models and inserted a threaded screw, allowing us to mount them on a force transducer.

To measure penetration forces of the physical models into granular media, we built a linear actuator run by a stepper motor that moved a force transducer on which we mounted each physical model (Fig. 1). We measured the force used to penetrate each medium by each model using a Kistler type 9203 force transducer, attached to a Type 5995 Charge Amplifier (Kistler Instrument Corp., Amherst, NY, USA). The custom-built linear actuator was run by a HT23-394 stepper motor (2.8 V, 2.0 A/phase), run by a 3540i Driver and PS150A24 power supply (Applied Motion Products Inc., Watsonville, CA, USA). We used Si Programmer v2.7.22 software (Applied Motion Products Inc.) to control the stepper motor from a PC, while standardizing the acceleration, velocity and penetration depth of the model.

We used silt, sand and gravel (ISO, 2002) media composed of glass beads and natural rock to test how particle size and shape affect the force required to enter the media (Table 1). We used Dragonite Solid Glass Beads (Jaygo Inc., Randolph, NJ, USA) of five different particle diameters. We obtained rock media commercially as play sand and aquarium gravel. We oven dried and sieved 64 kg of play sand using a series of screen sieves (Hubbard #548, Forestry Suppliers Inc., Jackson, MS, USA), giving us four different particle diameters. We previously characterized the physical properties of these media, shown here in Table 1 (Bergmann et al., 2017).

First, we ran an experiment to determine how the dimensions of the container used to hold the media affected penetration force. We wanted to eliminate edge effects from the walls and floor of the container because edge effects lead to increased jamming of granular media, resulting in higher and biased penetration forces (Brzinski et al., 2013; Seguin et al., 2008; Stone et al., 2004). We determined the minimum container diameter and depth needed to eliminate edge effects by manipulating these two aspects of container size and selecting the size that resulted in the minimum force for substrate penetration. To accomplish this, we constructed containers with diameters of 5, 7, 10, 15 and 20 cm out of PVC, and filled them to depths of 6, 8, 10, 12, 14, 16 and 18 cm. When doing experiments manipulating container diameter, we filled the containers to a depth of 18 cm, and when manipulating substrate depth, we used a container that was 15 cm in diameter. We used one fine and one coarse substrate composed of glass beads (particle diameter midpoints: 0.250 and 4.0 mm) and a coarse and fine substrate composed of natural rock (midpoints: 0.118 and 1.250 mm), for four substrates in total. We used a single physical model with 16 mm diameter, 0.8 slope and 0.6 pointiness because we expected a thick and blunt model to require high force to penetrate substrates and be most subject to edge effects (Stone et al., 2004).

We filled the PVC containers for all trials in a controlled manner by pouring the media slowly and evenly into the container, similar to the procedure of Seguin et al. (2008). When the medium reached the desired depth, we gently smoothed out the surface with a brush. The packing of a granular medium can be measured as the volume fraction, which is the volume of the medium relative to the volume of the occupied space (Newhall and Durian, 2003). The volume fraction can be controlled precisely with a fluidized bed (Maladen et al., 2009), which we did not have. Our approach resulted in standardized, relatively packed media that are more biologically relevant by mimicking natural dry granular substrates that are packed by gravity and vibration due to wind (Brzinski et al., 2013; Dickinson and Ward, 1994).

We conducted both edge effect trials, and trials with all 28 physical models and all nine substrates following the same procedure. We mounted each model to the force transducer, which was mounted on the linear actuator, with the container filled with substrate underneath the physical model. The computer was used to move the linear actuator until the tip of the model was 5 mm above the substrate and positioned above the center of the container. The computer then drove the model into the substrate at a velocity of 50 mm s⁻¹ to a depth of 40 mm, and recorded the maximum force. The velocity was near the upper limit of our motor, but was selected because it was comparable to published lizard fossorial locomotion speeds (Bergmann et al., 2020; Sharpe et al., 2013). The depth was also biologically relevant in being comparable to depths that many fossorial lizards bury themselves to (Sharpe et al., 2013, 2015), and it ensured that the conical portion of all heads was submerged while standardizing depth. We conducted 10 trials per model and substrate combination for a total of 2800 trials, re-pouring the substrate before every trial to avoid increased packing of the medium from multiple trials (Albert et al., 2000). All force data are available from the Dryad digital repository (https://doi.org/10.5061/dryad.mkkwh710d).

We conducted all statistical analysis in R v3.5.2 (http://www.R-project.org/). Plots of the data showed that variance of penetration force increased with its magnitude (Figs 3 and 4; Fig. S3), and this was confirmed by plots of residuals of linear models. Hence, we used generalized linear models (GLMs) with the quasi family to define the model link function and error variance (Crawley, 2012). Specifically, we used the identity and log link functions, along with error variances that were constant, as well as proportional to the response variable mean (μ), its square (μ²) and its cube (μ³). The identity link with constant error variance represents a typical linear model that assumes normality of the residuals, as implemented in the function ‘lm’. We fitted this set of models to the glass bead and natural rock data separately, and analyzed container size, physical model shape and the average model data separately. For each set of analyses, we used the residual deviance to evaluate the relative fit of the models, and calculated a Δ value similar to ΔAIC (Burnham and Anderson, 2002) to select the best model.

We conducted all analyses with penetration force as the response variable. We also normalized the penetration force by the cross-sectional area, calculated as πr² of each model, and repeated the main analyses that included all models and substrates with normalized force as the response variable because drag force scales quadratically with object diameter (Albert et al., 2001a). In analyses of container dimension data, we used particle size (fine versus coarse), container diameter and depth of the substrate, as well as two-way interactions between particle size and diameter or depth as explanatory variables. In analyses of the data from the systematic manipulation of physical models and substrates, we used particle size of the substrate; diameter, slope and pointiness of the head models; all two-way interactions, and the three-way interaction between model diameter, slope and pointiness as explanatory variables. Particle diameter was the sole explanatory variable in analyses of data using the average physical model.

We found substantial variation in head shape characteristics of burrowing lizards. Considering both lateral and dorsal aspects of the head, diameter ranged from 1.4 to 24.4 mm, rostral slope ranged from 0.225 to 0.846, and rostral pointiness ranged from 0.156 to 0.605 (Fig. 2; Fig. S1). Evolutionary correlations showed that head width, slope and pointiness were relatively strongly integrated, but when phylogeny was taken into account, these relationships were not significant (Fig. 2; Fig. S1). In particular, animals with narrower heads tended to have snouts that more gently tapered and were less pointy, but there were many species that did not fit this pattern (Fig. 2A,B). Slope and pointiness were negatively related, where species that had steeper slopes had pointier snouts (Fig. 2C). These patterns were more evident from a dorsal view of the head (Fig. 2) than from a lateral view (Fig. S1).

There were also some differences in head shape among clades. Notably, the Amphisbaenia had significantly greater head diameters than the Gymnophthalmidae and Sphenomorphinae, and also more steeply sloping and blunter snouts than some clades (Table S1). The Scincinae had significantly pointier and more steeply sloping snouts than many of the other taxa, while the Gymnophthalmidae had significantly more gently sloping snouts than most other clades (Table S1; Fig. 2).

The mechanical models that we constructed covered the full range of evolved morphologies, although some combinations of diameter, slope and pointiness were not combinations represented by evolved lizards (Fig. 2). Our model representing the average evolved morphology fell within the densest sampling of evolved head shapes (Fig. 2).

Container size data for both glass bead and rock substrates were best fitted by GLMs with an identity link function and error variance proportional to the cube of penetration force (Table S2). Although penetration forces for rock substrates were higher than for glass bead substrates, the patterns were similar (Table S3, Fig. S2). Container diameter had a large effect size for both types of substrates, with wider containers yielding the lowest penetration forces, showing a leveling off above 10 cm diameter (Fig. S2A,C). Container depth did not affect penetration force in glass bead substrates (Fig. S2B), but did in rock substrates, with the lowest force occurring at a substrate depth of 16 cm (Fig. S2D). Coarser particles required higher penetration forces among glass beads, but the difference decreased with increasing container diameter (Fig. S2A,B). Coarser substrates also often, but not always required higher penetration forces among rock substrates (Fig. S2C,D), leading to significant particle size by container diameter interactions on both substrate types (Table S3). Hence, for subsequent trials, we used a 15 cm diameter container filled with 16 cm of substrate to minimize any edge effects.

Both glass bead and rock data were best fitted by a GLM with an identity link function and error variance proportional to the cube of penetration force (Table S2). All main effects and many interactions were highly significant as a result of the high sample sizes (Tables 2 and 3), and so we relied on effect size for interpretation.

Particle size, head diameter and their interaction had the largest effect on normalized penetration force in both types of substrate (Tables 2 and 3). In general, larger particles required higher absolute and normalized penetration forces (Figs 3 and 4). However, glass beads of an intermediate size (0.625 mm) required lower absolute and normalized forces to penetrate than both smaller and larger glass beads (Figs 3 and 4). This particle size is consistent with coarse sand. This pattern was not evident with rock substrates, but this may be because we tested only four rock substrates. Both absolute (not shown) and normalized penetration forces were 2–3 times higher on rock substrates than on those composed of glass beads (Fig. 4B versus A). Head models with larger diameters required more absolute force to penetrate any of the substrates that we used (Fig. 3A), but when we normalized force by cross-sectional area, they required the lowest force per square centimeter (Fig. 3B).

Slope and pointiness of the physical models, as well as their interaction, also had medium to large effects on normalized penetration force, with the effects being more pronounced with rock substrates than with glass beads (Tables 2 and 3). These effects were most evident when visualization was restricted to only the largest diameter (16 mm) heads (Fig. 4). In general, normalized penetration forces were lower for heads that were more gradually tapered (lower slopes) and pointier (Fig. 4). On glass bead substrates, slope had a larger effect than pointiness on normalized penetration force (Fig. 4A). On rock substrates, this pattern was still evident, but less pointy and more tapering heads often required more normalized force than pointier, less tapering heads (Fig. 4B), indicating multiple ways of attaining a given normalized penetration force. This more complicated pattern for rock substrates is evident from the much larger effects of pairwise interactions between head model diameter, slope and pointiness with these substrates (Table 3) than with glass bead substrates (Table 2).

Our analyses provide some evidence that aspects of head shape evolved in a correlated manner, but there are also some important clade differences. Although head diameter, slope and pointiness were only weakly integrated in fossorial lizards, pointier snouts tended to be more steeply sloping, contrary to our predictions (Fig. 2C). Nevertheless, head shape was highly variable, and there were various clade differences. Specifically, the amphisbaenians tended to have thicker, blunter and more steeply sloping snouts, while gymnophthalmids had relatively pointy and steeply sloping snouts, and the Scincinae had relatively blunt yet gently sloping snouts (Fig. 2).

It is likely that different clades evolved different strategies for moving through granular substrates that are adaptive for particular phenotypes. For example, the blunt, robust snout with a large profile, and the thick head and body seen in the Amphisbaenia may be needed for effective lateral or vertical movement of the substrate (Gans, 1974; Hohl et al., 2014). Many amphisbaenians live deeper underground than other fossorial squamates and in moister, more compacted soils, and so their girth may be important for generating sufficient force to burrow in these environments (Gans, 1968; Navas et al., 2004). A pointy snout may be adaptive for sand diving and fast burial (Arnold, 1995; Maladen et al., 2009). Finally, a thin body with an intermediate snout may be adaptive for slower burial by secretive animals that also rely on hiding in pre-existing tunnels and crevices or digging in wet substrates (Attum et al., 2007; Sharpe et al., 2015). This last combination of traits has also been found in terrestrial caecilians (Ducey et al., 1993; Herrel and Measey, 2010).

The physical model that achieved the lowest absolute penetration forces was narrow, gradually tapering and pointy, but the model that achieved the lowest penetration forces normalized by cross-sectional area was the thickest (16 mm diameter), and also gradually tapering and pointy (Figs 2–4). Although a thicker model requiring more absolute force is expected from past physics experiments (Albert et al., 1999), this difference between absolute and normalized forces provides a basis for explaining the large diversity in the girth of fossorial species. A fossorial animal must generate the absolute force needed to penetrate a substrate, and must overcome the maximum forces experienced as a result of jamming (Albert et al., 2001a,b, 2000), based on its phenotype. Most fossorial squamates are narrow, most likely to minimize this force requirement (Fig. 2A). However, exceptions exist, with relatively thick species (>1 cm diameter) evolving in the Amphisbaenia, Acontinae and Scincinae (Fig. 2A). Our results indicate that these thicker head (and body) shapes are likely more efficient, with less force needed per unit of cross-sectional area (Fig. 3B). The cost of transport during fossorial locomotion in squamates has been estimated to be 350 times that of surface locomotion, irrespective of substrate coarseness (Wu et al., 2015). Amphisbaenians are burrowers, the most expensive form of fossorial locomotion (Navas et al., 2004), and inhabit more compacted and moist soils further underground than other fossorial squamates (Gans, 1968). Amphisbaenians are able to generate high forces during burrowing because they have axial muscles with large cross-sectional areas used in burrowing, such as the m. longissimus dorsi (Hohl et al., 2014; Navas et al., 2004). Thus, it is likely that although a larger girth requires absolutely more force to move through a substrate, it is more efficient in at least certain types of substrates.

Although thin, our best model at substrate penetration did not match a combination of slope and pointiness that has evolved (Fig. 2C). This may be because snout slope and pointiness had smaller effects on penetration force than diameter, resulting in a range of phenotypes having relatively similar performance, and, therefore, animal diameter is likely under more stringent natural selection than other aspects of head shape. Our average physical model, which does represent an evolved phenotypic combination (Fig. 2), never had the lowest penetration force, but it consistently achieved absolute penetration forces only slightly greater than those of the best models, and normalized forces comparable to those of other thin models, except in the coarsest of substrates (Fig. 3, triangles).

Slope and pointiness also influenced maximum penetration force significantly and accounted for considerable variation in the force data (Table 3). That different combinations of slope and pointiness achieved similar substrate penetration forces (Fig. 4) suggests that there are multiple ways to optimize the phenotype adaptively for fossorial locomotion. This appears to be yet another system where many-to-one mapping is at play (Alfaro et al., 2005), and may also explain why the snout slope and pointiness are related in evolved heads (Fig. 2C). However, the evolved phenotypes do not appear equivalent in performance, with steeply sloping pointy snouts requiring more force to penetrate a substrate than gently sloping snouts that are blunt (Fig. 2C). This observation may be informed by a possible tradeoff between form drag and friction drag, where a short snout may decrease friction drag as a result of decreased surface area exposed to the substrate, but increase form drag as it pushes substrate particles aside (Biewener and Patek, 2018; Zhang and Goldman, 2014).

The observed diversity of head shapes may also be explained by different kinematic strategies being required to best penetrate a particular substrate, and by the characteristics of the substrate that different animals live in. Our experiments plunged the physical models tip-first into the substrate, which is a strategy that only some animals use (Benesch and Withers, 2002; Hosoi and Goldman, 2015). Even a seemingly highly adapted snout for substrate penetration, like that of the sandfish, S. scincus, is aided by kinematic strategy. In this species, high frequency undulations of the body help to fluidize the substrate around it, lowering the necessary penetration force (Baumgartner et al., 2008). Our experiments also used a vertical entry into the substrate, yet the few fossorial squamates that have been studied enter the substrate at angles 8–31 deg below the horizontal (Sharpe et al., 2013, 2015). As penetration force increases with depth (Albert et al., 1999; Katsugari and Durian, 2013), shallow angles of substrate entry likely facilitate penetration by requiring lower forces or a slower increase in required penetration forces. The depth to which we drove the models compares well with depths that animals bury themselves to, with S. scincus submerging to depths of 2–3 cm, depending on substrate compaction (Maladen et al., 2009; Sharpe et al., 2013), and C. ocellatus submerging to depths of 1–5 cm, depending on substrate moisture content (Sharpe et al., 2015). In nature, absolute depth can also be very important in terms of escaping extreme temperatures (Benesch and Withers, 2002). It is worth keeping in mind that our physical models were conical and, although they represented the range of evolved variation, they did not reflect actual, evolved head shapes. Ultimately, it is the interplay of head shape and kinematic strategy that may best explain the diversity of evolved head shapes, and this requires more study.

Why some phenotypic combinations, particularly ones that perform well, have never evolved is another question that our data highlight. Why are there no burrowing lizards with snouts that gently taper to a sharp point or are abrupt and blunt? The latter corresponds to a performance minimum (Fig. 2C), so may be explained based on performance. Such a snout approximates the end of a cylinder, a shape associated with high penetration forces that is not adaptive for fossorial movement (Brzinski et al., 2013; Clark and Behringer, 2013). However, the former corresponds to the performance maximum, requiring the lowest penetration forces. Likewise, although thick fossorial lizards have evolved, none of these have gently tapering snouts (Fig. 2A). In both of these cases, it is possible that such a snout would be too elongate to be structurally sound. Our current data cannot address this, but finite element analysis may be able to elucidate how forces are transmitted to the skull when the length of the snout is varied.

Substrates composed of rock particles required higher absolute and normalized forces to penetrate than glass bead substrates with similar sized particles (Fig. 4). This was expected because irregularly shaped particles are more prone to jamming and have higher load-bearing capacities and angles of stability as a result of higher internal friction (Albert et al., 1997; Bergmann et al., 2017; Stone et al., 2004). Although glass beads are an artificial substrate, natural granular substrates that are highly wind or water eroded tend to be more spherical, and natural substrates exhibit considerable variation in particle sphericity (Cho et al., 2006). Therefore, our comparison of substrates is biologically relevant, and predicts that locomotion through granular substrates that are highly eroded should have a lower cost of transport, and be more highly used by fossorial animals.

We found that substrate penetration force increased with particle diameter, but we observed the lowest penetration forces on glass beads of 0.625 mm diameter, consistent with coarse sand (Fig. 3A). This same intermediate particle substrate had lower load-bearing capacity and angle of stability than finer or coarser substrates (Table 1; Bergmann et al., 2017). In addition, several species of non-fossorial lizards run faster on intermediate particle sands (Bergmann et al., 2017; Li et al., 2011). Very fine particle substrates experience cohesive forces among particles and, therefore, may be more difficult to move (Li et al., 2013), while larger particles are heavier, so more likely to jam, and also provide an irregular surface for running (Collins et al., 2013). Because fossorial locomotion is energetically expensive (Wu et al., 2015), we would expect that fossorial animals would select fine or intermediate particle substrates that are highly eroded.

The variance in penetration forces among replicate trials increased with substrate particle size (Figs 3A and 4A). Substrates with high force variance were all in the gravel range of particle sizes (ISO, 2002). It is likely that with such large particles, there are fewer spatial permutations of the relative positions of the particles, how they contact the penetrating object, and how force chains develop and collapse between them (Mueth et al., 1998). Therefore, the position of each particle has a relatively large effect on the penetration force. How relevant this is in a natural system is unclear, and may be a further reason for fossorial animals to avoid such substrates.

Although some aspects of the findings presented here on how substrate properties relate to penetration force are expected, other aspects show the complexity of the behavior of these substrates. Granular substrates vary in many ways in addition to particle size and shape, including moisture content, packing and frictional properties (Cho et al., 2006; McInroe et al., 2016; Mehta and Barker, 1994). Yet, how these factors affect how animals interact with them is largely unexplored. Various animals have habitat preferences defined by these substrate characteristics (Attum et al., 2007; Ducey et al., 1993; Greenville and Dickman, 2009; Grizante et al., 2012), suggesting that substrate variation will affect fossorial locomotion. For example, even a small amount of moisture can dramatically increase the resistive forces of a substrate by generating cohesive forces between particles (Albert et al., 1997), limiting the velocity and depth an animal can burrow (Sharpe et al., 2015). Despite this, many fossorial animals, including lizards, use moist or even wet substrates (Bergmann et al., 2020; Hosoi and Goldman, 2015; Mazouchova et al., 2010). Thus, considering how substrate variation and head shape variation interact to affect fossorial locomotor performance is key to understanding the evolution of fossorial animals and their radiation into diverse substrate niches.

We are grateful to Molly Provost for help with preliminary data collection, Arshad Kudrolli for advice on instrumentation and Joel Norton for construction of physical models. We thank José Rosado at the Museum of Comparative Zoology and Paul Doughty at the Western Australian Museum for access to specimens, as well as the California Academy of Sciences, Field Museum of Natural History, US Museum of Natural History, American Museum of Natural History and Museum of Vertebrate Zoology for the loan of specimens. We thank Henry Astley for extensive discussion and advice with the project, and Amy Cheu, Hannah Guss, Quincy Milton, Mayte Torres Zelaya and Isabel Tonelli-Sippel for feedback on earlier versions of the manuscript.