Adaptive or non-adaptive? Cranial evolution in a radiation of miniaturized day geckos

Abstract

Lygodactylus geckos represent a well-documented radiation of miniaturized lizards with diverse life-history traits that are widely distributed in Africa, Madagascar, and South America. The group has diversified into numerous species with high levels of morphological similarity. The evolutionary processes underlying such diversification remain enigmatic, because species live in different ecological biomes, ecoregions and microhabitats, while suggesting strikingly high levels of homoplasy. To underscore this evolutionary pattern, here we explore the shape variation of skull elements (i.e., cranium, jaw and inner ear) using 3D geometric morphometrics and phylogenetic comparative methods on computed tomography scans (CT-scan) of a sample encompassing almost all recognized taxa within Lygodactylus. The results of this work show that skull and inner ear shape variation is low (i.e., there is high overlapping on the morphospace) across geographic regions, macrohabitats and lifestyles, implying extensive homoplasy. Furthermore, we also found a strong influence of allometry shaping cranial variation both at intra and interspecific levels, suggesting a major constraint underlying skull architecture, probably as a consequence of its miniaturization. The remaining variation that is not allometric is independent of phylogeny and ecological adaptation and can probably be interpreted as the result of intrinsic developmental plasticity. This, in turn, supports the interpretation that speciation in this group is largely concordant with a non-adaptive hypothesis, which results mainly from vicariant processes.

Introduction

Understanding the processes underlying morphological variation has long been one of the major aims of evolutionary biology. Historically, phenotypic variation has been associated with selective pressures acting both on morphology and function. However, recent studies have shown that evolutionary constrains (e.g., phylogenetic inertia, allometry, integration/modularity) can yield major effects on morphological variation [1,2,3], underscoring both adaptive and non-adaptive radiations [4]. Consequently, evolutionary patterns like homoplasy may result from multiple evolutionary scenarios in which function is involved (convergent evolution), stochastic changes that result in similar forms [5], and the reuse and influence of deeply shared developmental pathways (parallelism; [6]). Additionally, the degree of functional similarity between different phenotypes may vary considerably, from complete convergence [7], to incomplete [8], or imperfect convergence [9]. Furthermore, some studies have also pointed that phenotypic plasticity can be a source of novel phenotypes, providing the grounds for speciation processes (known as ‘plasticity-first evolution’, PFE; [10]). All this diversity of evolutionary scenarios highlights the complexity and multidimensionality of the patterns and processes underlying morphological similitude across taxa/across groups.

Geckos are among the most species-rich and geographically widely distributed group of all reptiles [11, 12]. Such an exceptional diversification has also resulted in one of the most morphologically diverse group of reptiles, harboring a large variety of ecomorphs and a multiple life-history traits [13,14,15]. Geckos range from miniature sizes (like Chatogecko, Lygodactylus, Microgecko and Sphaerodactylus), to large forms like the New Caledonian giant gecko (Rhacodactylus leachianus), some member of the Malagasy geckos of the genus Phelsuma, or the extinct Gigarcanum delcourti [16]. Interestingly, however, morphological evolution of geckos, as in most groups of vertebrates, has resulted in a complex combination of plesiomorphic, apomorphic, and homoplastic characters [13, 17, 18]. While some gecko groups display ecological and morphological variation that is correlated with lifestyle [13], other groups have retained striking levels of morphological similarity, with conserved morphological patterns that are better explained by phylogenetic conservatism than by ecological adaptability [19].

The skull of gekkotans presents a unique “gecko architecture” with remarkable kinetic mobility [20, 21], which results mainly from the reduction or even loss of some of the circumorbital bones (and hence the postorbital bar) [20, 22]. Furthermore, geckos exhibit a large variety of cranial apomorphies that have been historically used in taxonomic and phylogenetic studies [24,25,26,27]. Recent studies have found a strong influence of allometry and phylogeny as evolutionary constraints shaping morphological diversity in some of the skull elements in geckos, such as the quadrate bone [28]. Additionally, some homoplastic characters have been identified in miniaturized geckos, including the increment of the overlap in the rostral portion of the craniofacial system, or the loss or reduction of temporal bones [29,30,31]. In addition, inner ear has been lately associated with ecological adaptation on lizards [32, 33] due to the relationship of the semicircular canals and vestibular systems with the locomotor activity, agility and speed [34].

The genus Lygodactylus represents a radiation of diurnal dwarf geckos that has successfully colonized different biomes and ecoregions, exploring both arboreal and rupicolous lifestyles [35, 36]. However, Lygodactylus species exhibit high levels of external morphological conservatism, with several cryptic taxa and without significant signal of morphological adaptability [19, 36,37,38,39]. Consequently, some authors have suggested that species-richness in this group is the result of non-adaptive radiation [19, 36]. Despite this, the skull of Lygodactylus is particularly interesting, given the fact that it has a unique architecture associated with a series of plesiomorphic and apomorphic character states that could be related to the miniaturization of the lineage. Such architecture includes a unique combination of characters, i.e. fused nasal bones, the reduction or absence of the postorbitofrontal, loss of the squamosal, and vestigial or lost jugal bones [31, 40].

Such morphological conservatism across a wide range of habitats and lifestyles raises questions such as if internal and external morphological characters have evolved in similar morphological directions among ecological features, or if downsizing affected similarly different species. Accordingly, here we explore skull shape diversity in Lygodactylus, testing the effect of different evolutionary drivers (i.e., allometry, geographic distribution, phylogenetic signal, macrohabitat and lifestyle) on the skull and inner ear, aiming to elucidate evolutionary patterns of phenotypic change and adaptation within this morphologically conservative lineage of dwarf geckos.

Methods

Taxon sampling and morphological data

We included representatives of 79 of the 96 described nominal taxa (including species and subspecies) within Lygodactylus, plus two candidate new species (Table S1). Also, aiming to explore intra- and interspecific variation, we included 11 additional individuals of the strictly arboreal species, Lygodactylus kibera (Table S2). Specimens selected for this work were in good condition, without broken bones or other artifacts (Fig. 1). We prioritized, when possible, type material to avoid misidentifications that can arise from the conservative features of this group (Table S1). We also recorded distribution and ecological traits for each sample, primarily from original descriptions and information recorded by [38, 39 and 41], as follows (Table S1): geographic distribution (Africa, America, Madagascar, or oceanic islands); macrohabitat (rainforest, dry forest [including dry forest/savanna/desert-dwelling], afromontane or transition zone); and lifestyle (arboreal, rupicolous or mix [opportunistic species that use arboreal and rupicolous habits indistinctly]).

Osteological data

To enable efficient sampling and optimize the quality of every scan, computed tomographic (CT) scans for every individual were generated using different sources, parameters, and protocols, at different institutions and using different CT scanners defined in Table S2. All 3D segmentation models and reconstructions were generated for articulated cranium, left mandibles and right inner ears in Avizo Lite 2020.2 (Thermo Fisher Scientific 2020). We excluded the inner ear of L. insularis and L. soutpansbergensis due to the high level of endolymphatic calcification that precluded a good reconstruction of the internal structures. To facilitate visualization, individual cranum and mandibles were colored following the same color palette as in 31. Cranium, mandible, and inner ear reconstructions were exported as mesh files (.ply). Finally, meshes were reduced in size and smoothed in MeshLab v.2022.2 [42] to facilitate posterior placement of landmarks and semi-landmarks. Tomogram series for all scans are available on MorphoSource.org (Table S2).

Landmarking

We used Stratovan Checkpoint v.2018.08.07 (Stratovan Corporation, Davis, CA, USA; https://www.stratovan.com/products/checkpoint) to digitize three 3D landmark datasets, namely, (a) the cranium (defined as the entire skull excluding the lower jaws and the inner ear), (b) the mandible (left lower jaw) and, (c) the inner ear (right). To remove asymmetry from the shape data we placed 66 fixed landmarks (Fig. S1) in homologous repeatable points on the right-side of the cranium. We then mirrored the right-side using “mirror2plane” function [43] of the R package Morpho v.2.11 [44] to create the left-sided coordinates (Fig. S3). We used 20 fixed landmarks in the mandible using also homologous and repeatable points (e.g., foramen or foramina, or sutures, Fig. S1). Finally, to record the different structures of the inner ear and the orientation and direction of the semicircular canals, we used a combination of 11 fixed landmarks and 48 sliding semi-landmarks (Fig. S2A) located along the semicircular canals. Sliding semi-landmarks were automatically subsampled by utilizing the “curve” option in Stratovan Checkpoint, i.e., adding two equidistant landmarks between every landmark, resulting in a total of 144 sliding semi-landmarks. Semi-landmarks were slid to minimize total bending energy of a thin plate spline (TPS) across all specimens using R package Morpho v.2.11 [44]. Semi-landmarks were accommodated on the different species to have the same number of landmarks equally distributed between the anterior, posterior, and horizontal semicircular canals (Fig. S2B). All landmark configurations were superimposed utilizing the Procrustes criterion [45, 46], by scaling, rotating, and aligning them to the mean using ‘gpagen’ function in R package geomorph v.4.0.1 [47] and using MorphoJ v.1.7. [[48]; Fig. S1]. Descriptions of the landmarks for the cranium and mandible are available in Tables S3 & S4.

Phylogenetic trees and divergence times calibration

To provide a phylogenetic framework for macroevolutionary analyses, we obtained a time-calibrated multilocus species tree using the supermatrix of 38 supplemented by new published data from the Lygodactylus gutturalis group [19] and names updated from [36, 39]. The time-calibrated species tree was generated in BEAST v2.7.4 following the same parameters (calibration points, partition scheme, and best-fitting model) defined in [38]. To check for convergence and stability of the analyses, two separate runs of 10⁸ MCMC generations were performed, and their results compared in Tracer v.1.7.2 [49].

Datasets

Different datasets were put together for the three skull elements (cranium, jaw and inner ear) aiming to optimize the data analysis (Table S2). Dataset 1 (D1) included all the species CT-scanned without including the intraspecific variation. Dataset 2 (D2) include the species from D1 pruned using the molecular phylogenetic tree and the morphological dataset to use only the species included in both datasets. Dataset 3 (D3) include all the intra and inter-specific variation.

Shape analyses

First, we explored skull shape variation within Lygodactylus (cranium, mandible, and inner ear) through different principal component analyses (PCAs) using the Procrustes coordinates in R package geomorph v.4.0.1 [47, 50] for each dataset. To visualize geographical and ecological tendencies of shape variation in morphospace, we plotted the different groups defined above (phylogenetic groups, geographic distribution, macrohabitat, and lifestyle; Tables S1), on the first three PC axes (PC1, PC2 and PC3). Additionally, we explored evolutionary morphological trajectories by projecting morphological (MST tree) phylogenies on the phylomorphospace of D1 and molecular (BI tree) phylogeny on the phylomorphospace of D2. The MST tree was constructed by using Procrustes coordinates to construct a matrix using the unweighted pair group method with arithmetic mean (UPGMA) in PAST4 [51], which was the basis for a complete linkage tree (UPGMA tree) and a range of minimum spanning trees (MST). We projected the branches of the trees onto the morphospace [52], using the R packages phytools v.1.5.1 and geomorph v.4.0.1.

Evolutionary allometry

To determine the effect of allometry on skull shape variation in Lygodactylus, we investigated the effect of size on the evolution of shape of the three skull elements (cranium, mandible, and inner ear). To this aim we tested the allometric component of the skull using ANOVA on D1 and D3 included in the ‘procD.lm’ function of the R package geomorph v.4.0.1 [47, 50], Procrustes-distance-based phylogenetic regressions [53] on D2 included in the ‘procD.pgls’ function of the R package geomorph v.4.0.1 [47, 50], using the Type II sum of squares and 10000 iterations. We considered the coefficient of determination (R²) as a measure of the strength of a relationship and used it to compare models. These models were performed based on multivariate phylogenetic regressions between skull shape (Procrustes shape data) and size (log centroid size; the log of the sum of the squared distances between the landmarks in real scale; [54]), including intra (D3) and interspecific (D1) data, accounting for the phylogenetic relationships among species [55] in MorphoJ v.1.7 [48].

Ecomorphology, geographic diversity and lifestyle

We tested the effect of the different selected evolutionary traits on the skull elements of Lygodactylus. To this aim, we also implemented multivariate analyses (ANOVAs) on D1 and D3, and phylogenetic regressions (PGLS) [53, [56]] on D2, to evaluate the effect of geographic distribution and ecological traits (macrohabitat and lifestyle) on the size (centroid size and snout-vent length ~ SVL [measured from tip of snout to anterior cloaca opening]) and geometries (shape) of the different skull elements (cranium, jaw, and inner ear). Statistical explanation for the evolution of shape of the three elements (cranium, mandible, and inner ear) were tested as using Procrustes-distance-based phylogenetic regressions [53] previously mentioned above. First, as an exploration of the different allometric components affecting the skull, we tested the relationships between size traits and skull shape (i.e., shape ~ size) and then, the additive effects of various size traits (i.e., shape ~ size 1 + size 2). In this way, we were able, for example, to separate the effect of the increase in cranium size to correctly assess the independent effect of labyrinth size. Secondly, we investigated individually the relationship between shapes and biogeography/ecology traits (e.g., shape ~ macrohabitat). Then, we implement combined models to explore interaction between different traits and non-independence among explanatory variables (e.g., shape ~ lifestyle + geographic distribution, and shape ~ lifestyle * geographic distribution), making it possible to determine whether there are differences in form related to the different evolutionary traits. Finally, as an exploration of the relationship between allometry and geographic distribution, macrohabitat or lifestyle, we implemented combined models evaluating the independent effects of combined traits, both allometric and ecological (i.e., shape ~ size 1 + macrohabitat), thus testing whether the influence on shape of the ecological trait is independent of that on size, trying to identify other than size-dependent evolutionary trends. Finally, for those variables that show differences in mean shape a post hoc pairwise comparisons were used to test how different traits are from one another and explore their influences on statistical differences, using ‘pairwise’ function [57] ‘procD.lm’ function of the R package geomorph v.4.0.1 [47, 50].

Disparity and phylogenetic signal

We estimated morphological disparity of the three sets of landmarks (cranium, mandible, and inner ear) and calculated statistical differences using the Procrustes variance by group means and 1000 iterations implemented in the ‘morphol.disparity’ function of the R package geomorph v.4.0.1 [47], dividing the sample into the following groups: (1) phylogenetic groups, (2) geographic distribution, (3) macrohabitat, and (4) lifestyle. We also calculated the disparity-through-time [58] for the three elements, applying the average squared Euclidean distance index and 1000 simulations, using the ‘dtt’ function of the R package geiger v.2.0.7 [59]. We calculated the phylogenetic signal (i.e., the extent to which closely related species show similar values compared to those of more distantly related species) for the shape data applying the ‘physignal’ function of the R package geomorph v.4.0.1 [47], and for the centroid size of the different elements with the ‘phylosig’ function of the R package phytools [60]. Phylogenetic signal was quantified using Blomberg’s K (or Kmult) [61, 62] with predicted values of one for a Brownian Motion model of evolution (diffusive evolution in which lineages are not strongly attracted to certain trait values), values smaller than one indicating weaker phylogenetic signal, and values close to one or larger than one indicating stronger phylogenetic signal.

Results

Cranium

Overall shape

(Figs. 1 and 2). The amount of shape variation between and among species (D1-D3 and D2 respectively), measured in Procrustes distance (> 0.1), spans a small morphospace. This means a homogeneous dispersion (no primary trend) over morphospace, whereby all taxa overlap, regardless their evolutionary trends, spanning a small Procrustes distance (< 0.08), entailing high similitude. The first two principal components (PC1 and PC2) of cranium shape variation explained 24.2% and 16.4%, respectively (Table S5) of the non-phylogenetic morphospace (Dataset 1). Species with a negative PC1 score possess a broader cranium, including relatively broader mandibular elements and neurocranium. Conversely, species with positive scores for PC1 have a narrower cranium with relatively medio-laterally compressed braincase. PC2, however, is related to variation in shape of the cranium the dorsoventral plane: namely, species with negative scores have a taller cranium, and those with positive scores have more depressed cranium. PC3 (6.6%) is related to the proportional size of the postorbitofrontal bone (Fig. S4).

Allometry

The allometric analyses found a significant correlation between the shape variation of the different skull elements (cranium, mandible and inner ear) and size (centroid size and SVL) at intra and interspecific levels, being the highest correlation the one between cranium and the centroid size (ANCOVA: R² = 0.1194, F_(1,62) = 8.271, p value = < 0.001; Fig. 3). This is also visible on the PCAs analysis, on which the specimens with larger centroid size (Cs) occupy the positive values of the PC2, being therefore, those animals that have shorter (or depressed) cranium (Fig. 4). Both centroid size and SVL are highly correlated, exhibiting similar effects in all the analyses (Table S8). In general, cranial shape in Lygodactylus follows what is commonly described in tetrapods as craniofacial evolutionary allometry (CREA; [63]), in which larger species have a larger snouts, dorsoventrally flattened crania, proportionally smaller orbits, and more flattened mandibles [64, 65; Fig. 3]. Conversely, smaller species have a more rounded dorsoventral profile of the cranium, with shorter snout, and proportionally larger eyes (Fig. 3), giving them a juvenile appearance.

Phenotypic trait covariation

The phylogenetic groups occupy different parts of the morphospace, but with high overlap between them (Fig. 4A). Some groups occupy a large morphospace area (like the L. picturatus-group, the L. pictus-group, and the L. capensis-group) whereas others occupy a more restricted area (like the L. bonsi-group and the L. somalicus-group). It is noteworthy that the morphospace, in turn, does not reflect phylogenetic signal, which is visible when mapping the molecular phylogeny (BI tree; Fig. 4E) and when using the K-statistic generalization on the three different skull elements (cranium: K = 0.089, p value = 0.501; mandible: K = 0.124, p value = 0.098; inner ear: K = 0.174, p value = 0.017). This result, revealing high morphological conservatism (homoplasy) is also confirmed by the MST and the UPGMA, which did not recover the same topology as the BI tree when examined for morphological similarity between cranial shapes (Fig. 4A–E).

In addition, the geographic distribution and macrohabitat show similar overall dispersion on the morphospace, with high levels of overlap (Fig. 4B–C). However, we found some segregation of the different lifestyles in the morphospace (Fig. 4D–E). Thus, rupicolous species are only distributed along the negative scores of PC1, with the exception of two species, L. graniticolus and L. rarus. On the other hand, it is noteworthy that the MST and the UPGMA recover a subgrouping within rupicolous species that has no-phylogenetic concordance (Fig. 5). Among the ecological variables, the singular models show that only lifestyle has a significant association with shape (Table 1). This relationship also shows a high fit (R² = 0.725, Z = 2.227, p < 0.001). When evaluating the independent influence of this variable in the combined models along with size, it is observed that the influence of lifestyle on shape is highly codependent on the allometric component, as the explanatory power of both decreases drastically (Table 1, S9 & S11). In contrast, the macrohabitat does not show a significant correlation with shape in simple models, but when considered independently in the combined model with the allometric influence, its explanatory power increases by 50% and gains a significant influence, although relatively low (R² = 0.0325, p < 0.001). This shows how the allometric component might obscure the adaptive component associated with macrohabitat, with the morphological transformations associated with it being completely independent of size. Finally, the geographical distribution does not show a significant influence on shape differences, both in singular models and independently in combined models (Table 1, S9 & S11). The results of the last two variables (i.e., macrohabitat and geographical distribution), confirm the interpretations derived from the results of the overlap in morphospace, as they have very little impact on the main shape differences identified for the Lygodactylus diversity. Finally, the statistical differences between means by groups confirm these results, finding no significant relevant differences between groups associated with different ecological variables (Table 2 & S11), highlighting the overlap between means, and therefore the morphological conservatism among the different groupings. Significant differences are only found between the mixed and arboreal lifestyle, the transitional-zone and rainforest macrohabitat (Table 2), but these differences may be due to only one species belonging to former variables in both comparisons.

Model	Variable	SS	MS	R2	F	Z	P value
Cranium
Single Models
shape ~ cranium size (cs)	Log(cs)	0.0365	0.0365	0.5869	89.5220	2.4633	< 0.001
shape ~ SVL	Log(SVL)	0.0250	0.0249	0.4013	42.2300	2.5663	< 0.001
shape ~ geographic distribution	distribution	0.0007	0.0002	0.0111	0.2284	-1.1489	0.8511
shape ~ macrohabitat	macrohabitat	0.0015	0.0004	0.0236	0.4913	-0.2775	0.6430
shape ~ lifestyle	lifestyle	0.0451	0.0226	0.7255	81.9190	2.2267	< 0.001
Combined Models
shape ~ ecology by size (lifestyle + cranium centroid size)	lifestyle	0.0104	0.0052	0.1671	20.7182	4.3416	< 0.001
	Log(cs)	0.0018	0.0018	0.0286	7.0837	4.6356	< 0.001
shape ~ ecology by size (macrohabitat + cranium centroid size)	macrohabitat	0.0020	0.0007	0.0325	1.7075	1.4131	0.0818
	Log(cs)	0.0371	0.0371	0.5958	93.9416	2.4698	< 0.001
shape ~ ecology by size (distribution + cranium centroid size)	distribution	0.0011	0.0004	0.0186	0.9417	0.1107	0.4564
	Log(cs)	0.0370	0.0370	0.5944	90.4096	2.4819	< 0.001
shape ~ geographic distribution + lifestyle	distribution	0.0007	0.0002	0.0119	0.8910	-0.2909	0.6100
	lifestyle	0.0452	0.0226	0.7263	81.5760	2.2289	< 0.001
shape ~ geographic distribution*lifestyle	distribution: lifestyle	0.0013	0.0013	0.0214	5.1492	3.8174	< 0.001
shape ~ geographic distribution + macrohabitat	distribution	0.0007	0.0002	0.0110	0.2200	-1.1855	0.8697
	macrohabitat	0.0015	0.0005	0.0235	0.4700	-0.3378	0.6497
shape ~ geographic distribution*macrohabitat	distribution: macrohabitat	0.0008	0.0003	0.01255	0.2414	-1.3329	0.9028
shape ~ macrohabitat + lifestyle	macrohabitat	0.0026	0.0009	0.0414	3.4962	4.3666	< 0.001
	lifestyle	0.0462	0.0231	0.7433	94.0721	2.2226	< 0.001
shape ~ macrohabitat*lifestyle	macrohabitat: lifestyle	0.0019	0.0010	0.0308	4.3467	4.5407	< 0.001
shape ~ ecolgeogr.model (distribution + macrohabitat + lifestyle)	distribution	0.0007	0.0002	0.0115	0.9704	0.0100	0.4937
	macrohabitat	0.0025	0.0008	0.0411	3.4590	4.2097	< 0.001
	lifestyle	0.0463	0.0231	0.7438	93.9979	2.2056	< 0.001

1. Significant results are bolded. To see the complete results of the PGLS analysis interaction see Table S8–S9 & S11. SVL = snout–vent length

Geographic Distribution					Macrohabitat					Lifestyle
Cranium
	1	2	3	4		5	6	7	8		9	10	11
1. Africa	0	0.0028	0.0002	0.0001	5. Dry forest	0	0.0002	0.0005	0.0036	9. Arboreal	0	0.0037	0.0002
2. America	0.067	0	0.0030	0.0029	6. Montane	0.68	0	0.0002	0.0038	10. Mix	0.036*	0	0.0036
3. Oceanic Island	0.772	0.085	0	0.0001	7. Rain forest	0.428	0.767	0	0.0040	11. Rupicolous	0.708	0.057	0
4. Madagascar	0.784	0.060	0.929	0	8. Transition Zone	0.082	0.071	0.039*	0
Mandible
	1	2	3	4		5	6	7	8		9	10	11
1. Africa	0	0.0016529	0.0003	0.0010	5. Dry forest	0	0.0013	0.0001	0.0043	9. Arboreal	0	0.0043	0.0010
2. America	0.149	0	0.0020	0.0026	6. Montane	0.031*	0	0.0012	0.0031	10. Mix	0.038*	0	0.0033
3. Oceanic Island	0.792	0.201	0	0.0007	7. Rain forest	0.876	0.064	0	0.0042	11. Rupicolous	0.076	0.072	0
4. Madagascar	0.093	0.066	0.474	0	8. Transition Zone	0.039*	0.093	0.042*	0
Inner Ear
	1	2	3	4		5	6	7	8		9	10	11
1. Africa	0	0.002	0.0026	0.0011	5. Dry forest	0	0.0014	0.0005	0.0090	9. Arboreal	0	0.0093	0.0011
2. America	0.489	0	0.0000	0.0013	6. Montane	0.311	0	0.0019	0.0076	10. Mix	0.066	0	0.0082
3. Oceanic Island	0.262	0.098	0	0.0014	7. Rain forest	0.762	0.223	0	0.0095	11. Rupicolous	0.428	0.084	0
4. Madagascar	0.412	0.071	0.054	0	8. Transition Zone	0.080	0.119	0.076	0

1. Displayed values represent significance levels (p-values) obtained via permutation (RRPP). Different traits are geographic distribution, macrohabitat and lifestyles. Distances above the diagonal represent differences in morphological disparity among the evolutionary traits. Asterisks (*) depict significant values (p ≤ 0.05)

Mandible

Overall mandible shape

(Fig. 2). The first three principal components of the PCA explain 40.8% of the variance, and PC1 and PC2 explain 20.7% and 11.1% of the variation in cranium shape, respectively (Table S6). PC1 is related to the shape of the coronoid bone, whereby positive scores possess a lower coronoid bone, and the negative scores present a more prominent eminence of the coronoid bone. PC2, however, is related to the length of the mandible, with positive scores showing larger mandible with proportionally larger compound bone + surangular (see [31], for comments on the compound bone + surangular structure), and negative scores relate to a shorter mandible with proportionally shorter compound bone + surangular (Fig. S5).

Allometry and phenotypic trait covariation

Results of the PGLS analyses are summarized in Table S8, S9 and S11. We found similar results for the mandible as for the cranium, with a significant effect of size (Mandible centroid size [Mcs] and SVL) on the mandible shape (PGLS (Mcs): R² = 0.592, Z = 2.477, p value = < 0.001). Looking at the geographic distribution, and ecological traits (macrohabitat and lifestyle) on the mandible, it seems clear that all the three groupings are similarly distributed in the morphospace with high overlap between them (Fig. 6A). These results agreed with the results of the cranium as reflected in the PGLS analyses (Table S9).

Inner ear

Overall inner ear shape

(Fig. 7). The first three principal components performed with the full dataset explain more than in other cases (53.3%), with PC1 and PC2 explaining 28.7% and 15.0% of the variation, respectively (Table S7). The inner ear in all species exhibits the common structure in reptiles, with a proportionally larger vestibule, a moderately large endosseous cochlear duct and three semicircular canals, which separate from the ampullae, and curve in different directions, joining in the common crus ([31]; Fig. 7]). However, small deviations from the average shape (0,0) are visible along the major principal components (see Fig. 6B & S6), mainly in the anterior and horizontal semicircular canals. Inner ears with positive scores on PC1 are more rounded overall, and the anterior semicircular canal elongates anteriorly, and then curves anteroventrally (downwards) including the anterior ampullary recess. Negative scores are associated with more elongated inner ears, and the anterior semicircular canal curves anterodorsally (upwards), extending to the common crus (Fig. 6B).

Allometry and phenotypic trait covariation

Results of the PGLS analyses are summarized in Table S8, S9 and S11. We found similar results for the inner ear as for the cranium and the mandible, with significative effect of size (Inner ear centroid size [Ecs] and SVL) on the mandible (PGLS (Ecs): R² = 0.56, Z = 2.471, p value = < 0.001). The examination of the different traits in the inner ear morphospace shows that all the three traits are widely distributed, entailing nearly no distinction between their inner ear shapes (Fig. 6). Of note, the PGLS results again show a significant relationship between lifestyle and shape, in agreement with the cranium and mandible (Table S9). Therefore, the three skull elements presented similar results on the different evolutionary analyses. Although, the cranium represents the skull element with higher difference on the lifestyle, either in the morphospace and the multivariate analyses when compared to the other two skull elements.

Disparity and tempo of shape evolution

Finally, the different skull elements (cranium, mandible, and inner ear), generally show an average subclade relative disparity over time (or ‘disparity-through-time, dtt) that resembles the distribution expected for a non-adaptive diversification, describing a pattern well-above (i.e., higher values of average subclade disparity, and thus morphological variations are almost equivalent between and among clades) in most of the phylogenetic nodes, than expected by Brownian Motion (neutral evolution). Morphological disparity indexes (MDI), which represent the overall difference between the observed variation among-clades disparity and the null distribution, show positive values (ranging between 0.089 for the cranium and 0.174 for the inner ear). This means that the disparity of the different subclades is evenly distributed across the morphospace (i.e., the cranial configurations diversify on already existing morphologies, highly overlapping in morphospace) (Fig. 5).

Discussion

Lygodactylus geckos represent one of the most unique cranial architectures within gekkotans, presumably related to their miniaturization [31]. These geckos, nevertheless, present variable ecological traits within the group [19, 38]. In this study, we have investigated the evolutionary trends in patterns of skull and inner ear shape variation at different taxonomic levels (i.e., generic, group and species) and how these trends may be correlated by allometric, phylogenetic, geographical, and ecological factors. The obtained results stress high levels of skull geometric similarity (homoplasy) across species, thus corroborating a strong morphological conservatism on Lygodactylus skull morphology. Notwithstanding this, the amount of shape variation between and among species, including the most basal group (i.e., the L. madagascariensis group) is small (measured in Procrustes distance). Accordingly, the skull shapes of Lygodactylus have undergone an evolutionary pattern of homogeneous variation. Given the consistency across the entire sample, where both the more basal and more derived clades maintain the same skull geometry, this could be interpreted as symplesiomorphy. However, the consistency of skull shape in Lygodactylus needs to be further tested in relation to the known evolutionary pattern of bone loss and/or non-ossification [31] in other related taxa (e.g., Phelsuma spp.), as well as its relationship with allometry and miniaturization.

Congruent with the assumption that size is an important factor constraining shape variation in miniatured geckos [66], we also found that most of the craniofacial shape variation of the Lygodactylus skull is significantly related to size (i.e., evolutionary allometry; [67]), in a similar way to most vertebrates, with size influencing the enlargement of the face [3, 63]. Moreover, we found that different ‘ecological morphs’ are equally distributed within morphospace, comparable to the variation seen at the intraspecific level, without any evidence suggesting ecological or functional diversification (namely, adaptability). Therefore, apart from the allometric component of the skull, the evolutionary hypothesis that explains the observed diversity in skull shape in Lygodactylus, contrary to the adaptive hypothesis, is that this variation is intrinsic to the species, likely explainable through developmental plasticity (68,69,70,71). Consequently, changes along ontogeny determine skull shape variation in Lygodactylus. In fact, the analysis of inner ear variation also supports the discrepancy between phylogenetic divergence, geographic distribution, and ecological trends, contradicting the expected adaptive signal observed in other reptiles at the genus level (e.g., Anolis [32]; Podarcis [72]; or turtles [34]). Perhaps the lack of an adaptive signal is a consequence of miniaturization in Lygodactylus compared to the larger size variation seen in other squamate groups, which may diminish such a functional signal due to a more constrained allometric range [73].

Intriguingly, the multivariate analyses detected subtle skull differences attributable to a dichotomy between different lifestyles (arboreal and rupicolous). It must be noted that a high proportion of rupicolous species are larger species and exhibit a more elongated and flattened cranium shape, which might be related with rupicolous habits (i.e., cracks and crevice-dwelling adaptation; [74]). However, other rupicolous species exhibit convergence with arboreal species, displaying a more rounded and shortened cranial shape. Additionally, pairwise results indicate that rupicolous species are equally different and/or similar to each other as they are to other arboreal groups. Specifically, certain phylogenetic groups that include species with different lifestyles (rupicolous and arboreal) exhibit no significant differences between them (e.g., the L. verticillatus-group or the L. bonsi-group), contradicting the adaptive hypothesis. Therefore, the rupicolous phenotype could be interpretable as a result of an incomplete convergent evolution, assuming these ecomorphotypes were undergoing ecological diversification [7, 8].

Interestingly, the intraspecific variation in L. kibera is nearly equivalent to that characterizing the entire genus. Therefore, the differences observed between different lifestyles may be biased by an incomplete sampling of intraspecific variation across Lygodactylus. However, this issue requires further study with more extensive intraspecific sampling of other species, emphasizing the need for increased sampling efforts in both rupicolous and arboreal morphs.

Conclusions

We provide a fine-scale study of cranial diversity in the lizard genus Lygodactylus, including an almost taxon-complete sample of an evolutionary radiation. This has allowed demonstrating that skull elements (cranium, mandible, and inner ear) show a high degree of morphological conservatism in Lygodactylus. Accordingly, a low phylogenetic signal contrasts with a strong allometric component in shaping cranial diversity between and among species, potentially associated with miniaturization. Craniofacial allometry at a generic level might have partially guided evolution along paths of least resistance [4, 75], which is consistent with what would be expected in a non-adaptative radiation. In addition, we only found a weak signal of lifestyle (arboreal vs. rupicolous) contributing to cranial diversity. However, this signal may be spurious, since not only the differences between and among different lifestyles are comparable across various macrohabitats and geographic distributions, but so are the observed cranial diversity differences in L. kibera compared to the entire genus. Accordingly, our results suggest that miniaturization in Lygodactylus has resulted in a unique cranial morphology, homoplastically maintained across various axes of diversity. Consequently, our results support the non-adaptive hypothesis previously proposed to explain this radiation of dwarf geckos.

Data availability

The specimens in this study are in the collections specified in the Material and Methods section. The CT data generated and analyzed are available at Morphosource (https:// www.morphosource.org/; Project ID: 00000C975) and Table S2 (Supplementary material). All other morphometric data are available in the Supporting Information, (Supplementary material Table S1).

Abbreviations

CT:

Computed tomography

CREA:

Cranio-Facial Evolutionary Allometry

SVL:

Snout-venter length

PGLS:

Phylogenetic generalized least squares

cs:

Centroid size

Mcs:

Mandible centroid size

Ecs:

Inner ear centroid size