Introduction

Locomotion is a key function for animal survival and crucial for most aspects of a motile animal’s life. Consequently, the locomotor apparatus is often optimised for performance in the structural microhabitat occupied and for the substrates animals encounter in nature (Garland & Losos, 1994; Irschick & Garland, 2001; Pillai et al., 2020), with the proviso that the degree of specialisation towards proportions of the available structural habitat space may vary along the generalist to specialist continuum. Thus, over evolutionary time, changes in the morphology (and performance) of the locomotor apparatus can serve as adaptations to locomotion in different microhabitats (Losos, 2011). Such processes, which are conceptualised as ecomorphology (Losos, 1994), can lead to adaptive radiations when different populations repeatedly and independently differentiate over the available structural microhabitat niche axis by optimising their morphology for distinct regions of the available morphospace (Stroud & Losos, 2016). In limbed vertebrates the appendages serve as lever arms for locomotion, and thus relative limb and digit dimensions are often associated with the locomotor performance and behaviour of vertebrates. Squamate lizards have been studied often as model clades in this context (e.g., Arnold, 1998; Losos, 1992; Melville et al., 2006).

Based on biomechanical considerations, two sets of predictions have repeatedly been suggested for limb proportions in lizards. On one hand, relatively long hindlimbs (and particularly longer femoral and pedal lengths) should allow for greater stride length and thus higher running speed, while the forelimb should be shorter relative to the hindlimb so as to not interfere with the hindlimbs at higher running speeds (Irschick & Jayne, 1999; Russell & Bels, 2001). This should be beneficial for cursorial species occupying open terrestrial habitats and employing flight behaviour as their main predator avoidance strategy (Arnold, 1998; Higham & Russell, 2010). Scansorial species, in contrast, are expected to have relatively shorter limbs to keep the centre of mass (CoM) closer to the substratum, and the fore- and hindlimbs should be more equal in length as they are predicted to be equally important for climbing (Zaaf & Van Damme, 2001). Alternatively, terrestrial species occupying closed habitats have been predicted to benefit from relatively shorter limbs as they may encounter more obstacles when navigating through dense low ground vegetation (Arnold, 1998), whereas scansorial species are predicted to benefit from longer limbs promoting greater limb spans and higher static stability (Foster et al., 2018).

Thus, there are two seemingly contradictory predictions for relative limb dimensions in lizards. For terrestrial species the contradiction is easily resolved by consideration of a second ecological factor – habitat openness. In this regard, the predicted pattern of long limbs in open habitats and short limbs in closed habitats has been confirmed repeatedly (Arnold, 1998; Herrel et al., 2002; Melville & Swain, 2000). For scansorial species, the picture is less clear. Scansorial microhabitats can vary in various parameters, including canopy cover, perch diameter (of arboreal habitats), or locomotor surface (e.g., arboreal or saxicolous). Within the arboreal realm, perch diameter has so far been a good predictor of relative limb dimensions. Among Anolis lizards, distinct arboreal ecomorphs have evolved repeatedly (Langerhans et al., 2006; Losos, 1990), and relative hindlimb length increases with diameter of the perches occupied in nature (Losos, 1994). In Tropidurus and Draco lizards, canopy species using narrow branches have relatively shorter hindlimbs (Kohlsdorf et al., 2001; Ord & Klomp, 2014). In Australian geckos, however, the correlation between limb proportions and perch diameter is reversed, with longer limbed species using narrower perches (Hagey et al., 2017). Thus, the correlation between perch diameter and limb proportions may be clade specific. Saxicolous species often have relatively longer limbs compared to their closely related terrestrial counterparts (Foster et al., 2018; Goodman et al., 2008; Norris et al., 2021), and sometimes also compared to arboreal species (Hagey et al., 2017). In terms of surface curvature, rock surfaces should be more comparable to broad tree trunks with wide perch diameter, but other factors such as perch angle (Riedel et al., 2020) or physical properties of the perch surface (Pillai et al., 2020; Russell & Johnson, 2014) can vary considerably between and within these two scansorial microhabitats.

Cyrtodactylus Gray, 1827 is the most speciose gecko genus (354 species) and the third most speciose vertebrate genus (Grismer et al., 2021c; Uetz et al., 2022). Most of these species have recently been assigned to one of 10 different microhabitat use categories (Fig. 1), including specialists found only on certain types of rocks (granite or karst) or occupying distinct strata within the arboreal niche (species restricted either to large tree trunks or to the narrower branches of the canopy) (Grismer et al., 2020, 2021b). With multiple independent transitions between these microhabitat categories that have occurred during its evolution, this genus is ideally suited for ecomorphological studies (Grismer et al., 2020). A phylogeny of the whole genus is now available (Grismer et al., 2021c), providing the foundation for such an analysis. We leverage this microhabitat and phylogenetic information to examine the evolution of limb proportions in relation to microhabitat use in Cyrtodactylus. Previous analyses of subclades within the genus have revealed that karst and cave specialists have relatively longer limbs than their closely related generalist congeners (Grismer & Grismer, 2017; Kaatz et al., 2021). Neither terrestrial nor arboreal Cyrtodactylus species have yet been included in such analyses, although the former are notably small and squat (Grismer et al., 2020) and seem to be more prevalent in closed habitats (Agarwal & Karanth, 2015), whereas the latter are among the largest species within the genus (Nielsen & Oliver, 2017; Oliver et al., 2014).

Fig. 1
figure 1

Representatives of the 10 ecotypes of the genus Cyrtodactylus a–k modified from Fig. 4 of Grismer et al., (2021b), and a distribution map of the genus Cyrtodactylus modified from Fig. 3 of Grismer et al., (2021b). (a) Cyrtodactylus yakhuna, Sri Lanka; (b) C. seribuatensis, Malaysia; (c) C. oldhami, Thailand; (d) C. payacola, Malaysia; (e) C. kulenensis, Cambodia; (f) C. trilatofasciatus, Malaysia; (g) C. sanpelensis, Myanmar; (h) C. honteensis, Vietnam; (i) C. serratus, Papua New Guinea; (j) C. elok, Malaysia. Photographs by (a) Suranjan Karunaratha, (b, c, d, f, g, h, j) L. Lee Grismer, (e) Peter Geissler, (g) Evan S. H. Quah, (i) Steve J. Richards

Full size image

We measured limb and body dimensions of Cyrtodactylus geckos covering 26 of the 31 subclades of the genus and all ecotypes. We assessed if, and how, relative limb dimensions are associated with microhabitat use. We predicted that terrestrial species will have relatively shorter overall limb proportions due to their prevalence in closed habitats, whereas saxicoline species are predicted to have relatively longer limbs in line with previous results (Kaatz et al., 2021). Within the different saxicoline ecotypes, we explored eventual differences and predicted cave dwelling species to be at the extreme end of greater relative limb length because we can expect them to experience, on average, steeper perch angles (Riedel et al., 2020). For arboreal species we predicted different patterns for trunk and canopy dwelling species due to differences in perch diameter. Because of contradictory results in the literature (Hagey et al., 2017; Irschick et al., 1997), however, we did not formulate an a priori prediction favouring one of the alternatives. For generalist species we predicted that they would exhibit limb proportions between those of the various specialists, overlapping with most of them.

Material and Methods

Study System

The geographic distribution of Cyrtodactylus ranges from India and Nepal across South-East Asia towards the Solomon Islands and northern Queensland (Australia), with a peak diversity in Indochina and Sundaland (Grismer et al., 2022a, 2022b) (Fig. 1k). Currently 243 species, including undescribed species, have been assigned ecotype status (Grismer et al., 2020, 2021b): most of these species have been assigned to the generalist microhabitat category (134 species), these being found in all available microhabitats, including the ground. Many species are closely associated with different rock formations and can thus more generally be defined as being saxicolous. Among these are karst (86 species), granite (32 sp.), sandstone (1 sp.), and cave (7 sp.) associated species. The latter are defined as exclusively occupying cave-like environments formed by large granite boulders. Other species are associated with various strata of the vegetation, comprising trunk (38 sp.), ‘arboreal’ (13 sp.), and swamp (5 sp.) species. As these three all fit the general definition of arboreal species, we rename the ‘arboreal’ category of Grismer et al., (2020, 2021b) as the crown ecotype to avoid confusion with the aforementioned categories because they are normally restricted to smaller diameter perches in the higher strata of trees and shrubs. The trunk ecotype is more common on trunks and large branches of large trees, whereas swamp species are restricted to smaller perches relatively low to the ground in swampy areas. Twenty-nine species are exclusively terrestrial, and one species is categorised as intertidal. Ancestral state reconstruction suggests that the common ancestor of the genus was a generalist species and that multiple independent transitions among these ecotypes have occurred (Grismer et al., 2020). For the current study we sampled 87 species (with 3–10 specimens per species) representing all ecotypes, using ethanol-preserved museum vouchers (Fig. 2a, Supplementary Table S1). Species status of older voucher specimens were checked using taxonomic literature (Agarwal & Karanth, 2015; Davis et al., 2019, 2021, 2023; Grismer et al., 2014; Termprayoon et al., 2023; Welton et al., 2010a, 2010b) and by correspondence about them with experts on the respective clades (L. Welton pers. com.; S. Mecke pers. com.). We included only adult specimens to minimise allometric effects. Furthermore, we did not differentiate between sexes because geckos are generally not sexually dimorphic (Norris et al., 2021). Supplementary Table S2 provides a detailed list of the voucher numbers and sources of all measured specimens. Ecotypes, for which we had only one representative species (cave, intertidal, sandstone, and swamp) are included in the morphospace analysis and ancestral state reconstructions but are excluded from statistical comparisons among ecotypes (see Data Analysis).

Fig. 2
figure 2

Phylogenetic relationships a of the species of Cyrtodactylus included in this study, based upon the phylogeny of Grismer et al., (2022a). Tip labels are colour coded by ecotype, while black lines and numbers indicate species groups. Major shifts in habitat use (ecotype classification) are indicated as coloured circles at the nodes, as reconstructed by Grismer et al., (2021b). b Length measurements of body and limb proportions recorded for the species examined in this study. SVL snout-vent length, AGL axilla-groin length, BrL brachium length, AbL antebrachium length, MaL manus length, ThL thigh length, CrL crus length, PeL pes length

Full size image

Measurements

The following parameters were measured externally in mm using a digital calliper: Brachium length (BrL), measured as the distance between shoulder and elbow on the dorsal surface of the forelimb; Antebrachium length (AbL), the distance from the tip of the elbow to the wrist; Manus length (MaL), the distance from the wrist to the tip of the 4th finger of the hand; Thigh length (ThL), the distance from the insertion of the hind limb into the body to the knee, measured along the long axis of the thigh; Crus length (CrL), the distance from the tip of the knee to the ankle; Pes length (PeL), the distance from the ankle to the tip of the 4th toe; Axilla-groin length (AGL), measured from the posterior margin of the insertion of forelimb into the trunk to the anterior margin of the insertion of the hind limb; Snout-to-vent length (SVL), measured from the tip of the snout to the anterior margin of the vent (Fig. 2b). AGL is the more effective measurement for assessing aspects of locomotor morphology (Bauer et al., 1996; Powell & Russell, 1992), but SVL is the more commonly employed trait for body size in ecomorphological studies. Therefore, we applied AGL as our primary measure of body size but included both traits for comparability. All measurements were repeated thrice to reduce measurement error and mean values were calculated for each specimen. For phylogenetic comparative methods requiring single measurements per species mean plus standard deviation was then calculated per species. Measurements were taken from the left side of the body whenever possible, but in cases of damaged or missing limb segments the right side was employed instead, on the assumption of bilateral symmetry.

Data Analysis

All statistical analyses were conducted in R version 4.2.2 and R Studio version RStudio 2022.12.0 + 353 (R Core Team, 2023). All measurements were log-transformed to improve normality and reduce heteroscedasticity. To control for phylogenetic signal in our data (Felsenstein, 1985), we used the phylogeny of Grismer et al., (2022a), pruned to contain only the species included in our dataset. Two species (Cyrtodactylus aaroni Günther & Rösler, 2003 and C. irianjayaensis Rösler, 2001) included in our study, are not in that phylogeny and have not previously been assigned to an ecotype. Cyrtodactylus aaroni is morphologically most similar to C. mimikanus (Boulenger, 1914) and the description of its microhabitat matches that of the crown ecotype (Günther & Rösler, 2003). We therefore assign C. aaroni to the crown ecotype and place it in the position of C. rex Oliver et al. (2016) retaining the branch length, one of the closest relatives of C. mimikanus for which we have data. Cyrtodactylus irianjayaensis is putatively a trunk species (Rösler pers. com.) and most closely related to C. novaeguineae (Schlegel, 1837) (Rösler et al., 2007). Thus, we assigned it to the trunk category and placed it as the sister species of C. novaeguineae, replacing C. zugi Oliver et al. (2008) in the original phylogeny while again retaining the branch lengths.

To account for variation in body size, we used AGL, as it is more closely related to the locomotor function of the limbs than SVL (Bauer et al., 1996; Powell & Russell, 1992), and applied and compared two different approaches. First we regressed each limb measurement individually against body size, applying phylogenetic generalised least square regression (PGLS) models (Martins & Hansen, 1997; Revell & Collar, 2009) and extracted the residuals for further analysis (Mahler et al., 2010). To implement the PGLS models we used the ‘gls’ function of the nlme package (Pinheiro et al., 2018) and the ‘corPagel’ function of the ape package (Paradis et al., 2004, 2018). We included a maximum likelihood (ML) estimation of the scaling parameter λ, where λ = 0 indicates no phylogenetic signal and λ = 1 indicates a strong phylogenetic signal converging to a Brownian motion model of trait evolution (Freckleton et al., 2002) An exception was the model for CrL, for which the maximum likelihood estimation of λ failed to converge on an optimum. Here, we compared three models with λ values of 0, 0.5, and 1 using the Akaike information criterion (Burnham & Anderson, 2004). To account for within species variation, we added standard errors as fixed weights. As a second approach, all measurements were converted into size corrected variables by applying an allometric growth model as implemented in the R-package ‘GroupStruct’ (Chan & Grismer, 2022). This model is based on the allometric formula of Thorpe (1976), but adjusted to differentiate between and within population differences in body size, with species being our grouping factor. The MANOVA, PCA, phylomorphospace analyses and the Ancestral state reconstructions were then repeated for both sets of values, and we refer to these as RES (residual) and GS (GroupStruct) respectively.

Initially, to test for overall morphological differences among ecotypes we performed a phylogenetic MANOVA, which tests for deviations from Brownian Motion, using the ‘anov.phylo’ function from the Geiger package run with 10,000 simulations (Harmon et al., 2008). To examine which ecotypes differed and how, we first conducted a principal component analysis (PCA) including all measurements (scaled and centred around 0) using the ‘prcomp’ function of the stats package (v4.2.2, R Core Team, 2023) and plotted the morphospace for axes that jointly explain greater than 90% of the variance. Next, we repeated the PCA analysis with a reduced dataset excluding ecotypes of which we had fewer than 3 species to conduct a phylogenetic ANOVA for the principal component axes that explain at least 5% of the variance, to test which ecotypes were separated by which combination of morphological parameters. The scores of the respective principal component (PC) axes were used as response variables and ecotype was used as the explanatory variable. This test was implemented with the ‘phylANOVA function (phytools package, Revell, 2012) computed on 10,000 simulations each, and followed by a post hoc test accounting for multiple p-values using the Benjamini & Hochberg method (Benjamini & Hochberg, 1995). To explore patterns of evolutionary change in the morphospace we additionally plotted the phylomorphospace for the same axes using the ‘phylomorphospace’ function of the phytools package (Revell, 2012). The phylomorphospace plot was overlain with an ancestral state reconstruction of the ecotypes implemented with the ace function of the ape package (Paradis et al., 2004, 2018) and an equal rates model (ER) was implemented following Grismer et al. (2021a, 2021b).

Because the different segments of the fore- and hindlimbs loaded almost equally in the PCA analyses (see results), we then recorded forelimb length (FLL) and hindlimb length (HLL) as the sum of the individual segments of the extremities and size corrected these as described above. Next, we examined the evolution of body and limb proportions by performing an ancestral state reconstruction of SVL, AGL, and size corrected measures of FLL and HLL. For these reconstructions we used the ‘fastanc’ algorithm with the ‘contMap’ function in phytools (Revell, 2012). Finally, we constructed four phylogenetic generalised least square regression (PGLS) models using the reduced dataset as described above. Here the model for HLL failed to converge on an optimum for λ; thus, we compared models as described before (Burnham & Anderson, 2004). The log-transformed original measurements were used as response variables and Ecotypes as explanatory variables, with log transformed AGL as additional response variable for FLL and HLL. The models were analysed with a Type II ANOVA (Langsrud, 2003) and post hoc comparison implemented with the emmeans package (Russell, 2019).

Results

Measurements and MANOVA

Measurements for all specimens examined are provided in Supplementary Table S2. Overall, trunk, granite, and karst dwelling species have relatively long limb proportions and the former two also have elongated bodies (indicating larger size), whereas terrestrial species have relatively low values for all measurements (Fig. 3). The phylogenetic MANOVA revealed significant morphological differences between ecotypes in both datasets (PGLS: F1,9 = 3.226, p < 0.001; GS: F1,9 = 3.281, p < 0.001).

Fig. 3
figure 3

Box plots of body length (SVL and AGL) and limb measurements (in mm), summarized by ecotype, of the species of Cyrtodactylus examined in this study. n = number of species sampled per category. For number of specimens per species see Supplementary Table 1

Full size image

PCA and Phylomorphospace

The PCA analysis differed between both datasets. For the RES data, the first four principal component axes explained 93.18% of the overall morphological variance, while for the GS data 91.74% was explained by the first two axes, but with PC2 explaining less than 5% (Table 1). Even so, we plotted the first two GS PC axes for better visualisation. Tarit loadings varied mostly with regard to AGL. On the RES PCA all limb segments had strongly negative loadings on PC1, while on GS PC1 all measurements, including AGL, loaded uniformly negatively (Table 1). RES PC2 distinguished between AGL plus PeL and ThL on the positive axis and Abl and BrL on the negative axis. GS PC2 showed the same pattern but with CrL additionally contributing to the positive axis. On RES PC3 AGL and MaL had strong negative loadings, while ThL and CrL loaded strongly positively. On RES PC4 AGL and BrL loaded strongly negatively, while PeL and MaL loaded strongly positively (Table 1).

Table 1 Loadings for the morphological traits on the first four PC axes of RES PCA and the first two axes of GS PCA covering 90% of the total variation in both analyses
Full size table

Ecotypes varied in axilla-groin (interlimb) length, and limb segment proportions. In the RES morphospace (Fig. 4a) terrestrial species clustered in the negative regions of PC1 and 2, opposing the granite, karst, and trunk species. Both saxicoline ecotypes showed great overlap although karst species occupied a broader and more central area. Trunk species clustered in the positive area of PC2 but relatively centrally on PC1. Crown species exhibit the narrowest range in the morphospace and lie in the negative area of both axes, and generalists occupy the greatest area, overlapping all other groups except the single swamp and cave species. The GS morphospace differed in that among-ecotype differences are restricted to PC1, with karst, granite and trunk ecotypes more in the positive area and terrestrial and crown species in the negative (Fig. 4b). The RES morphospace for PC3 and 4 segregated trunk species from karst and terrestrial species mostly on PC3, with granite, crown, and particularly generalist species overlapping with one another and the other ecotypes (Supplementary Fig. S1c).

Fig. 4
figure 4

Ordination plots (a, b) and phylomorphospace (c, d) of the first two axes of a principal component analysis (PCA) performed employing the RES data (a, c) and the GS data (b, d) for seven continuously varying morphological dimensions for 87 species of Cyrtodactylus. Morphological measurements: AGL axilla-groin length, BrL brachium length, AbL antebrachium length, MaL manus length, ThL thigh length, CrL crus length, PeL pes length

Full size image

Species of different ecotypes do show significant differences in their scores on RES PC1 (F1,5 = 5.180, p = 0.004), RES PC2 (F1,5 = 6.862, p = 0.008), RES PC3 (F1,5 = 5.449, p = 0.028), and GS PC1 (F1,5 = 11.499, p < 0.001) after accounting for phylogeny (Table 2). For RES PC1, karst (1.303 ± 2.018) species had lower values than crown (− 2.567 ± 1.552) species. The differences between karst and terrestrial species approached significance, with the former showing somewhat higher values (Supplementary Table S3a). For RES PC2, differences among all ecotypes were insignificant, although some differences approached significance with karst showing somewhat higher values than terrestrial species, as did granite and trunk compared to terrestrial and generalist species (Supplementary Table S3b). Trunk (− 0.972 ± 1.207) species had significantly lower values than karst (− 0.596 ± 0.726) ones on RES PC3, with the differences towards terrestrial species again approaching significance (Supplementary Table S3a). On GS PC1, terrestrial species (− 3.68 ± 1.64) have significantly lower values than karst (0.65 ± 1.47), granite (1.95 ± 1.85), and trunk (2.45 ± 2.82) species, but generalist (− 0.87 ± 2.06) and crown (− 1.46 ± 1.57) species have significantly lower scores than granite and trunk species (Table 2). Differences between terrestrial and generalist species approach significance, with the former exhibiting somewhat lower values than the latter (Supplementary Table S3).

Table 2 Phylogenetic ANOVA and post hoc comparison for PC axes 1 to 4 of RES PCA and 1 of GS PCA
Full size table

Species from similar ecotypes cluster in morphospace relatively independently of their phylogenetic placement (Fig. 4c, d, Supplementary Fig. S1b, c). Due to the similar direction and relatively equal loadings of ThL, CrL, and PeL, and BrL, AbL, and MaL, respectively, on all relevant axes except for RES PC3 & 4, we combined these two sets of traits as hindlimb length (HLL) and forelimb length (FLL) for downstream analysis. Both traits were calculated as the sum of the respective segmental measurements, and then log-transformed and size corrected as described previously.

Ancestral State Reconstructions

The ancestral state reconstructions of the relative length of the extremities differs between the RES and the GS data. While the RES data showed differential patterns for relative front and hind limb length, which also differed from the pattern of both body length measures (AGL and SVL) (Fig. 5), the GS data showed a far more uniform pattern among all four traits (Supplementary Fig. S2). This result aligns with the equal loading of all traits on GS PC1 (Table 1), indicating that the GS measures of the extremities might not be completely decoupled from body size (see discussion). We therefore focus here on the results of the RES data.

Fig. 5
figure 5

Ancestral state reconstruction of log-transformed a snout-vent length (SVL), b axilla-groin length (AGL), c RES forelimb length (FLL), and d RES hindlimb length (HLL) for the species of Cyrtodactylus examined in this study. For GS FLL and GS HLL see Supplementary Fig. S1. Phylogeny based on Grismer et al., (2022a). The reconstruction was conducted applying the ‘fastanc’ algorithm in the R package phytools (Revell, 2012). Black lines indicate phylogenetic groups. Coloured circles before the species names indicate ecotype (colour code as in Fig. 2a). Black circles at basal nodes indicate shifts in habitat use at that node analogously to Fig. 2a. 1: generalist ancestor of the genus, 2, 11 & 13: shift from generalist to karst habit; 3, 5, 9, 10: shift from generalist to trunk habit, 4 & 12: shift from generalist to crown habit., 6: shift from generalist to terrestrial habit, 7: shift from generalist to granite habit, 8: shift from granite to karst habit

Full size image

Most ancestrally generalist species are characterized by ancestrally intermediate values for all four traits, although some ancestrally generalist species (e.g., C. papuensis (Brogersma, 1934), C. cattienensis Geissler et al., 2009, C. rubidus (Blyth, 1861)) evolved a shorter AGL (Fig. 5b) and a moderately shorter SVL (Fig. 5a). Secondarily generalist species, such as C. zebraicus Taylor, 1962 (karst ancestor) and C. mimikanus (trunk ancestor) exhibit reduced body size and relative hind limb proportions compared to reconstructed scansorial ancestors, while C. fraenatus (Günther, 1864), evolving from a terrestrial ancestor, exhibits longer hind limbs. A noticeable exception to this overall pattern is the ancestrally generalist philippinicus group, in which some derived generalist species show a moderate increase in all four traits. Terrestrial species have evolved short bodies and trunks (Fig. 5a, b). Limb proportions are reduced in terrestrial species, except for in the triedrus group (Fig. 5c).

Within the two arboreal ecotypes, crown species evolved overall shorter bodies, trunks, and limb proportions (Fig. 5). This trend is more pronounced in relative limb proportions in the brevipalmatus group compared to C. aaroni or C. lateralis (Werner, 1896). Trunk clades evolved overall greater values for body length, trunk length, and moderately elongated hind limbs (Fig. 5). Here, the philippinicus group again deviated from this overall pattern in that it’s trunk species retained medium trait values except for a reduction of AGL in C. ingeri Hidika, 1990.

Saxicoline clades and species evolved longer hindlimbs (Fig. 5d). Within the saxicoline species, the pulchellus group stands out as additionally exhibiting a strong increase in body and trunk length (Fig. 5a-b). This trend is stronger in granite species than in the karst or generalist species. A similar trend is present in the karst-associated linnwayensis group, while the karst-associated chauquangensis group shows longer front limbs (Fig. 5c).

PGLS Models

For HLL, model comparison revealed a λ value of 0 as optimal (Supplementary Table 3). All other traits showed strong phylogenetic signal, with lambda values ranging from 0.617 (SVL) to 0.811 (AGL) (Table 3). The ANOVAs for all models revealed significant differences among ecotypes (Table 3), and a significant positive effect of AGL on both limb measures (FLL: χ21,1 = 235.137, p < 0.001, HLL: χ21,1 = 706.711, p < 0.001).

Table 3 Results of the PGLS model, showing p and support values for Ecotypes
Full size table

For SVL, the post hoc test showed that karst species are relatively shorter than generalist, granite and trunk species, with the latter two also being longer than terrestrial species, while crown species overlap with all other ecotypes (Fig. 6a). The pattern for AGL was similar to that observed for SVL, except that terrestrial species also had lower values than generalists while the latter had lower values than granite and trunk ones (Fig. 6b). FLL did not differ significantly among ecotypes, although terrestrial and crown species tentatively had lower values than the other ecotypes (Fig. 6c). Crown species had shorter hind limbs than all other scansorial ecotypes, while karst species had longer hind limbs than all ecotypes except for trunk ones (Fig. 6d).

Fig. 6
figure 6

Results of the post hoc comparison among ecotypes of a SVL, b AGL, c FLL, d HLL for the species of Cyrtodactylus examined in this study. Dots indicate predicted mean values, with bars showing standard deviations. Letters a-d indicate significant differences (groups not sharing the same letter differ significantly)

Full size image

Discussion

This study provides the first genus-wide examination of ecomorphology in Cyrtodactylus. Our results supported our prediction that morphological differences in body size and relative limb dimensions exist between the terrestrial and scansorial ecotypes. Terrestrial species tend to have smaller body size compared to trunk or granite species, and their hind limbs are shorter compared to karst species (Fig. 6). The crown ecotype, however, notably deviates from this pattern by sharing trait values overlapping with terrestrial and generalist species. Our prediction for generalist species was confirmed in that they show the greatest amount of variation in morphospace (Fig. 4), but they can still be discriminated from karst species based on shorter hind limbs (Fig. 6c, d), and from granite, trunk, karst and terrestrial species by a torso length which lies between the former two and the latter two (Fig. 6b).

Karst ecotypes are smaller, shorter bodied, but have relatively longer hindlimbs than granite ecotypes. The single cave adapted species included in this study overlaps only with karst species in morphospace in a position indicating long relative limb dimensions, thus tentatively supporting our predictions of showing more extreme trends in limb elongation compared to other saxicoline species (Fig. 4). This aligns with previous results presented by Grismer and Grismer (2017). The sandstone species C. kuluensis Grismer et al., (2021a) shares its morphospace with the karst, generalist and terrestrial species, by not showing pronounced changes in body size or limb proportions.

For the two arboreal ecotypes (trunk and crown), we also found that trunk species had proportionally longer limbs than crown species, although the differences are not significant for the front limbs. This indicates that in this radiation, perch diameter is positively correlated at least with hind limb length, sharing similarities with Anolis, Tropidurus, and Draco (Kohlsdorf et al., 2001; Losos, 1994; Ord & Klomp, 2014) but differing from Australian geckos (Hagey et al., 2017). Both the single intertidal species C. seribuatensis Youmans & Grismer, (2006) and the lone swamp species C. pantiensis Grismer et al., (2008) have overall smaller body size and limb proportions (particularly forelimbs). Swamp species were previously found to be similar to generalist species in their morphology (Kaatz et al., 2021), while in our analysis C. pantiensis falls outside of the morphospace occupied by the generalists (Fig. 4c, d).

Body Size Corrections

Although the overall results were comparable, the two different methods of correcting for body size yielded different results. In the PCA analysis, the GS data provided a single axis explaining more than 90% of the variance, with all traits, including body size, loading almost equally on this axis (Fig. 4b). Furthermore, the pattern of limb dimension evolution reconstructed with the GS data almost mirrors the reconstructed evolution of body size (SVL and AGL), while the RES dataset recovers a more differentiated pattern (Fig. 5a, b, Supplementary Fig. S2). This indicates that the limb dimensions have not been fully decoupled from body size. The allometric equation for the ‘GroupStruct’ function, based on Thorpe (1976), revealed the standardized variables of interest to be proportional to the species-specific mean values of the body size metric (Chan & Grismer, 2022). Thus, body size is only removed as a variable within species. Consequently, the standardized limb measurements tentatively come out larger for larger-bodied species, even if they are proportionally smaller. In contrast, our PGLS approach standardizes variables proportional to the overall variation in body size in the dataset, i.e., across species (Mahler et al., 2010). An alternative approach using the allometric equation of Thorpe (1976) would use a global mean instead of species-specific means, but that would incorporate the untested built-in assumption that all species in the genus have similar allometric growth trajectories, as in the PGLS approach, which, if wrong, would occlude signals of differences in limb dimensions (McCoy et al., 2006; Reist, 1986). Overall, this indicates that examining differences in allometric growth trajectories in Cyrtodactylus might be a promising avenue for further research.

Biomechanical Implications

Terrestrial Cyrtodactylus species appear to be prevalent in closed habitats (Agarwal & Karanth, 2015). Therefore, the short and stocky habitus of the terrestrial ecotype corroborates the hypothesis that this bauplan is beneficial in closed habitats where animals need to navigate an increased density of obstacles, such as stones, twigs, vines, or branches in the undergrowth beneath the herbaceous and shrub-layer vegetation (Arnold, 1998). These obstacles also limit maximum running speed that can be attained, given the need for frequent manoeuvres; thus, selection pressure may be less intense on running performance than it is on manoeuvrability. Quantitative data on habitat openness are absent for terrestrial Cyrtodactylus species but, given the gradual reduction of relative hind limb length in the triedrus group, further testing of the association between this ecological parameter and relative limb length in terrestrial lizards could be illustrative of a phenotype–environment correlation, and a structure–function relationship.

Biomechanical predictions for scansorial species indicate that relatively longer limbs could be beneficial in increasing limb span and overall static stability (Foster et al., 2018), whereas shorter limb spans could be beneficial in keeping the CoM closer to the substratum (Zaaf & Van Damme, 2001). Importantly, the latter prediction has, to our knowledge, not been explored experimentally, and can at least partially be questioned on theoretical grounds in that brachium and thigh length should not contribute much to the distance of the CoM from the substratum because of the sprawling posture of lizards (Birn-Jeffery & Higham, 2014; Russell & Bels, 2001). Our results support the former hypothesis, in line with many previous studies (Arnold, 1998; Foster et al., 2018; Goodman et al., 2008). In Anolis, limb length in arboreal species is positively correlated with perch diameter (Losos, 2011), although this relationship is reversed for Australian geckos (Hagey et al., 2017). Wide perch diameters should more easily allow for the maintenance of the CoM close to the substrate for taxa with greater limb spans, whereas narrow diameter perches should have greater potential for the torso to be raised farther away from the perch as limb length increases (Goodman et al., 2008). Thus, a positive correlation between limb span and perch diameter might be in line with aforementioned biomechanical predictions. The differences between the crown ecotype and the trunk and saxicoline ecotypes in Cyrtodactylus may be associated with the differential biomechanical demands imposed by the perch diameter of the locomotor surface.

Chameleons provide another exception to the norm of long limbs being associated with small diameter perches. But they possess a highly derived locomotor morphology compared to the typical lizard bauplan (Fischer et al., 2010; Peterson, 1984) including zygodactylous autopodia, which allow for a multidirectional strong grip on narrow perches (Molnar et al., 2017). The opposite correlation in Australian geckos might relate to differences in digit morphology. Anolis, for example, possesses adhesive toe pads located on the basal and central parts of the digits (Garner et al., 2019), whereas most Australian geckos have terminal, leaf-like adhesive toepads with toepad area being mainly located distally on the digits (Russell & Gamble, 2019). This latter configuration might allow for a better multidirectional grip on small branches, countering the effects of the greater distance of the CoM from the perch. Cyrtodactylus is traditionally described as a genus of padless “bent-toed, clawed geckos,” but some scansorial species have repeatedly been described as being intermediates in a morphological series from padless to basal pad-bearing forms (Russell & Gamble, 2019). Therefore, it might be expected that Cyrtodactylus ecotypes should follow similar evolutionary trajectories to Anolis or padless lizard clades. Also, lizards could alter the distance of the CoM from the substratum through locomotor behaviours, such as changing between a more sprawling and a more erect posture (Fuller et al., 2011; Irschick & Jayne, 1999). Therefore, experiments comparing the climbing performance and kinematics of crown and trunk Cyrtodactylus species on different perch diameters could prove highly informative for advancing our understanding of the biomechanical consequences of longer versus shorter relative limb length. The gradual overlap between karst and granite ecotypes in their limb morphology indicates that different rock types do not impose drastically different demands upon locomotor morphology in terms of limb dimensions or body size.

Geographic Patterns of Variation

Despite the ecomorphological signal recovered in this study, we detected some clade specific differences within ecotypes, consistent with the strong phylogenetic signal recovered for most measurements. All saxicoline species have evolved proportionally longer hind limbs and, in the case of the granite ecotype, have also evolved longer bodies and torsos (Fig. 6). This trend is, however, differentially expressed in different clades. The increased in body length is more pronounced in the pulchellus group than in any other saxicoline clade or species, and the karst associated pulchellus species also show increased body and torso lengths, contrary to the overall trend in this ecotype (Fig. 5, 6). The pulchellus group is part of the major Papuan clade which has dispersed from Indochina to Papua and subsequently diversified (A15 in Grismer et al., 2022a). From Papua, the ancestor of the pulchellus group dispersed to Sundaland (Grismer et al., 2022a). This clade also contains many of the trunk radiations such as the loriae, louisiadensis, and novaeguineae groups (Grismer et al., 2020), which also evolved large body sizes and relatively long hindlimbs (Fig. 5, 6). Previous authors have already noted that Cyrtodactylus shows a west-to-east gradient of increasing body size, with the largest forms occurring in the Papuan region (Oliver et al., 2014). As such, we consider it possible that the large body sizes characteristic of the Papuan radiation, including the pulchellus group, evolved in relation to different ecological opportunities in Papua as opposed to those of other bioregions (Losos, 2010). More generally, this geographically regional pattern implies that body size evolution in Cyrtodactylus is perhaps not only driven by adaptations associated with microhabitat differences but also by biogeographic patterns related to different ecological opportunities in different regions (Stroud & Losos, 2016).

The only noticeable contradiction from the ecomorphological patterns elucidated in this study is evident in the species of the philippinicus group. In this clade, many species that retain their ancestrally generalist habit have evolved somewhat longer extremities and increased body size, while the trunk species have retained the moderate body size and limb proportions typical of generalists (Fig. 5). The philippinicus group initially radiated in Sundaland, most likely in the area that is now Borneo, in parallel with other, closely related species groups including the malayanus group (Grismer et al., 2020). The malayanus group is an ancestrally trunk-associated radiation, whereas the trunk-associated species of the philippinicus group evolved later, from a generalist ancestor (Grismer et al., 2020). Therefore, niche space for large bodied trunk species may already have been occupied in the region of origin, causing descendants to follow a different evolutionary trajectory due to inter-specific competition (Tejero-Cicuéndez et al., 2022). Most of the large and long-legged generalist species of the philippinicus group belong to the clade that dispersed to the Philippines and radiated within the archipelago (Grismer et al., 2020). Thus, changes in morphology might also be explained by different ecological opportunities (Losos, 2010). Future studies should also take biogeography into account when considering the within- versus extra-regional evolution of morphology in Cyrtodactylus.

The habitat of Philippine species of the philippinicus group is described as being predominantly ‘riparian’; thus, individuals are often found in forest patches along rivers (Oaks et al., 2019; Siler et al., 2010; Welton et al., 2009, 2010a, 2010b). A riparian corridor specialization could offer an alternative explanation for their deviation from the generalist ecomorphology, in that this microhabitat may differ in some unknown parameter from the habitat occupied by other generalist Cyrtodactylus species in other bioregions (Grismer et al., 2020). As for many other gecko lineages, microhabitat data for Cyrtodactylus are only available for a few species (mentioned in passing in taxonomic literature) and are otherwise restricted to regional summaries based upon ‘expert’ experience (Grismer et al., 2020, 2021b). Reliance on subjective, non-quantitative characterization of habitat types could lead to either oversimplification or incorrect assignment of microhabitat categories for some species (Riedel et al., 2020). Thus, additional fine-scale, quantitatively characterized microhabitat data for Cyrtodactylus species could help to improve our understanding of the evolutionary ecomorphology of the genus (Riedel et al., 2020).

Cyrtodactylus Ecotypes

We conclude that most of the ecotypes identified by Grismer et al., (2021a, 2021b) differ in their locomotor morphology and that these differences evolved multiple times, independently, in association with changes in structural microhabitat preferences. In addition to their ecological characterisation (Grismer et al., 2021a, 2021b), therefore, we further define ecotypes for the genus Cyrtodactylus here as follows: The “Terrestrial” ecotype is defined by a short body length and short relative limb proportions; the “Trunk” ecotype, by a large body size and elongated torsos and hind limbs; the “Generalist” ecotype, by an intermediate body length and limb proportions that are substantially variable; the “Crown” ecotype by similar morphological characteristics to the Generalist, but with far less variation among species. Both, “Karst” and “Granite” ecotypes have elongated hind limbs, but karst-dwelling species have longer hind limbs than granite-dwellers, while the latter are larger with more elongated trunks than the former (Fig. 6). Additional differences in traits not considered here but recovered in previous studies include an overall depressed body form for saxicoline species (Grismer & Grismer, 2017; Kaatz et al., 2021), or a prehensile tail in crown species (Grismer et al., 2021a, 2021b). Whether the “Swamp” and “Cave” ecotypes show distinct locomotor morphologies awaits clarification in future studies. Such studies should also consider additional morphological traits which could be associated with adaptations to different microhabitats, to test if at least some of the ecotypes may constitute ecomorphs that could be distinguished solely on morphological grounds. Incorporating three dimensional geometric morphometric data could also prove informative in this context (Yuan et al., 2019).

Cyrtodactylus is clearly a well-suited genus for analysis of the relationships among morphology, ecology, and biomechanics. We have illustrated the close relationships between ecology and morphology, as well as some of the possible biomechanical consequences/drivers. However, functional experiments will be important for connecting morphology to performance under different ecologically relevant contexts. Furthermore, components of morphological diversity might also be driven by biogeographic dispersal patterns (Oliver et al., 2014), rendering Cyrtodactylus an important clade and model system for the examination of evolutionary patterns of form-function relationships under different biogeographic conditions and ecological opportunities (Stroud & Losos, 2016).

Overall, our results indicate that Cyrtodactylus might constitute an adaptive radiation. Ecotype diversity in Cyrtodactylus matches the diversity of ecomorphs in Anolis lizards, the classical model for adaptive radiations in lizards (Huie et al., 2021; Losos, 2011), and our study indicates that at least some of these ecotypes show morphological adaptations of the locomotor system. If the ecotypes can be established as ecomorphs this would further support sympatric speciation and subsequent niche partitioning as drivers of diversity in Cyrtodactylus (Stroud & Losos, 2016). Concurrently, biogeographic pattern also contributes to this diversity (Grismer et al., 2022a). This combination of drivers of biodiversity indicates that Cyrtodactylus might be a promising model for the study of adaptive radiation in regions with a complex geological history. Notably, some scansorial Cyrtodactylus have putatively evolved incipiently expressed adhesive toepads (Russell & Gamble, 2019). Adhesive toepads have been proposed as a putative key innovation for the evolution of arboreality (Miller & Stroud, 2021), and key innovations are suggested to be important drivers of adaptive radiations, along with geographic colonisation among other factors (Losos, 2010; Stroud & Losos, 2016; Yoder et al., 2010). Thus, while most previous studies on adhesive toepads as key innovations have treated these simply as a presence absence trait (e.g., Garcia-Porta & Ord, 2013; Miller & Stroud, 2021), Cyrtodactylus would provide the opportunity to examine the role of gradually evolving toepads for the diversification of this putative adaptive radiation, and to disentangle the differential contribution of a putative key innovation and geographic dispersal as drivers of adaptive radiations under the influence of differential ecological opportunities (Losos, 2010; Stroud & Losos, 2016; Yoder et al., 2010).