Five specimens representing three different extant Gallotia species were analysed (Table 1), including individuals from giant and non-giant species within the genus, with ontogenetic variation present in G. stehlini.

All five specimens were recovered already dead in the field as a consequence of bird predation (shrikes for G. atlantica and juvenile G. stehlini specimens and kestrels for the G. galloti specimen), except for the adult G. stehlini specimen, which was provided already dead by the Tafira Wildlife Rehabilitation Centre (Gran Canaria Island, Spain), with an unknown cause of death. The G. atlantica and juvenile G. stehlini specimens were exposed before collection (post-mortem interval) for a maximum of 2 weeks, since regular culls were carried out. For the G. galloti specimen, the post-mortem interval was approximately 3–4 days, while for the adult G. stehlini specimen, collection occurred a few hours after death. The sex of the specimens is unknown.

Once recovered, all the specimens were immersed in water for 1 h to facilitate the extraction of skin and soft tissue. These individuals were subsequently used exclusively for nanoindentation analyses, with material sampling carried out independently from the specimens used in biomechanical modelling to preserve their integrity.

Due to the complex geometry of the skull and the need to access specific anatomical regions, each skull was fully embedded in a transparent cold-curing epoxy resin (Struers Epofix Resin). This avoided exposure to heat or pressure and ensured the production of a planoparallel block. The embedding process required approximately 8 h for complete curing of the resin.

Skull samples from the different species were fully embedded, and the targeted regions were accessed through sequential polishing. Polishing was performed starting from the upper surface of the skull, progressively exposing the areas of interest for analysis. Polishing was conducted using Struers abrasive paper with a grit size of P240 until the region of interest was exposed. Subsequently, standard polishing was performed with 1200-grit SiC paper, followed by fine polishing using diamond suspensions of 3 µm and 1 µm. Finally, an OP-S colloidal silica suspension with a particle size of 40 nm was applied to obtain a mirror-like finish. Although surface roughness was not quantitatively measured, polishing was continued until the desired surface appearance was achieved, ensuring that any residual surface irregularities did not compromise the reliability of the nanoindentation results (Suppl. material 9: fig. S1).

To investigate material properties, extant specimens were analysed using nanoindentation to assess the evolution of the contact response during loading and unloading, from which Young’s modulus (E) was obtained following the Oliver–Pharr method (Oliver and Pharr 1992; Oliver and Pharr 2004), providing an indication of the elastic response of the material. The nanoindenter, equipped with a Berkovich pyramidal tip, was used to obtain several measurements across each selected region of interest. Nanoindentation arrays were subsequently performed on 11 areas of interest across the cranium and mandible, as described in Table 2, including seven regions of bone (B1–B7), two cranial sutures (S8, S9), and two tooth-related areas (D10, dentine, and E11, enamel). Tests were performed under load-controlled conditions, with the maximum load set to 2 mN. A high-speed indentation system capable of generating material property maps was used. Nanoindentation mappings were performed using a Berkovich tip at a test frequency (load–hold–unload) of 100 kHz. At the specified load, the indentation depth and contact stiffness were recorded, enabling the calculation of Young’s modulus using the Oliver–Pharr method (Oliver and Pharr 1992; Oliver and Pharr 2004).

Given the size of the regions analysed, the spacing between indentations within each array was set to 8 µm. This ensured complete coverage of the area while preventing interaction between the stress fields generated by adjacent indentations. Considering the Berkovich tip geometry used in the tests, the relationship between penetration depth and indentation size was evaluated to confirm that the selected spacing avoided mutual influence between neighbouring indents, thereby ensuring the reliability of the measurements. Achieving the 8 µm spacing required an applied load of 2 mN. In each area of interest, the dimensions of the arrays were adjusted to the specimen size (see Suppl. material 1), and the average E values were derived from the resulting maps (Suppl. materials 2–6). A Shapiro–Wilk test was conducted on the E measurements, rejecting the assumption of normality (p < 0.05). Consequently, median E values and mean absolute deviation (MAD) were computed for each location to provide a robust representation of both central tendency and variability (Figs 2, 3, Suppl. material 9: fig. S3; Table 4).

Nanoindentation analyses were conducted under a legal permit issued by the Canary Government (Exp. 33-2023-0524072603).

Young’s modulus of the cranial bone areas of interest (E) was analysed to test whether there were differences in E among the seven cranial areas (B1–B7) and to assess whether dorsal cranial length (DCL) was a reliable predictor of E for extrapolation to fossil specimens of the same genus, Gallotia. Because E measurements constitute pseudoreplicates (multiple measurements per cranial area) nested within cranial areas and specimens, linear mixed-effects models (LMMs) were applied to account for the hierarchical, non-independent structure of the data.

Prior to conducting the LMMs, outliers for each cranial area and specimen were identified visually using Q–Q plots and boxplots generated with ggplot2 v4.0.0 in R (Wickham 2016; R Core Team 2018) and quantitatively using Tukey's 1.5 × IQR method (Tukey 1977; Hoaglin et al. 1986). Descriptive statistics were calculated with psych v2.5.6 in R (Revelle 2025) on datasets both including and excluding these outliers to ensure that the removal of these atypical values would not negatively affect the analyses. The 95% confidence intervals (CIs) for the mean E value of each extant specimen were calculated from the observed data using standard t-based intervals (the dataset is available in Suppl. material 7, and the full scripts for descriptive statistics and LMMs are available as Suppl. material 8).

First, to test for differences in E among cranial areas (B1–B7), only adult individuals (n = 3) were included in the LMM (LMM 1) to ensure a balanced design and to allow the fixed effect to be estimated independently. In this model, “E” was the dependent variable, and “Cranial Area” was included as a fixed effect, whereas the random effect “1|Specimen:CranialArea” accounted for the non-independence and clustering of repeated measurements, resulting in the following model structure: E ~ CranialArea + (1|Specimen:CranialArea).

To study the relationship between E and DCL, an additional LMM (LMM 2) was fitted, including DCL, Cranial Area, and Ontogeny as fixed effects and Specimen as a random effect, resulting in the model structure “E ~ DCL + CranialArea + Ontogeny + (1|Specimen:CranialArea)”. This second LMM included juveniles of G. stehlini (n = 2) alongside all adults to introduce sufficient variation in DCL, allowing estimation of its fixed effect on E while controlling for ontogeny, adjusting for any variation in E among cranial areas, and accounting for specimen-level random variability and nested repeated measurements of E.

The inclusion of juveniles increased the overall variability in DCL between specimens, since cranial length varies across species in adults and across ontogenetic stages in G. stehlini in this dataset. This greater between-specimen variation improves the estimation of the fixed effect of DCL and helps mitigate potential collinearity with other predictors. However, due to sample limitations, the fixed effect Ontogeny is partially confounded with Specimens/Species; therefore, the ontogenetic effect in this model must be interpreted with caution, as it may also reflect interspecific or specimen-specific differences rather than true developmental effects, a common challenge in studies of ontogenetic allometry when sample sizes and ontogenetic data are limited (Klingenberg 1998; Brown and Vavrek 2015). Because of the limited sample size and the purpose of the LMM, interaction terms between fixed effects (e.g. Cranial Area × Ontogeny) were not included in this LMM to avoid overparameterisation.

Model assumptions were verified by inspecting the residuals of the LMMs, and model diagnostics were performed using the R packages lme4 v1.1-37 (Bates et al. 2015) and performance v0.15.2 (Lüdecke et al. 2021). Normality was assessed using Q–Q plots and histograms of residuals, whereas homoscedasticity was evaluated by plotting the residuals against fitted values. Variance components for random effects indicated meaningful variation among specimens and cranial areas; therefore, the random effects in both LMMs correctly captured the nested structure of the data and pseudoreplicates. Collinearity among the fixed effects was inspected in LMM 2, which included three fixed effects, and the results of the variance inflation factors (VIFs) indicated low collinearity (VIF < 1.2). Since no substantial deviations from model assumptions were detected, no data transformations or additional procedures were required to fit the models. Both LMMs were adjusted using restricted maximum likelihood (REML) with the lme4 v1.1-37 package in R. To evaluate the significance of fixed effects in each LMM, type III ANOVAs with Satterthwaite approximation for degrees of freedom were performed using the lmerTest v3.1-3 package in R (Kuznetsova et al. 2017). In LMM 1, this ANOVA tested for differences in E among cranial areas, accounting for the random effect “Specimen:CranialArea”, while in LMM 2, it assessed the relationship between E and DCL while accounting for Cranial Area and Ontogeny as fixed effects and the random effect “Specimen:CranialArea”.

To estimate Young’s modulus (E) of fossil Gallotia specimens, the second LMM (LMM 2) fitted to the extant Gallotia dataset was used to predict E values for each cranial area based on the fossil’s DCL and ontogenetic stage (adults). Predictions were obtained using only the fixed effects of LMM 2, thereby obtaining E values representing a typical Gallotia specimen with ggeffects v2.3.1 (Lüdecke 2018) and the stats v4.5.1 base package in R. This procedure allowed the extrapolation of E values for extinct Gallotia taxa (G. goliath and G. auaritae) from the allometric relationship observed between DCL and E in extant Gallotia specimens. The 95% CIs for extrapolated E values in fossils were estimated via parametric bootstrap (1,000 iterations) using the bootMer function from the lme4 v1.1-37 package in R. In each bootstrap iteration, model parameters were resampled from their estimated distribution, and E values were subsequently predicted for each cranial area and averaged to obtain a mean cranial E value for each fossil specimen. The 95% CIs were calculated for both individual cranial areas and mean cranial E values as the 2.5^th and 97.5^th percentiles of the corresponding bootstrap distributions.

Finite element analysis (FEA) was conducted to evaluate biomechanical performance during biting, with a primary focus on comparing the effects of different mechanical properties of the bone. Models were generated for three extant Gallotia species and one extinct species (G. goliath), resulting in a total of seven models: four models using homogeneous material properties and three other models using heterogeneous properties, except for G. goliath, for which only homogeneous properties could be assigned due to the lack of region-specific material property data. The models were designed to include muscular forces as intrinsic scenarios, which enabled the generation of biting models of the lizards.

Specimens were micro-CT scanned at the Centro Nacional de Investigación sobre la Evolución Humana (CENIEH, Burgos, Spain) using a GE Phoenix V|Tome|Xs240 mCT scanner (see Table 3 for specimen and scanning details). Scanning parameters were adjusted to maximise the resolution of the CT images according to specimen size. In particular, the fossil specimen G. goliath DZUL24553 required different scanning settings due to its large size and fossil preservation compared with the extant specimens. Raw data were processed in Avizo v.2023 (Thermo Fisher Scientific) to segment and generate 3D surfaces. After segmentation, a CAD model of each specimen was created using Geomagic Wrap. During this last step, surface irregularities were repaired using refinement and smoothing tools, and the models were imported into ANSYS 2024 R1 software to perform FEA (Marcé-Nogué et al. 2011). See Suppl. material 9: text S1 for additional information. A static structural analysis was performed on the skulls of the four different Gallotia specimens using the finite element package ANSYS 2024 R1 on an HP Z6 G4 Workstation with 96 GB (8 cores × 12 GB) and 1.90 GHz. All the models were meshed with an adaptive mesh of tetrahedral elements (Marcé-Nogué et al. 2011) to appropriately capture the distribution patterns of stress and strain. Depending on the model, the meshes contained between 1.5 and 2.5 million nodes (Table 3).

Mechanical response was assessed using von Mises stress and principal strains (maximum and minimum). Von Mises stress was selected as an adequate criterion for predicting the yield of ductile materials when isotropic material properties are used in cortical bone (Doblaré et al. 2004). Principal strains, on the other hand, were used to evaluate the deformation response of the skull under mechanical loads resulting from differences in material stiffness. Therefore, they serve as a good indicator for assessing the influence of heterogeneity and homogeneity in the distribution of E (Keyak et al. 1998; Cowin 2001; Gómez-Benito et al. 2005).

Elastic and linear material properties were assumed for the bone. As mentioned above, two different approaches were used in all the models to compare the differences between a homogeneous material and a heterogeneous material. In the first approach, homogeneous properties were assumed for the bone using the following values for a reptile skull: E as the average value of all the bony locations (Table 4) and Poisson’s ratio v = 0.3. The second approach used the different values from the nanoindentations to recreate a heterogeneous material (Table 4). The E values were assigned to the same points of the models where they were analysed. Using the thermal diffusion method implemented in ANSYS, the experimentally measured E values were interpolated across the skull to generate a continuous spatial distribution of material properties (Krings et al. 2021). In this context, the thermal field was not used to represent heat transfer but rather as a numerical proxy to smoothly diffuse the elastic values between adjacent elements, minimising abrupt transitions that could cause artificial stress concentrations. This approach allows a more realistic representation of the gradual variation in bone stiffness across cranial regions. Sutures, dentine, and enamel were excluded from the diffusion process to facilitate comparison with the homogeneous models and to specifically assess the influence of interspecific differences in bone material properties.

A simple loading scenario was analysed, considering an intrinsic case of bilateral biting on both the left and right sides of the skull in the rostral position. To simulate biting, fixed vertical displacements at the tooth positions were included to replicate prey capture. Movements of the skull were restricted using joints that allow rotation but not displacement, simulating, first, the jaw-joint contact and, second, the condyle and the presence of the vertebral column (Fig. 1).

Muscle forces were estimated using the standard biomechanical relation F = σ × PCSA, where PCSA (physiological cross-sectional area) represents the area perpendicular to the muscle fibres and correlates with the maximum isometric force a muscle can produce. PCSA was used instead of the anatomical cross-sectional area (ACSA), as it better reflects the muscle’s force-generating potential. Because detailed architectural parameters (fibre length and pennation angle) were unavailable, PCSA was approximated as the mean value between the cranial (origin) and mandibular (insertion) attachment areas for each muscle, a simplification commonly applied in comparative and palaeobiomechanical modelling (Alexander 1992).

Moreover, a uniform specific tension of σ = 0.3 MPa was assigned to all muscles, following McNeill Alexander (1992). This value is widely adopted in vertebrate biomechanics as a representative estimate of maximal contraction pressure, and it provides a consistent basis for comparing muscle forces across taxa. Although some variation may occur among species or developmental stages, specific tension is generally considered a conservative property of skeletal muscle, and using a uniform value helps to minimise additional assumptions about species-specific physiology.

Muscle forces were applied as distributed loads over their attachment areas, with directions defined by vectors connecting the centroids of each muscle’s origin and insertion sites. Ten major jaw adductor muscles were included, representing the principal contributors to biting performance in lizards (see Suppl. material 9: tables S2–S5).

Institutional abbreviations: DZULL, Department of Zoology, Universidad de La Laguna (Tenerife, Canary Islands, Spain); PCCRULL, Palaeontology Collection – Carolina Castillo Ruiz, Universidad de La Laguna (Tenerife, Canary Islands, Spain).

Anatomical abbreviations: mAMEM: m. adductor mandibulae externus medialis; mAMEP: m. adductor mandibulae externus profundus; mAMES: m. adductor mandibulae externus superficialis; mAMP: m. adductor mandibulae posterior; mDM: m. depressor mandibulae; mLPt: m. levator pterygoideus; mPPt: m. protractor pterygoideus; mPSTp: m. pseudotemporalis profundus; mPSTs: m. pseudotemporalis superficialis; mPt: m. pterygoideus.

Region E11 (enamel) demonstrates higher values and greater variability in comparison with other regions (Figs 2, 3; Table 4), indicating its rigidity and reduced deformability. Following E11, several bone regions, namely B7 (inner part of the dentary), B2 and B4 (central part of the parietal and central part of the pterygoid, respectively), and B3 (quadrate), generally exhibit higher E values than D10 (dentine), whereas B6 (coronoid) and B1 (premaxilla), along with B5 (orbital osteoderm), are among the bone regions with the lowest E values. In contrast, the sutural regions S8 (contact area of the frontals and parietal) and S9 (vomer–palatine suture) consistently exhibited the lowest E values across all specimens (Figs 2, 3, Suppl. material 9: fig. S3; Table 4), indicating that they are the most flexible cranial elements.

Despite the observed differences in stiffness values among the bone regions, the range between the largest and smallest values (13.87–28.46 GPa; Table 4) suggests moderate variability overall, which can still be considered relatively low in the context of biological tissues and thus supports the assumption of acceptable homogeneity.

At the specimen level, mean E values ranged from 16.88 GPa to 24.31 GPa (Table 4). The smallest specimen, with a skull length of 18.56 mm (Table 1) and identified by the label code 020322-4 (G. atlantica), exhibits the highest stiffness, with an E value of 24.31 GPa. Following an observed trend between stiffness and size, the specimen with the label code G.galloti/21-1 (G. galloti) ranks next, with an E value of 18.97 GPa. Subsequently, the adult and largest G. stehlini, with a skull length of 61.10 mm, has a stiffness value of 16.96 GPa, while juveniles of the same species display an average of 17.04 GPa. The smallest juvenile, G. stehlini (040322-3), with a skull length of 21.66 mm, is the most elastic specimen, with a value of 16.88 GPa (Table 4).

The summary of the first LMM, which tested whether Young’s modulus (E) differed among cranial areas (B1–B7) in adult Gallotia specimens, indicated that no individual cranial area differed significantly from the reference level (B1). The Type III ANOVA with Satterthwaite approximation confirmed that there were no significant differences in E among cranial areas (F = 0.54, p-value = 0.77). These results indicate that E values can be considered homogeneous across the examined cranial regions in adult specimens of Gallotia. A full summary and the results of LMM 1 are available in Suppl. material 7.

Results from the summary of LMM 2, which tested the relationship between Young’s modulus (E) and dorsal cranial length (DCL) while controlling for other fixed effects (“Cranial Area” and “Ontogeny”) and random effects, indicated that DCL was significantly negatively associated with E (DCL: β = −0.12, SE = 0.03, t = −3.53, p-value = 0.002). The Type III ANOVA with Satterthwaite approximation confirmed the significance of DCL (F = 12.44, p-value = 0.002). This finding suggests that specimens with larger crania exhibit lower E values than smaller specimens, supporting the use of DCL to extrapolate E values to extinct and closely related Gallotia species. LMM 2 provided both an average E value for the three fossil specimens (G. goliath DZULL 2453, G. goliath DZULL 2454, and G. auaritae PCCRULL1169) and predicted E values for the seven cranial bone regions of each specimen (Fig. 4). G. goliath DZULL 2454 had the highest predicted E value (16.19 GPa), whereas the other two specimens had similar, lower values (14.67 and 14.63 GPa). Predicted E values for each cranial bone region (B1–B7) are listed in Table 5.

The fixed effect “Cranial Area” was not significant according to the results of the Type III ANOVA (F = 2.20, p-value = 0.077), showing that there were no significant differences in E between the analysed cranial areas (B1–B7), in accordance with the previous findings from LMM 1.

As previously explained, juvenile specimens were included in LMM 2 solely to increase variation in DCL between specimens, enabling a more reliable estimation of this fixed effect, without aiming to draw strong conclusions regarding the effect of ontogeny on E, which could not be fully assessed due to the limited sample size and the restricted number of ontogenetic stages represented. Exploratory results indicate that the fixed effect “Ontogeny” on E was significant in the model (β = −2.822, SE = 1.09, t = −2.60, p-value = 0.015), according to the Type III ANOVA results (F = 6.76, p-value = 0.015), suggesting that adult individuals tend to exhibit higher E values than juveniles. However, a larger sample of both juvenile and adult specimens is required to draw robust conclusions regarding the effect of ontogenetic stage on E. The full summary and results of LMM 2 are available in Suppl. material 7.

Both heterogeneous and homogeneous material models exhibited similar patterns of stress and strain distribution, as evidenced by the von Mises stress results (Fig, 5, Suppl. material 9: fig. S4). Although slight variations in stress concentration values were observed between the two approaches, these differences were not substantial enough to indicate different mechanical behaviours based on the two materials.

There is a clear distinction between species when examining the von Mises stress and strain distributions, with G. galloti having the most resistant skull (Fig, 5, Suppl. material 9: fig. S4). G. galloti develops lower stress and strain values in comparison with the other species, with fewer areas of the cranium exhibiting high stress values. G. stehlini exhibits higher stresses and strains than G. galloti, specifically in the central region of the parietal and premaxilla in the dorsal part of the cranium, but it also has one of the stiffest skulls. However, G. stehlini experiences less stress in the pterygoid than its smaller relative, while it shows a higher range of values in the lateral regions of the parietal and neurocranium.

The fossil specimen G. goliath exhibits lower stress and strain in the premaxilla and anterior part of the skull, with high stress values concentrated in the vomer, palatal bones, and lateral pterygoids, distinguishing it from G. stehlini. In contrast to its relatives, G. atlantica experiences higher overall stress, particularly in the anterior and posterior skull regions, suggesting a comparatively greater susceptibility to failure under equivalent loading conditions and making it the most fragile of the species analysed.

Nanoindentation has proved to be an effective methodology for gathering data to test the aforementioned hypotheses regarding whether a homogeneous or heterogeneous material should be used in FEA models. Moreover, the results highlight the importance of understanding material features to assess their influence on functional morphology analyses.

The studied insular species of the genus Gallotia exhibited a distinct pattern indicating a negative correlation between E and skull size. Specifically, smaller lizards showed higher median E values in the bone areas (24.31 ± 4.08 GPa), while lizards with larger skulls exhibited comparatively lower median E values (16.96 ± 3.21 GPa). Gallotia taxa demonstrate distinct skull morphological variations despite close phylogenetic relationships. These variations in the correlation between skull stiffness and individual body size may be explained by the particular environmental conditions to which each species has adapted on the various islands of the archipelago (Brown et al. 2016; Molina-Borja and Bohórquez-Alonso 2023). Each island provides a unique ecological setting where species diversify to occupy vacant ecological niches, leading to morphological and functional adaptations in response to insular conditions (Benítez-López et al. 2021). This aligns with the island rule (Whittaker and Fernández-Palacios 2007), which proposes that island species often exhibit shifts in body size, with smaller species tending to become larger and larger species becoming smaller due to ecological pressures. Such conditions may include factors such as predation level, resource availability, competition level, weather fluctuations, and geological changes such as volcanism, all of which contribute to shaping the evolutionary trajectory of each species (Vitt 2000; Meiri 2008).

Results from nanoindentation analyses suggest a relationship between size and bone stiffness. This pattern may be explained by dietary differences observed among Gallotia species, as larger species tend to maintain a more plant-based diet throughout the year, while smaller species rely more on animal-based sources, which may be scarcer and seasonally variable. Such differences in resource availability could contribute to the observed size variation among island species, in line with general patterns of insular adaptation (e.g. Masson et al. 2023; Pinho et al. 2024). However, these results also indirectly support a relationship between size and strength. Whereas stiffness is a measure of a material’s resistance to deformation, strength is its ability to withstand applied stress before failing or yielding (Cowin 2001; Currey 2006). A high E does not necessarily indicate a higher capacity to resist fracture or yielding. However, in cortical bone, different studies have suggested that strength is generally proportional to E because both properties are influenced by the bone's composition and microstructure (Doblaré et al. 2004; Currey 2006). Therefore, results from nanoindentation analysis can be interpreted as a proxy for analysing both strength and stiffness.

The stiffness of Gallotia skull bone is related to the dietary preferences of the species, as the different species under study also exhibit variations in their diet. In general terms, Gallotia taxa are classified as omnivorous (Valido and Nogales 2003), but the proportion of plant-based and animal-based food intake differs among the three analysed species. The largest species, G. stehlini, is mainly found on the island of Gran Canaria. Although considered omnivorous, its diet is primarily herbivorous, with a particular preference for leaves and seeds and a low intake of animal prey (Mateo and Lopez-Jurado 1992; Carretero et al. 2006). The results also display a consistent pattern that has been observed in previous studies (Hennekam et al. 2023), indicating that larger individuals tend to consume a higher proportion of plant matter. On the other hand, G. galloti, which has its origins on Tenerife Island but is also present on La Palma Island (Cox et al. 2010), occupies an intermediate position in terms of size when compared with its sister species. It is also classified as an omnivorous lizard, although its plant intake differs from that of G. stehlini, being more fruit- and nectar-based (Rodríguez et al. 2008; González-Castro and Siverio 2024). In addition to the effects of insularity, dietary habits are currently also shaped by agricultural practices on Tenerife. The availability of food resources introduced by plantations over recent centuries, such as tomatoes, strawberries, and grapes (Roca et al. 2005), provided a dietary alternative to seeds and leaves. Furthermore, the establishment of these plantations led to an overall reduction in food availability during winter, when animal-based sources became more important (Valido and Nogales 1994, 2003).

Lastly, the smallest-sized taxon, G. atlantica (Valido and Nogales 1994; Salvador 2015), has adapted its diet to the more arid conditions of the islands of Lanzarote and Fuerteventura. Consequently, it exhibits a higher intake of invertebrates, including hard-shelled invertebrates belonging to the orders Coleoptera, Hemiptera, and Gastropoda. Thus, this species has the least plant-based diet and consumes the toughest food items among the three studied species (Salvador 2015).

The dietary habits of vertebrates are known to directly influence cranial morphology and body size (Herrel et al. 2001; McBrayer 2004; Hennekam et al. 2023). Consequently, it is expected that a diet based on hard food resources would require stronger bone in comparison with a soft-food diet (Marcé-Nogué et al. 2017). Thus, it is not surprising that the smallest-sized species, G. atlantica, with a diet based on hard-shelled invertebrates such as snails or beetles, presents the highest E values, followed by G. galloti, which has a softer diet that includes fruit but requires additional robustness to hold and crush harder invertebrates to survive through the winter, not only because of the need to crush and ingest harder structures but also because of the need to support stronger jaw muscles. Finally, the largest living species, G. stehlini, presents an almost vegetarian diet throughout the year and presents the lowest E values. Considering the relationship between E and bone strength described above (Doblaré et al. 2004; Currey 2006), these lower values suggest a mechanically weaker skull compared with the other analysed species, whereas higher E values are consistent with stronger skulls. However, this relationship should not be interpreted as a direct indicator of biting performance, as bite force depends on multiple factors, including cranial morphology, muscle architecture, and body size (Saulnier Masson et al. 2023).

Compared with continental lizards, G. stehlini and G. galloti follow an insular trend towards increased herbivory, which has morphological implications. This trend is associated with an increase in body size, as larger lizards facilitate vegetation digestion by having a longer intestinal tract that helps absorb nutrients and process tougher foods such as seeds and leaves (Schluter 1984; Zimmerman and Tracy 1989). As body and skull size increase, the need for a robust skull decreases as the strength required to process food diminishes.

In contrast, G. atlantica, inhabiting the islands closer to the continent, exhibits a distinct pattern. This is likely due to its initial adaptation to a new niche with limited edible vegetation and a shift towards a more insectivorous diet, with higher resource availability. This dietary shift demanded the ability to handle larger prey with minimal handling time to obtain greater energetic content (McBrayer 2004; Herrel et al. 2010).

The comparatively less robust skull and higher stress concentrations observed in G. goliath suggest a dietary shift towards herbivory with increased body size related to the island rule (see below), allowing effective herbivorous feeding despite a less reinforced skull.

Dietary preferences are the main factors driving the mechanical behaviour of the jaw. In contrast, other factors also influence the protective role of the brain and sensory structures (Tseng and Flynn 2018). Therefore, the material properties obtained from nanoindentation suggest a trend more strongly associated with size and diet, although other functional influences cannot be ruled out.

Smaller species, such as G. atlantica, face higher predation pressure from a greater variety of predators. Not only is G. atlantica preyed upon by common predators shared with G. stehlini and G. galloti – such as common ravens or cats – but its smaller body size also makes it more vulnerable to smaller predators. These include small mammals such as hedgehogs and shrews, as well as smaller bird species such as shrikes or Berthelot’s pipit (Valido and Nogales 1994, 2003). Smaller species of this genus are also found on islands closer to the continent, where they face more predators and competitors due to easier island colonisation through time. The higher immigration rate of species from the mainland further intensifies these pressures.

The increased predation risk faced by smaller lizards might further contribute to the evolution of a stronger cranium by enhancing the protection of vital areas where predation pressure is high. In contrast, the results obtained for the extinct giant species G. goliath are possibly a consequence of the island rule (Palkovacs 2003; McNab 2010): these animals had very high resource availability and low or no predation and competition levels before humans settled on these islands (Castillo et al. 1994). The latter was possibly related only to birds, without the presence of cats, humans, or snakes, as occurs today.

In this regard, results from the nanoindentation tests would indicate that G. atlantica has the stiffest skull among the four species studied. However, the conclusion that G. atlantica has the stiffest skull does not necessarily imply that it is the most resistant overall, as the material properties of the skull were the sole parameters analysed, while geometry plays a crucial role in the overall mechanical resistance of the cranium and the distribution of stress, making a skull more or less resistant to failure. Accordingly, FEA reveals that, under identical biting conditions, G. galloti exhibits lower stress values than the other analysed species. Within the context of the simulated feeding scenario, this suggests a greater capacity to withstand biting-related loads without approaching structural limits. Although the models were not designed to evaluate resistance to predatory compression directly, this mechanical behaviour may also reflect comparatively greater overall cranial robustness.

This is not what would be expected if dietary preferences were the only factor shaping skull design, as G. galloti is not the species with the highest proportion of hard food in its diet. Therefore, it is suggested that protection against predation may have been an additional factor influencing the morphology of G. galloti, resulting in a comparatively more robust skull, with bony structures that may be more resistant to potential fractures. This may contribute to explaining why this species is widespread across different islands and occupies a broader range of ecological niches. In contrast, G. atlantica exhibits overall higher stress and strain values under the tested loading conditions, suggesting lower structural resistance than the other taxa. This may have implications beyond feeding performance, indicating that its evolutionary history has been shaped by multiple interacting factors, as previously discussed.

From a methodological perspective, the comparisons between homogeneous and heterogeneous material assignments showed very similar patterns of stress and strain distribution. Given that no significant differences in E were found among the cranial areas analysed, the homogeneous approach appears to be the most appropriate for the present case.