Comparative cranial biomechanics in two lizard species: impact of variation in cranial design

Lepidosaurs, and particularly lizards (i.e. non-ophidian squamates), exhibit a remarkable anatomical and ecological diversity and have been used as a model to investigate the drivers of morphological and functional variation during evolution (Evans, 2008; Herrel et al., 2007; Stayton, 2006, 2008; Watanabe et al., 2019). The considerable diversity of skull forms in lizards has been well described (e.g. Evans, 2008; Watanabe et al., 2019), but comparative data on cranial biomechanics remains limited compared with that available in mammals. Unlike the mammalian skull, where the neurosensory organs are enclosed in a shell-like bony capsule, the skull of most lizards is an open framework of bars and struts. These architectural characteristics are likely to result in important differences in the mechanical behaviour of the cranium between mammals and lepidosaurs (Curtis et al., 2011a; Porro et al., 2014; Preuschoft and Witzel, 2002; Ross et al., 2018). Investigation of the biomechanics of the cranium in lepidosaurs thus provides an alternative perspective on skull function, which is important for formulating general principles regarding the factors driving skull shape diversity across tetrapods.

Previous studies have suggested that feeding behaviour and diet have a strong influence on the evolution of cranial shape in lizards (Herrel et al., 2007; McCurry et al., 2015; Metzger and Herrel, 2005; Stayton, 2011; Watanabe et al., 2019). This link probably reflects the response of the cranium to feeding loads in some ways, as the structural organisation and material properties of bones often vary to withstand muscle loads and external forces (Meakin et al., 2014). Yet, the skull performs many functions other than feeding, including housing and protecting the brain and sensory organs, supporting the respiratory tract, and providing ornaments for sexual display. Consequently, the evolution of a complex system such as the skull appears to be driven by diverse, and potentially conflicting, demands.

Data collected on other amniotes suggest that the overall shape of the skull is not optimally designed (i.e. maximum strength for minimum material) for resisting feeding loads (Hylander and Johnson, 1997; Hylander et al., 1991; Ross and Metzger, 2004). Therefore, bone shape and bone mass distribution in the cranium do not necessarily reflect an adaptation to feeding loads (Ross, 2001). In lepidosaurs, biomechanical simulations have demonstrated the importance of certain components of the skull, such as the lower temporal bar and the quadrate–pterygoid joint, in the structural resistance of the whole system (Moazen et al., 2009a,b; Wilken et al., 2019). By contrast, other features of the lepidosaur cranium appear to have no effect on its structural resistance to external loads. For instance, the chondrocranium has little influence on the strain regime of the surrounding cranial bones in the South American tegu (Salvator merianae) during simulated bites, suggesting that this structure serves to support the brain, eyes and olfactory organs rather than to absorb and redistribute feeding loads (Jones et al., 2017).

The postorbital bar is formed by the dorsal extension of the jugal bone that connects to the postorbital or compound postorbitofrontal (Fig. 1) (Evans, 2008) and is present in most non-burrowing lizards, but has been reduced independently in two clades, Gekkota and Varanidae. In gekkotans, the postorbital bar and the supratemporal bar have been lost completely, with the jugal reduced to a small remnant in the ventral orbital margin. In varanids, the jugal is larger and extends roughly halfway up the postorbital margin, but it fails to meet the dorsal postorbitofrontal, leaving a gap of variable size. In gekkotans, loss of the bar has been linked primarily to constraints of space imposed by the increase in the size of the eye for nocturnal vision (Werner and Seifan, 2006), rather than a functional demand associated with feeding. Nonetheless, the loss of the postorbital bar has consequences for skull function during feeding by allowing a pronounced mesokinesis – movement of the snout relative to the postorbital region of the cranium – in different gecko species (Gekko gecko, Phelsuma madagascariensis, Lialis burtoni) (Herrel et al., 1999a, 2000, 2007; Montuelle and Williams, 2015; Patchell and Shine, 1986). The nature and amplitude of the intracranial movements are more contentious in varanids, and probably vary across species and during ontogeny (Metzger, 2002). Some varanid species (Varanus bengalensis, Varanus exanthematicus, Varanus niloticus) are reported to be mesokinetic (Frazzetta, 1962; Rieppel, 1978; Smith and Hylander, 1985), but others are probably not (Herrel et al., 2007; Metzger, 2002). Therefore, the presence or absence of a complete postorbital bar does not appear to be directly related to the pattern of intracranial kinesis in lizards.

How variation in cranial architecture and kinesis in lizards relates to the structural behaviour of the cranium in response to feeding loads remains unclear. It has been suggested that a complete postorbital bar increases the rigidity of the cranium by anchoring the sides of the snout to the back of the cranium (Evans, 2008; Jones et al., 2011; Porro et al., 2014; Ross et al., 2018). Therefore, the reduction of the postorbital bar might result in higher strain magnitudes in the cranial bones as a result of the bending of the snout during biting (Ross et al., 2018). Alternatively, taxa with intracranial kinesis might experience lower bone strain magnitudes, as feeding loads are dissipated by the more flexible components of the cranium (Ross et al., 2018). In this case, the structural integrity of the cranium would be preserved despite the reduction of the postorbital bar. In both varanids and gekkotans, the absence of a complete postorbital bar is also associated with an unusual frontal morphology in which the subolfactory laminae meet in the ventral midline (Fig. 1). The frontal plate forms a cylinder-like shape that might strengthen the skull while allowing the postorbital bar to be reduced (e.g. Evans, 2008). Comparison of in vivo strain gauge data obtained in different lizard species (Porro et al., 2014; Ross et al., 2018) suggests that variation in the distribution and magnitude of cranial bone strain is not obviously related to the degree of cranial kinesis, or the presence or absence of complete postorbital and supratemporal bars. As such, how the overall cranial shape and the postorbital bar, when present, affect strain magnitude and distribution in the cranial bones of lizards remains to be tested.

Computer-based biomechanical simulation techniques offer the opportunity to test in silico hypotheses of the function of biological structures. These approaches can further be used to create artificial morphologies (Gröning et al., 2013a; Lautenschlager et al., 2013; Moazen et al., 2009b; Nakashige et al., 2011; Sharp and Rich, 2016) and change the material properties of the tissues (Jones et al., 2017; Moazen et al., 2009a; Reed et al., 2011; Wilken et al., 2019) to assess the effect of a given structure in different scenarios. In the present study, we investigated cranial mechanics during feeding in two lizard species, the Argentine black and white tegu (Salvator merianae Duméril and Bibron 1839) and the African Nile monitor (Varanus niloticus Fitzinger 1826), by combining in vivo measurements with two in silico modelling techniques: multibody dynamic analysis (MDA) and finite element analysis (FEA). We used inverse dynamics in our MDA to calculate muscle activity, muscle forces, bite force and joint-reaction forces based on high-speed video recordings of the jaw movements during feeding. We then used the MDA results to define physiologically realistic boundary conditions for the FEA to calculate the strain pattern and magnitude generated by the feeding loads.

We chose the Argentine black and white tegu (S. merianae) and the Nile monitor (V. niloticus) as model organisms for the following reasons. First, both species are large (1–2 m in length) and, in the wild, are active, omnivorous hunters and scavengers. They both eat a wide variety of food materials including insects, eggs and small vertebrates, and both employ inertial feeding with larger prey items (Colli et al., 1998; Luiselli et al., 1999; Montuelle et al., 2009; Schaerlaeken et al., 2011). They thus occupy similar niches, albeit on different continents (South America, Africa). In fact, Daudin (1802) originally placed both species in the genus Tupinambis, although they are only distantly related (Tonini et al., 2016), and their respective lineages (Teiidae, Lacertoidea; Varanidae, Anguimorpha) diverged at least 150 million years ago (Burbrink et al., 2020; Pyron, 2017). Salvator merianae and V. niloticus display clear differences in their cranial shape and architecture, yet neither show any measurable mesokinesis (Herrel et al., 2007). The cranium of S. merianae is shallow and broad, with a short snout and a flat unpaired median frontal, whereas V. niloticus has a lighter, narrower cranium, an elongated snout, and a paired frontal with subolfactory laminae meeting in the ventral midline (Fig. 1). The frontals are separated along the midline by the interfrontal suture, which can fuse in old individuals. Importantly, a complete postorbital bar is present in S. merianae, but not in V. niloticus, where it is replaced dorsally by a short postorbital ligament (Fig. 1).

In this study, we used S. merianae and V. niloticus as model organisms to address the following questions. (1) Is the variation in cranial architecture observed between S. merianae and V. niloticus associated with differences in the loading regime (i.e. magnitude of bite force and muscle forces) of the cranium? (2) Are differences in the cranial architecture between S. merianae and V. niloticus related to differences in the deformation regimes of the cranial bones during biting? (3) More specifically, do the postorbital bar, present in S. merianae, and postorbital ligament, present in V. niloticus, have an impact on the pattern and magnitude of cranial bone strain? (4) Beyond obvious species-specific differences, are the overall deformation patterns of the cranial bones similar in these two species?

In vivo bite forces were measured on 63 wild and captive specimens of S. merianae and V. niloticus (Table S1). This sample includes the two specimens of S. merianae (ID: 000621516C) and V. niloticus (ID: 000617D5F1) used for biomechanical modelling (see below). The measurements were taken with a piezoelectric isometric Kistler force transducer (9311B; range ±5000 N) at the front of the jaw (Herrel et al., 1999b). The measurements at each bite position were repeated 5–10 times, and the highest measured force from those trials was retained as a measure for maximum bite performance.

We used multiple linear regression models to test the null hypothesis that mean bite force does not differ significantly between S. merianae and V. niloticus. We performed multiple linear regressions, each with one of the head dimensions (head width, head length and head depth) and the associated species as independent variables, and bite force as a dependent variable. All data were log₁₀-transformed prior to statistical analyses, which were carried out in R (http://www.R-project.org/).

The specimen (ID: 000621516C) of the Argentinean black and white tegu S. merianae (formerly Tupinambis merianae) was an adult female with the following dimensions: snout–vent length (SVL) 360 mm, head length 80.05 mm, head width 56.73 mm, head depth 45.83 mm. The specimen (ID: 000617D5F1) of the Nile monitor Varanus niloticus (formerly Varanusornatus; see Dowell et al., 2016) was an adult male with the following dimensions: SVL 435 mm, head length 81.07 mm, head width 44.34 mm, head depth 36.60 mm. Both animals were obtained through commercial dealers and housed in the Functional Morphology Laboratory, Department of Biology, University of Antwerp, Belgium, in conditions described by Ross et al. (2018). All experimental procedures were approved by the University of Antwerp Ethics Committee (reference 2006/18).

High-speed video records of the feeding events were made at the University of Antwerp, Belgium. The specimens of S. merianae and V. niloticus were filmed in lateral view while feeding. A Redlake Motion Pro 2000 digital high-speed camera (Integrated Design Tools, Inc., Tallahassee, FL, USA) attached to a Philips 14-inch image intensifier (Philips, Amsterdam, The Netherlands) was used to record the feeding events at 250 Hz. X-rays were generated using a Philips Optimus M200 X-ray generator (Montuelle et al., 2009; Schaerlaeken et al., 2011). The position of the tip of the upper and lower jaw was manually tracked in the software Tracker 5.1 (https://physlets.org/tracker/), and the gape (i.e. distance between the upper and lower jaw) calculated for each frame.

Animals were euthanised by an intramuscular injection of pentobarbital. The heads of S. merianae and V. niloticus specimens were dissected (from defrosted cadavers) and individual muscles separated. Muscles were immediately weighed after their dissection (wet mass), then placed into a 20% aqueous solution of nitric acid for 4–6 h to separate the individual muscle fibres. Nitric acid was replaced by a 50% aqueous solution of glycerol to stop the digestion, and 10–20 muscle fibres were randomly selected and photographed. The length of each fibre was then measured using the software Fiji (Schindelin et al., 2012) to calculate the average fibre length of each muscle (Table S2).

Before dissection, the defrosted heads of the specimens were scanned at the University of Hull using X-Tek HMX 160 μCT system (Nikon, X-Tek Systems Ltd). The S. merianae head was scanned to obtain an isometric voxel size of 0.1112 mm with the following parameters: beryllium target, 113 kV, 25 μA, 1000 projections, 0.1 mm copper filter. The V. niloticus head was scanned to obtain an isometric voxel size of 0.1178 mm with the following parameters: 70 kV, 17 μA, 973 projections. After reconstruction, the image stacks were saved as .tiff files and imported in Avizo 9.2.0 (FEI Visualization Sciences Group, Hillsboro, OR, USA) for segmentation. For the finite element models, four materials were manually segmented based on their density: cortical bone, trabecular bone, sutures and teeth.

3D reconstructions of the skulls obtained from the segmentation were saved as .stl files and imported into Meshmixer (Autodesk, San Rafael, CA, USA) to be altered artificially. To test the role of the postorbital bar in cranial biomechanics, it was separated from the rest of the cranium in S. merianae, and digitally sculpted and inserted for V. niloticus. The artificial bar in V. niloticus had a surface area of 178 mm² and a maximal cross-sectional area of 14.02 mm². The dimensions of the artificial bar in V. niloticus were therefore similar to those of the S. merianae postorbital bar (surface area of 155 mm², maximal cross-sectional area of 9.38 mm²). For both species, the postorbital bar was modelled as a separate segment from the rest of the cranium to test the effect of its presence and reduction on the cranial biomechanics by using the same mesh. For V. niloticus, the artificial postorbital bar was sculpted so that its extremities smoothly connected to the adjacent bones of the cranium. The surface of the postorbital bar was then imported into Avizo and converted into a 2D label that was added to the initial set of labels obtained from segmentation. This approach ensures that no artefacts are present at the boundaries between the artificial postorbital bar and the adjacent structures. The new set of labels was then used to generate a new 3D surface and then a finite element mesh of the cranium. In V. niloticus, the ventral lamina of the left and right frontal was separated from the rest of the bone and modelled as a separate material. Rendering of the surface models (Fig. 1) was performed in Blender 2.82 (https://www.blender.org/).

MDA was performed in Adams 2015 (MSC Software, Newport Beach, CA, USA). The multibody dynamic models of S. merianae and V. niloticus comprised four and six moving parts, respectively. In both models, the cranium was fixed at the level of the foramen magnum, so that the other parts could move relative to it. In the S. merianae model, the two quadrates and the two hemi-mandibles moved independently and were connected to each other by different types of joints: the hemi-mandibles were connected at the mandibular symphysis by a spherical joint with 3 rotational degrees of freedom; the quadrate–mandibular joint was defined as a hinge joint with 1 rotational degree of freedom; the quadrato-squamosal joint was defined as a spherical joint with 3 rotational degrees of freedom.

In the V. niloticus model, the two quadrates, the two hemi-mandibles and the two pterygoids could move independently. The joints were modelled as follows: the hemi-mandibles were connected at the mandibular symphysis by a spherical joint with 3 rotational degrees of freedom; the quadrate–mandibular joint was defined as a hinge joint with 1 rotational degree of freedom; the quadrato-squamosal joint was defined as a spherical joint with 3 rotational degrees of freedom; the pterygoid-basipterygoid process joint was defined as a translational joint with 1 degree of freedom. For both models, the joint types and constraints were chosen based on the joint mobility assessed during the dissection of the modelled individual, and in vivo observations available for the same species. Moving parts were imported in Parasolid format to allow for the calculation of mass and inertial properties, the latter being calculated automatically in Adams 2015 using a bone density of 1.05 g cm⁻³ (Sellers and Crompton, 2004).

A maximum fibre strength of 40 N cm⁻² was chosen for both species (Gröning et al., 2013a). The maximum muscle force was then divided by the number of strands representing the muscle in the multibody model and assigned to each strand of the muscle.

Inverse dynamic analysis was performed to calculate muscle force, the joint-reaction forces and the bite force based on kinematic data obtained from high-speed video records. The dynamic geometrical optimisation algorithm (Curtis et al., 2010) was used to simulate muscle activation dynamics during rigid-body motion. Muscle forces, joint-reaction forces and bite force corresponding to the maximum bite force for both anterior and posterior bite were exported in a format directly readable by the finite element software ANSYS v17 (ANSYS Inc., Canonsburg, PA, USA).

The finite element meshes for S. merianae and V. niloticus consisted of 7,125,144 and 7,957,252 tetrahedral elements, respectively. Adaptive meshes were generated in Avizo to guarantee the modelling of small structures, such as the sutures, while limiting the number of elements and the size of the files. Meshes were then converted to .txt format using a custom-written R script (http://www.R-project.org/) and imported into ANSYS, where the linear 4-node tetrahedral elements were converted into higher-order 10-node tetrahedral elements (ANSYS SOLID 187). Bite force, muscle forces and joint-reaction forces calculated in the multibody model were then applied to the mesh. We chose to apply the joint-reaction forces calculated with the multibody model at the tip of the quadrate instead of constraining it to avoid erroneous contact stresses in this region. The resultant sum of all applied forces was close to zero (<0.1 N), confirming equilibrium of the applied loading, but to prevent any rigid body motion of the finite element model, the neurocranium was constrained at three nodes, in all degrees of freedom, around the foramen magnum. We chose to constrain the neurocranium as it is the fixed component of the cranium in the MDA model, with respect to which the other bones are moving. FEA was run for two loading cases: an anterior bilateral bite, and a posterior unilateral bite located at roughly 70% of the out-lever length (Lappin and Jones, 2014).

Both models consisted of four materials all modelled as homogeneous, isotropic and linear elastic. Materials were assigned a Poisson's ratio of 0.3 and the following elastic modulus values: cortical bone, 17,000 MPa; trabecular bone, 560 MPa; teeth (dentine), 5500 MPa; sutures, 20 MPa. Material properties were measured with a nanoindenter (CSM Instruments S.A., Peseux, Switzerland) on the defrosted V. niloticus specimen. Because of the limitation in the scan resolution, and the computational power needed to mesh the models, sutures were artificially enlarged but remained less than 0.4 mm thick. In varanids, the anterodorsal margin of the temporal fascia is thickened and forms the postorbital ligament spanning between the postorbitofrontal and the jugal (Fig. 1). We simplified this complex morphology and modelled the postorbital ligament with a 3D spring element (ANSYS LINK 180) spanning between the postorbitofrontal and the jugal. This spring element had uniaxial tension-only capability, and an assigned cross-sectional area of 2 mm² based on measurements from the specimen used. Analyses were run for an elastic modulus of 50, 250 and 500 MPa and Poisson's ratio of 0.4. In the absence of data for lizards, we chose these values because they fall within the range of elastic modulus values reported for different ligaments in mammals (Munns et al., 1994; Nakagawa et al., 1996; Shetye et al., 2009; Stäubli et al., 1999; Vafek et al., 2018).

Rendering of the contour plots was performed in Paraview (Ahrens et al., 2005). Strain magnitude in individual bones was calculated by averaging the values collected from all the surface nodes forming the external surfaces of the bones (Tables S4 and S5). Comparisons were made between FEA results and in vivo strain gauge measurements published for Anolis equestris, Iguana, Gekko gecko and S. merianae (including the specimen here used for biomechanical modelling) (Ross et al., 2018). With respect to S. merianae, the average and maximal nodal strain recorded at each location of the strain gauge location were calculated from a series of analyses run with different loading conditions that replicated the in vivo transducer biting measurements published by Ross et al. (2018). The same was done with V. niloticus based on unpublished data that will be used for a future study. Maximal shear strain (γ_max) at each gauge location was calculated by subtracting the first and third principal strain.

Head dimensions were significant predictors of bite force (Table 1), and taxonomic group contributed significantly to the multiple linear regression model only when associated with head length. This probably reflects the difference in the relative length of the snout between the two species (Fig. 1). Our data thus support the null hypothesis that bite force does not differ between V. niloticus and S. merianae. The specimens selected for biomechanical modelling had similar maximal in vivo bite force magnitudes (Table 2), although V. niloticus showed a slightly higher bite force relative to skull width than S. merianae. The total adductor muscle mass (Table S2) was higher in S. merianae (24.90 g) than in V. niloticus (20.18 g), with the m. pterygoideus largely accounting for this difference, being about 1.6 times larger in S. merianae. When scaled to the same head width, the two species showed a similar total adductor muscle mass.

The MDA results are summarised in Table 2. The maximum bite force calculated with the MDA showed a good agreement with the in vivo bite force measurements collected for both species at the same location along the jaw, with less than 5% difference in each case. The geometry of the skull of V. niloticus appears to be better suited for transmitting muscle force to the bite point as its biting efficiency (i.e. bite force/total adductor muscle force) was 13% higher than in S. merianae during an anterior bite.

The overall strain distribution and magnitude across the entire cranium was similar between the two species but notable differences were observed (Figs 2 and 3). Average bone strain magnitude across the cranium of V. niloticus was slightly greater during an anterior bite (Table S3). Differences in bone strain magnitude were, however, particularly marked in the snout (sections 1–4, Fig. 2A), where principal strain magnitude was 49% higher in V. niloticus, and between individual bones of the cranial roof (Fig. 2; Tables S3 and S4). The two species showed similar strain distribution across the cranium (Figs 2 and 3), with an anteroposterior gradient in tensile and compressive strain and greatest magnitudes in the posterior half of the cranium. Considering the individual bones, strain magnitude in the parietal was greater than in the maxilla, and the pterygoid experienced by far the greatest strain magnitude among the bones sampled in both species. The strain gradient was lower across the cranium and between the individual bones sampled in V. niloticus compared with S. merianae (Fig. 2). Salvator merianae experienced lower strain magnitude in its snout compared with V. niloticus and showed a sharper strain gradient between the antorbital–interorbital and postorbital regions of the cranium (sections 6–7, Fig. 2A). This was reflected in the strain magnitudes in the individual bones: strain in the frontal of S. merianae was about half that in the parietal, and tensile strain magnitude was 1.6 times lower in the frontal than in the maxilla, while compressive strain magnitudes were similar (Fig. 2B; Table S4). By contrast, in V. niloticus, strain magnitudes in the frontal were similar to those in the parietal and higher than in the maxilla. However, it is important to note that fusing the interfrontal suture in V. niloticus revealed a similar pattern between the two species: lower strain magnitudes in the frontal than in the parietal and the maxilla (Fig. 2B; Table S4).

The two species showed similar deformation regimes across the cranial bones sampled (Fig. 2B). Tensile strain, however, appeared more dominant (higher |ε₁:ε₃| ratios) in the cranial bones of V. niloticus (Fig. 2B). The maxilla and the pterygoid were predominantly under tensile strain (|ε₁:ε₃|>1; Fig. 2), whereas the frontal was mostly under compressive strain (|ε₁:ε₃|<1; Fig. 2). The parietal experienced predominantly compressive strain during an anterior bilateral bite in both species, and during a posterior unilateral bite in S. merianae only. The predominant tensile strain in the parietal of V. niloticus during a posterior bite resulted in a |ε₁:ε₃| ratio close to 1 when the two loading cases are averaged (Fig. 2B).

In both species, the greatest tensile strain (>2000 με) during an anterior bilateral bite was distributed in the bones forming the cranial floor (vomer, palatines, pterygoids) (Fig. 3). Peak compressive strain (<−1500 με) was in the premaxilla, prefrontal, parietal, the supratemporal bar, the pterygoids and palatines. Varanus niloticus differs from S. merianae in having large areas of compressive strain in the maxilla and on the lateral sides of the frontals, and greater strain magnitude in its elongated premaxilla and at the base of the parasphenoid. Salvatormerianae displayed large peak compressive strain areas in the jugal and the postorbital bar.

During a posterior unilateral bite, the cranial floor of both species experienced greater tensile strain than the dorsal side of the cranium (Fig. 3). Dorsally, peak tensile strain in S. merianae was observed in the maxilla, in the frontal, the working-side supratemporal bar and in the contralateral side of the parietal, whereas in V. niloticus, peak tensile strain was solely located on the posterior half of the nasal and the left frontal. In both species, large compressive strain was found in the working-side antorbital arch, the balancing-side supratemporal bar, the contralateral side of the parietal and the working-side pterygoid. Varanusniloticus, however, displayed larger compressive strain areas in the prefrontal and in the working-side frontal. In S. merianae, the postorbital bar was predominantly under compression in all the loading cases simulated (Table 3). Compressive strain in the working-side postorbital bar was 2.5 times greater than on the balancing side during a posterior unilateral bite. Tensile strain in the left postorbital bar was, however, more dominant during an anterior bilateral bite than a posterior unilateral bite (Table 3).

The effect of the postorbital bar on the magnitude and distribution of strain in the cranial bones was more marked in S. merianae than in V. niloticus (Figs 2–4). In both species, the impact of the postorbital bar on bone strain was more marked during a posterior loading case.

In S. merianae (Figs 2–4), removing the postorbital bar during an anterior bilateral bite increased absolute peak tensile and compressive strain in localised areas, such as the supratemporal bar (tensile strain) and the prefrontal (compressive strain). In the frontal, the absence of the complete postorbital bar resulted in an important increase (>50%) in tensile strain magnitude, and a decrease in compressive strain magnitude (Figs 3 and 4). This suggests that the postorbital bar is important for resisting bending during an anterior bite. Regions of the cranium (e.g. the snout) where low absolute strain values were recorded experienced a moderate increase in strain magnitude when the postorbital bar was added. Thus, the postorbital bar not only decreases peak strain but also redistributes strain over the whole cranium in S. merianae. The effect of the postorbital bar on cranial strain was clearer during the unilateral posterior biting case (Figs 2–4). Removing the postorbital bar increased the magnitude of tensile strain in the cranial roof bones by at least 75% as well as the size of high tensile strain regions in the left premaxilla and prefrontal, the frontal and the parietal (Figs 3 and 4). Areas of high compressive strain magnitude were larger in the prefrontal, frontal, parietal, postfrontal and palatine (Figs 3 and 4). On the balancing side, the premaxilla, prefrontal and palatine experienced much greater compressive strain when the postorbital bar was absent.

In V. niloticus, the impact of a complete postorbital bar on cranial bone strain was more limited than in S. merianae (Figs 2–4). During an anterior loading case, the addition of a complete postorbital bar reduced tensile strain magnitude in the anterior aspect of the parietal by more than 50%, and more moderately in the snout and the back of the skull (Figs 3 and 4). Compressive strain magnitude was also lower in the parietal, but slightly higher in the frontal when the postorbital bar was present (Fig. 4). During a posterior loading case, the inclusion of a complete postorbital bar clearly decreased tensile strain magnitude in the nasal, prefrontal, anterior aspect of the parietal in the dorsal skull, and left vomer, palatine and pterygoid ventrally (Figs 3 and 4). The inclusion of the postorbital bar also decreased compressive strain magnitude in the dorsal surface of the frontal and in its subolfactory process (Figs 3 and 4). This suggests that the frontal morphology in V. niloticus may increase the structural resistance of the cranial roof in the absence of a complete postorbital bar.

Varanus niloticus lacks a complete postorbital bar but a postorbital ligament spans the dorsolateral gap between the jugal and the postorbitofrontal (Fig. 1). The effect of the postorbital ligament on the overall strain in the skull is minor compared with that of the postorbital bar (Fig. 5A). Increasing the stiffness of the postorbital ligament slightly decreased peak tensile and compressive strain magnitude in the frontal, nasal and anterior parietal during a posterior bite (Fig. 5B). Therefore, it is unlikely that the postorbital ligament alone fulfils the mechanical role of a postorbital bar in V. niloticus.

Altering the morphology of the frontal affected the strain magnitude in this bone but had a limited impact on the strain regime in the rest of the cranium (Fig. 6). Removing the subolfactory processes of the frontals markedly increased strain in the ventral surface of the frontal during a posterior bite (Fig. 6A), with larger peak tensile strain (>2000 με) areas on the working side frontal and larger peak compressive strain (<−1500 με) areas on its counterpart. The subolfactory processes of the frontals therefore appear to increase the structural resistance of the interorbital region of the cranial roof.

Fusing the interfrontal suture decreased tensile and compressive strain magnitude in the frontal of V. niloticus by about a half without really affecting the adjacent parietal (Figs 2 and 6B; Table S4), but strain magnitude remained higher than in the frontal of S. merianae. Strain was more evenly distributed in the frontals (Fig. 6B), and compressive strain became more dominant (|ε₁:ε₃|=0.69; Table S4). Notably, this is reflected during a posterior unilateral bite by large peak compressive strain in the balancing side subolfactory processes (Fig. 6B), which further highlights the importance of these structures for resisting feeding loads as the interfrontal suture can fuse in large adults.

We observed that bone strain magnitude was more uniform along the cranium in V. niloticus than in S. merianae. Strain magnitude was higher in the maxilla and other cranial bones of V. niloticus and in the entire anterior portion of the cranium (Fig. 2; Table S3). These differences were irrespective of the presence of a complete postorbital bar or postorbital ligament (Figs 2–5), but were most probably linked to the distinct snout form observed in these two species. Previous FEA predicted that archosaurs with fenestrated and flattened snouts (i.e. platyrostral cranium) experience higher strains and stresses than their tall and domed counterparts (i.e. oreinirostral cranium) (McHenry et al., 2006; Rayfield and Milner, 2008; Rayfield et al., 2007). With respect to squamates, FEA performed on Iguana (Simões et al., 2016), Sphenodon (Curtis et al., 2011a), Gekko (Cost et al., 2020) and Uromastyx (Moazen et al., 2009a) predicted low stress/strain in the snout. In these taxa and S. merianae, the shorter, broader and somewhat domed snout probably maximise the second moment of area and thus better resist bending (Ghavami, 2015) than the long, narrow, flatter and fenestrated snout of V. niloticus. This might be reflected by the more dominant tensile strain in the maxilla of V. niloticus. However, the elongated snout of V. niloticus and other varanids (McCurry et al., 2015) probably increases the rotational velocity at the tip of the jaw for capturing elusive prey (Herrel et al., 2007; Metzger, 2002; Metzger and Herrel, 2005; Stayton, 2005). In addition, the lighter cranium and taller postorbital region in V. niloticus may maximise the rotational velocity of the head generated by the cervical muscles as varanids, with their highly specialised lingual apparatus, appear to rely more on inertial feeding than does S. merianae (Elias et al., 2000; Montuelle et al., 2009).

The two species studied here also differ markedly in their degree of sexual dimorphism and mating behaviours. Whereas varanids show little sexual dimorphism and males engage in ritualised, wrestling-like, combat (Khan et al., 2018; Murphy and Mitchell, 1974; Tsellarius and Tsellarius, 1997), S. merianae shows strong sexual dimorphism in head form and muscle size (Fabre et al., 2014a; Naretto et al., 2014), and males engage in combat involving biting. Moreover, in male S. merianae, bite force scales disproportionately with head width compared with that of females and is associated with more aggressive behaviours in males (Herrel et al., 2009), which, as in other dimorphic lizard species (Lailvaux and Irschick, 2007; Lappin and Husak, 2005; Lappin et al., 2006), might favour success in male combat, and resource and mate defence (Naretto et al., 2014). Accidents or antagonistic interspecific and intraspecific interactions can be associated with relatively higher bone strain and injury (Cooper and Vitt, 1987; Jurmain, 1997), which impose a greater demand for bone resistance. Strain magnitude has an important impact on the mass and distribution of bone, and variation in bone form reflects the loading regime experienced during a range of habitual and infrequent behaviours and events (Ehrlich and Lanyon, 2002; Frost, 2003; Gröning et al., 2013b; Meakin et al., 2014). Although we used a female S. merianae specimen, intraspecific differences were far less pronounced than those between S. merianae and V. niloticus. Therefore, it is also possible that the agonistic behaviours between males may impose a greater demand on certain regions of the cranium in S. merianae, resulting in lower bone strain magnitudes in the frontal and the snout, and larger gradients across the cranium compared with those in V. niloticus during feeding.

Our biomechanical simulations provide a mechanistic interpretation of the pattern of co-variation between bite force, cranial shape and muscle cross-sectional area observed between males and females in S. merianae (Fabre et al., 2014a,b). Regions of the cranium, such as the postorbital portion of the cranium and the palate, whose shape strongly co-varies with bite force and muscle cross-sectional area, are predicted to experience greater bone strain in our FEA. By contrast, cranial regions that showed little shape variation, such as the nasals, are those that experienced low bone strain. The differences in the shape of certain cranial regions between males and females in S. merianae might therefore represent a response to the loading regime of the whole skull, similar to the intraspecific variation shown within marsupial species due to masticatory loading (Weisbecker et al., 2019).

We assessed the accuracy of our finite element models by comparing bone strain magnitudes calculated in a FEA series (Table 4) at individual strain gauge sites with in vivo bone strain measured in other squamates (Ross et al., 2018). Note that the values obtained from these additional analyses (Table 4) are not directly comparable with the strain magnitudes for the entire bones (Fig. 2B; Table S4). Strain calculated in our FEA falls within the range of values measured experimentally. Mean principal strain values measured in Anolis equestris, Gekko gecko, Iguana iguana, Uromastyx geyri and S. merianae ranged from 102 to 1004 με (ε₁) and −147 to −1195 με (ε₃) (Porro et al., 2013; Ross et al., 2018). Bone strain magnitudes calculated at strain gauge sites for S. merianae were underestimated when compared with published strain gauge records (Table 4) made on the same individual (with the exception of the tensile strain in the frontal), but fell within the range of those obtained on two other specimens (Ross et al., 2018). This discrepancy might be due to the fact that our simulations did not capture the whole range of loads that the cranium experiences during biting (e.g. tearing forces caused by the pull back of neck muscles, side to side shaking of relatively large prey), and/or that the restraint of the animal during transducer biting might have caused higher bone strains. With respect to V. niloticus, tensile strain (ε₁) magnitudes in the frontal bone calculated in our model (Table 4) are within the range of magnitudes (100–600 με) collected by Smith and Hylander (1985) in V. exanthematicus during feeding. For both species, most strain magnitudes obtained from the entire bones are greater than those from strain gauge sites when similar loading cases are compared (Table S5). Together with the good match between the MDA results and experimental data, these comparisons suggest that reasonable biological interpretations can be drawn from the present biomechanical models. A validation study is currently being undertaken using unpublished in vivo data to further assess the accuracy of our models and determine the key parameters that affect their output.

Whether a complete postorbital bar was included or not, we observed similarities in the overall pattern and magnitude of bone strain between S. merianae, V. niloticus and other lizards. Our results thus support the hypothesis that bone strain magnitude and distribution do not radically differ between lizard species with and without a complete postorbital bar. Bone strain is not homogeneous across the cranium, and the highest strain magnitudes are located in the circumorbital and postorbital regions, and the palate (Figs 2 and 3). The pterygoid is more highly strained than any other bone, probably because it serves as an attachment area for the m. pterygoideus – the largest jaw-closing muscle. Interestingly, these regions also show high disparity and rate of evolution in lizards (Watanabe et al., 2019). Hence, biomechanical demands appear to be reflected in the variation in the form of the cranial regions at both interspecific and intraspecific levels. The parietal experiences higher strain magnitudes than the maxilla (Fig. 2) and displays large areas of peak strain whose distribution varies with the location of the bite point along the tooth row (Fig. 3). Thus, the parietal does not simply serve as a muscle attachment area to withstand muscle loads but also resists the loads transferred from the bite point to the back of the cranium. Areas of peak strain on the parietal are reduced when sutures are fused (Jones et al., 2017; Fig. S1), which further underlines the role of the parietal in resisting biting loads and the importance of sutures for load transfer across the cranium (Curtis et al., 2013; Moazen et al., 2009a). However, we do not know the relative contribution of each suture in this phenomenon and whether the fronto-parietal suture plays as prominent a role in S. merianae and V. niloticus as it does in Uromastyx hardwickii (Moazen et al., 2009a).

Consistent with in vivo strain gauge measurements, we found higher strain magnitudes in the parietal than in the maxilla when the entire bones were sampled (Fig. 2). However, the similar or higher strain magnitudes in the entire parietal compared with the entire frontal contrast with strain gauge measurements, which consistently found the opposite pattern in Iguana, A. equestris and G. gecko (Ross et al., 2018). This discrepancy between experimental and finite element results might be because the regions on the parietal that experience the highest strain magnitudes are covered by muscles and thus cannot be sampled with in vivo strain gauges (Figs 2 and 3). Although the good match between our simulations and experimental data (Tables 2 and 4) gives us confidence in our models, the effect of modelling approximations and potential artefacts still cannot be ruled out; our ongoing validation study will thus be important to clarify this point. The higher strain magnitudes in the frontals of V. niloticus are caused by the interfrontal suture, and it is important to note that the two species showed a consistent strain pattern across the cranial bones when this suture was fused (Fig. 2). Finally, despite marked differences in the cranial shape of the two species, the individual bones sampled across the cranium showed similar deformation regimes. Beyond species-specific differences, these results therefore suggest that lizard species with different cranial shapes may share a common deformation regime, something also seen among anthropoid primates (Ross et al., 2011).

In S. merianae, the postorbital bar is among the regions of the cranium that experience the highest strain magnitudes. We found that the postorbital bar was predominantly under compressive strain, whereas strain gauge measurements suggest that tension is the dominant loading regime in the postorbital bar of U. hardwickii (Porro et al., 2014). This difference might be because the jugal serves as an attachment area for the external bundle of the m. pterygoideus in this species. It is also worth noting that, in our models of S. merianae, peak tensile strain was located on the lateral side of the postorbital bar, whereas larger, peak compressive strain was on the medial side (Fig. 3). Therefore, the different deformation regimes of the postorbital bar between S. merianae and U. hardwickii could also be explained by the placement of the strain gauges on the lateral side rather than the medial of the jugal (Porro et al., 2014). Unfortunately, FEA results reported for U. hardwickii cannot provide clear answers to this question (Moazen et al., 2009a). Strain magnitude in the postorbital bar of S. merianae increased during posterior biting compared with anterior unilateral biting, as in U. hardwickii (Moazen et al., 2008; Porro et al., 2014), and its removal increased strain magnitude in the bones of the cranial roof and palate. When present in lizards, a complete postorbital bar therefore appears to be important for maintaining the structural integrity of the cranium by reducing its bending (during anterior biting) and twisting (during unilateral biting), and by providing an anchoring strut for the muscle attachment areas in the postorbital region.

The variation in the cranial struts of the frame-like skull of lepidosaurs was hypothesised to be tightly linked to the evolution of bite performance and feeding function (Rieppel and Gronowski, 1981). Geckos, which have lost the postorbital and supratemporal bars, have a relatively lower adductor muscle mass and bite force, and a lighter cranium compared with other squamates (Herrel et al., 2007). Although based on a limited sample, we did not find significant differences in bite force between specimens of S. merianae and V. niloticus. The loading regimes (Table 2) of the crania of the two specimens used for modelling were also similar, and V. niloticus skull geometry was slightly more efficient at transmitting muscle force to the bite point. Yet, the postorbital bar, digitally added in V. niloticus, was less efficient in decreasing peak bone strain than in S. merianae. Thus, neither biting performance nor the ability of the cranium to withstand high feeding loads appears to be impaired by the reduction of the postorbital bar in varanids. This suggests that morphological and/or behavioural changes during varanid evolution might compensate for the absence of a postorbital bar or reduce the importance of a previous role. We think it unlikely that the postorbital ligament alone can fulfil the role of the postorbital bar as it appears to play a minor role in strain absorption. However, our finite element model may also not represent the soft tissue anatomy adequately enough to fully exclude this possibility. In V. niloticus, the postorbital ligament is the thickened free anterodorsal margin of a sheet of temporal fascia that stretches across the temporal region (upper and lower fenestrae) enclosing the supratemporal and postorbital bars and attaching to the rictal fold ventrally and the quadrate posteriorly. The tensioning of muscle fascia by muscle bulging was shown to decrease peak bone strain in macaques and Sphenodon (Curtis et al., 2011a,b). A similar effect might occur in V. niloticus but including the entire temporal fascia would have greatly increased the complexity of the finite element model and analyses.

From an evolutionary perspective, the drivers of the reduction of the postorbital bar in varanids remain unclear. It is possible that the reduction of the bar was originally associated with mesokinesis that has been secondarily lost in large varanids or that the reduction of the bar is a by-product of accelerated growth of the postorbitofrontal during development (Werneburg et al., 2015). The acquisition of an active foraging life-style (McBrayer, 2004) and an inertial feeding mode (Herrel et al., 2000), which entails important accelerations of the head, have been suggested to be linked to the lengthening and lightening of crania in varanids. In this regard, the frontals of V. niloticus notably appear to be better optimised for maximum strength with minimal material during biting than in S. merianae (Figs 1 and 2). However, the available data appear to contradict this hypothesis as varanids were not found to have a lower skull to body mass ratio compared with other lepidosaurs (Metzger, 2002). Male–male interactions and mating behaviour might represent another important potential driver for cranial evolution in varanids and other lepidosaurs, but the influence of these factors on cranial mechanics has never been directly assessed.

Put in a broader context, our results bring additional insights into the factors underlying the evolution of the cranial design in amniotes. During maximal biting, the maximal shear strain magnitude recorded at the working-side postorbital bar of S. merianae was at least 2.5 times higher than values reported for Eulemur, Otolemur, Aotus and Macaca (Nakashige et al., 2011; Ross and Metzger, 2004; Ross et al., 2011). In Macaca, the postorbital bar and septum appear to have little role in the structural resistance of the cranium to biting loads (Nakashige et al., 2011; Ross et al., 2011) and might instead serve for oculomotor stability (Cartmill, 1980; Heesy et al., 2007; Nakashige et al., 2011; Ravosa et al., 2000; Ross and Hylander, 1996), whereas preliminary finite element results suggest that digitally removing the postorbital bar does impact cranial bone strain in Eulemur (Strait et al., 2014). Bones forming the circumorbital region in S. merianae and V. niloticus also experienced higher peak strain magnitude than in the homologous region in mammals (Bright, 2012; Ross and Metzger, 2004; Ross et al., 2011). These differences might reflect the greater biomechanical role of the circumorbital region during biting in lizards compared with mammals.

Strain distribution across the cranial bones of S. merianae and V. niloticus is more homogeneous that in mammals. Consistent with previous observations made on lizards, the calvarial bones (parietal and frontal) of S. merianae and V. niloticus have higher strain magnitudes than those recorded in mammals. The parietal is more loaded than the maxilla in lizards, whereas it experiences lower strain than the facial bones in mammals (Behrents et al., 1978; Bright, 2012; Cox et al., 2012; Herring and Teng, 2000; Ross and Metzger, 2004; Thomason et al., 2001). Previous studies (Cox et al., 2012) and our preliminary results also suggest lower strain gradients and magnitudes in the palate of rodents and rabbits compared with values in the two lizards studied here. However, a more detailed comparison with finite element results obtained on mammals is difficult as previous analyses did not necessarily incorporate the same level of details or employ the same metrics. The combination of anatomical, developmental and biomechanical data in an explicit phylogenetic framework will be essential to better understand the determinants of the variation in the skull form during tetrapod evolution.

We used in vivo bite force measurements, high-speed X-ray videoradiography and computer-based biomechanical simulation techniques to investigate the cranial biomechanics of S. merianae and V. niloticus. The differences in the strain regimes of the cranial bones were not related to the presence of a complete postorbital bar, but rather to the distinct overall cranial architecture observed between these two species (tall and broad snout in S. merianae, long and narrow snout in V. niloticus). The postorbital bar is important for the structural resistance of the cranium to feeding loads in S. merianae, and potentially during antagonist male–male interactions, whereas the postorbital ligament probably does not have a substantial biomechanical role in V. niloticus. Our results suggest that the reduction of the postorbital bar in V. niloticus impaired neither its biting performance nor the structural resistance of the cranium to feeding loads. Beyond differences related to species-specific variation in morphology, the two species share a similar strain and deformation regime of the cranium during biting. Strain magnitude is greater in the postorbital region (specifically in the parietal and pterygoid) and the circumorbital region, which appears to be important for resisting feeding loads. This suggests that common mechanical behaviour might underlie the frame-like cranium of lizards.

We thank Sue Taft (University of Hull, UK) for her assistance in scanning the specimens and measuring the bone material properties, and the Viper High Performance Computing facility of the University of Hull and its support team for their help and assistance in running the FEAs. Alexander Blanke (University of Bonn, Germany) and Manuel Pinheiro (Ghent University, Belgium) are thanked for the valuable discussions during this project.