Convergent evolution of karst habitat preference and its ecomorphological correlation in three species of Bent-toed Geckos (Cyrtodactylus) from Peninsular Malaysia

By studying ecomorophology in the context of phylogeny, researchers can parse out similarity due to common ancestry versus that due to convergence. This is especially true among relatively closely related species where both phylogenetic and environmental constraints may be operating simultaneously. We explored these issues among three karst-associated species from two lineages of Cyrtodactylus—the sworderi group from Peninsular Malaysia and the swamp clade from Peninsular Malaysia and western Indonesia of the agamensis group. A stochastic character mapping analysis using five different habitat preferences corroborated a larger previous analysis in recovering a general habitat preference as an ancestral condition for all habitat preferences and a karst habitat preference in C. guakanthanensis and C. gunungsenyumensis of the sworderi group and C. metropolis of the swamp clade as convergently evolved. Multivariate and univariate analyses of 10 morphometric characters revealed that the ecomorphological similarity between C. guakanthanensis and C. gunungsenyumensis of the sworderi group was also convergent. The ecomorphology of C. metropolis of the swamp clade was intermediate between a karst-adapted ecomorphology and a swamp-generalists ecomorphology. Of the 10 morphometric characters employed in this analysis, only three—head length, head width, and forelimb width—showed any signs of phylogenetic signal. Cyrtodactylus metropolis is hypothesized to be a recently refuged swamp-dwelling species that frequented the Batu Caves environments prior to urbanization of the surrounding swamp habitat to which it is now confined.


Introduction
The concept that an animal's form has evolved in response to the way it navigates its habitat underpins the study of ecomorphology-the intersection of organismal morphology, life history, and adaptation ( Van der Klaauw 1948;Wainwright and Reilly 1994). Studying ecomorphology in a large monophyletic group, enables researchers to decipher among morphological similarities based on common ancestry versus those generated independently from similar selection pressures in similar environments (e.g., Baxter et al. 2008;Gross et al. 2009;Losos 2009;Chan et al. 2010;Mahler et al. 2013;Toyama 2017;Grismer et al. 2020a). Within diurnal lizards, ecomorphology has been well-studied in a number of lineages of skinks, anoles, and tropidurines (Melville and Swain 2000;Bergmann and Irschick 2010;Losos 2010;Lee et al. 2013;Pincheira-Donoso andMeiri 2013, 2015;Grismer et al. 2018;Toyama 2017). These studies demonstrated that suites of morphological characteristics have not only adapted their constituent species to their respective habitats, but contributed to diverse, adaptive radiations within these lineages-even though in tropidurines, changes in some of these characters were constrained more by phylogeny than habitat (Toyama 2017).
Only a few such studies (e.g., Grismer and Grismer 2017;Grismer et al. 2015Grismer et al. , 2020aNielson and Oliver 2017) have been done on nocturnal Bent-toed Geckos (Genus Cyrtodactylus) despite their ecological and morphological diversity. Cyrtodactylus is the third largest vertebrate genus in the world with well over 306 nominal species (Sitthivong et al. 2019;Utez et al. 2021;Grismer et al. 2021) that collectively range from South Asia to Melanesia. Across its vast distribution, species of Cyrtodactylus inhabit a multitude of environments such as karstic and granitic landscapes, swamps, intertidal zones, caves, forest floors, and various arboreal microhabitats (Grismer et al. 2020b). This ecological diversity is mirrored in their morphological diversity which includes small, stout, short-limbed terrestrial species; large, robust, trunk-dwelling species; and flattened, elongate, gracile cave-dwelling species to name a few (see Grismer et al. 2020b). In Peninsular Malaysia, there are 32 species of Cyrtodactylus that also occur in a wide variety of environments and exhibit a broad diversity of body shapes (Grismer 2011;Utez et al. 2021). Previous studies have indicated that a number of habitat preferences evolved multiple times in a number of species and that karstic landscapes in particular, were invaded multiple times by species from various clades (Grismer et al. 2014a, b). Moreover, there have even been multiple invasions of karst habitats by different species within the same clade (Grismer et al. 2012(Grismer et al. , 2014bQuah et al. 2019;Wood et al. 2020). However, the morphology of these karst-associated species has never been examined, and thus, it is not known if their convergence in habitat preference also coincides with a convergence in morphology. To address this issue, we conducted a comparative phylogenetic analysis on two non-reciprocally monophyletic lineages from western Sundaland-the C. sworderi group and the C. agamensis group (sec. Grismer et al. 2021)-both of which contain karst-associated species that are endemic to Peninsular Malaysia. The goals of this study are to assess the relationship between habitat preference and the morphology of karst-associated species within and between each clade in order to ascertain 1) if morphology is or is not correlated to habitat preference in each, 2) if there are ecomorphological similarities that have convergently evolved among species within and between each clade, and 3) if so, is the evolution of the ecomorphology in these lineages influenced by phylogeny.

Methods
We ran a Bayesian Evolutionary Analysis by Sampling Trees (BEAST) (see below) using the genetic data set of Grismer et al. (2021: Table 2 and A1), containing 1469 base pairs of the mitochondrial gene NADH dehydrogenase subunit 2 (ND2) and its flanking tRNAs and 310 described and undescribed species of Cyrtodactylus (see below). This analysis was run in order to test for topological consistency among the Maximum likelihood (ML) and Bayesian Inference (BI) trees of Grismer et al. (2021) with the BEAST analysis herein. The BEAST tree was subsequently employed in a stochastic character mapping (SCM) analysis (see below) to test for the corroboration of habitat preference evolution between it and that of a smaller data set of 243 species from Grismer et al. (2020b). A section of the BEAST tree was then used here as our ingroup sample that included 31 species (one individual per species) from four species groups (sec. Grismer et al. 2021) that comprised the most exclusive monophyletic lineage from western Sundaland that contained the three karst-associated species that were being compared: the marmoratus group (nine species), the lateralis group (three species), the sworderi group (five species), and the agamensis group (14 species). Within the ingroup, two non-reciprocally monophyletic lineages-the sworderi group from Peninsular Malaysia and Sumatra and swamp clade of the agamensis group from Peninsular Malaysia, Singapore, and Natuna Besar, Indonesia-contained the three karst-associated species of interest (Grismer et al. 2021: Fig. 3).

Morphological data and statistical analyses
Data were taken from all species in the sworderi group (Cyrtodactylus sworderi (n=5), C. gunungsenyumensis (n=10), C. guakanthanensis (n=22), C. tebuensis (n=12), C quadrivirgatus (n=10)) and four available species in the swamp clade (Cyrtodactylus majulah (n=3), C. metropolis (n=6), C. payacola (n=7), and C. pantiensis (n=13)). Measurements were taken on the left side of the body when possible to the nearest 0.1 mm using Mitutoyo dial calipers under a Nikon SMZ 1500 dissecting microscope and follow Grismer and Grismer (2017). Measurements taken were: snout-vent length (SVL), taken from the tip of the snout to the vent; axilla to groin length (AXG), taken from the posterior margin of the forelimb at its insertion point on the body to the anterior margin of the hind limb at its insertion point on the body; head length (HL), the distance from the posterior margin of the retroarticular process of the lower jaw to the tip of the snout; head width (HW), measured at the angle of the jaws; eye diameter (ED), the greatest horizontal diameter of the eye-ball; eye to snout distance or snout length (SNT), measured from anterior most margin of the bony orbit to the tip of snout; pelvic width (PW), distance between the lateral surfaces of the dorsal tips of the ilia; pelvic height (PH), distance from the dorsal tip of the ilium to the ventral surface of the pubis; hind limb length (HDL), measured from a point equidistant between its anterior and posterior insertion points on the body to the tip of the fourth toe; forelimb width (FLW), measured from the anterior to the posterior margins of the brachium immediately adjacent to their insertion points on the body; and forelimb length (FLL), measured from a point equidistant between its anterior and posterior insertion points on the body to the tip of the fourth finger.
Analyses of variance (ANOVA) were conducted on characters (see below) with normalized data and statistically similar variances (i.e., p values ≤ 0.05 in a Levene's test) to search for the presence of statistically significant mean differences (p < 0.05) among species across the data set. Characters bearing statistical differences were subjected to a TukeyHSD test to ascertain which species pairs differed significantly from each other for those particular characters. Boxplots were generated in order to visualize the range, mean, median, and degree of differences between pairs of species bearing statistically different mean values. All statistical analyses were performed in R [v3.4.3].
The morphospatial clustering of the sampled individuals was visualized using principal component analysis (PCA) from the ADEGENET package in R (Jombart et al. 2010). To remove potential effects of allometry in the mensural characters, only adults were used (determined by minimum size of gravid females or gonadal inspection) and variation in adult size was normalized using the following equation: X adj =log(X)-β[log(SVL)-log(SV L mean )], where X adj =adjusted value; X=measured value; β=unstandardized regression coefficient for each population; and SVLmean=overall average SVL of all populations (Thorpe 1975(Thorpe , 1983Turan 1999;Lleonart et al. 2000). The metrics of each species were normalized separately so as not to conflate intra-with interspecific variation (Reists 1986). The data were scaled to their standard deviation to ensure they were analyzed on the basis of correlation and not covariance. The adjusted data and their summary statistics are in the supplementary material (Tables 1S, 2S).
A discriminant analysis of principal components (DAPC) from the ADEGENET package in R was also performed. DAPC relies on scaled data calculated from its own PCA as a prior step to ensure that variables analyzed are not correlated and number fewer than the sample size. Dimension reduction of the DAPC prior to plotting, is accomplished by retaining the first set of PCs that account for approximately 90% of the variation in the data set (Jombart and Collins 2015) as determined from a scree plot generated as part of the analysis.
Related species in a clade can resemble one another because of common ancestry (phylogenetic signal), evolutionary convergence (no phylogenetic signal), random evolution, or various degrees of each. To test for the presence or absence of phylogenetic signal in the morphological data set, a mean value for each scaled trait was calculated and tested independently by employing the phylosig() command in the R package phytools (Revell 2012) and generating Blomberg's K and Pagel's lambda (λ) statistics (Blomberg et al. 2003). Theoretically, the K and λ statistics compute values from 0 to infinity where a K = 1 indicates that species resemble one another as much as would be expected under a model of Brownian motion (BM)-something similar to random drift. A K > 1 indicates that species are more similar to one another than would be expected under a BM model-more phylogenetic signal due to closeness of relationship. A K < 1 indicates that there is less phylogenetic signal in the data and that species are more similar to one another than would be expected under a BM model which could be due to adaptive evolution not correlated with phylogeny (Blomberg et al. 2003)-convergence. A λ close to zero indicates the phylogenetic signal in the data is equivalent to that expected if the data arose on a star phylogenythat is, no phylogenetic signal and species resemble one another due to convergence. λ = 1 corresponds to a BM model where there is phylogenetic signal in the data and species resemble one another based on closeness of relationship. 0 < λ < 1 is somewhere in between. Significant differences of both statistics to the values of zero (K or λ = or close to 0) or one (K or λ = 1) were also tested.

Phylogenetic analyses
A BEAST 2 analysis version 2.4.6 (Drummond et al. 2012) was implemented in CIPRES (Cyberinfrastructure for Phylogenetic Research; Miller et al. 2010). An input file was constructed in BEAUti (Bayesian Evolutionary Analysis Utility) version 2.4.6, a lognormal relaxed clock with unlinked site models, linked trees and clock models, and a Yule prior. Model choice was done through bModelTest (Bouckaert and Drummond 2017) implemented in BEAUti, was used to numerically integrate over the uncertainty of substitution models while simultaneously estimating phylogeny using Markov chain Monte Carlo (MCMC). MCMC chains were run for 350,000,000 generations and logged every 35,000 generations. The BEAST log file was visualized in Tracer v. 1.6.0 (Rambaut et al. 2014) to ensure effective sample sizes (ESS) were above 200 for all parameters. A maximum clade credibility tree using mean heights at the nodes was generated using TreeAnnotator v.1.8.0 (Rambaut and Drummond 2013) with a burn in of 1000 trees (10%). Nodes with Bayesian posterior probabilities (BPP) of 0.95 and above were considered strongly supported (Hulsenbeck et al. 2001;Wilcox et al. 2002).
Habitat preference (general, granite, karst, arboreal, and terrestrial, sec. Grismer et al. 2020b) was mapped onto the tree using stochastic character mapping (SCM) implemented in the R package Phytools (Revell 2012) in order to derive probability estimates of the ancestral states at each node in the tree. A transition rate matrix was identified that best fit the data by comparing likelihood scores among alternate models using the Akaike Information Criterion (AIC) in the R package APE [v.3.4.3] (Para dis and Schilep 2018). Three transition rate models considered were: a 12-parameter model having different rates for every transition type (the ARD model); a six-parameter model with equal forward and reverse rates between states (the symmetrical rates SYM model); and a single rate parameter model that assumes equal rates among all transitions (ER). Lastly, a MCMC approach was used to sample the most probable 1000 trait histories from the posterior using the make.simmap() command and then summarized them using the summary() command.
Museum acronyms are LSUHC -La Sierra University Herpetological Collection, Riverside California, USA.

Results
The BEAST 2 analysis recovered the same relationships among the species groups as those in the ML and BI analyses of Grismer et al. (2021: Fig. 1A). However, the sister group relationship between the marmoratus group and the lateralis group + sworderi group and the monophyly of the agamensis group had weak support. However, an expanded BEAST data set with 346 species recovered the same topology as that herein with strong ML and BI support at all nodes (Grismer et al. 2021). The likelihood scores for the three transition rate models in the SCM analysis were SYM = -5.593, ER = -3.555 and ARD = -4.074. Based on the best likelihood score of the ER model, the SCM analysis corroborated Grismer et al. (2020b) in that a general habitat preference is ultimately the ancestral state for all other habitat preferences as well as for all species groups. Furthermore, it corroborated that the evolution of a karstic habitat preference evolved independently three times in the two focus groups in Peninsular Malaysia (Fig. 1)-twice in the sworderi group in C. gunungsenyumensis and C. guakanthanensis and once in the swamp clade (a subclade within the agamensis group) in C. metropolis (Fig. 1). Additionally, the general habitat preference in C. jarakensis evolved from a swamp habitat preference and the arboreal habitat preference of C. durio  and C. lateralis of the lateralis group, the granite habitat preference of C. tiomanensis, and the swamp habitat preference in C. semenanjungensis in the agamensis group, also evolved from a general habitat preference (Fig. 1). The morphological analyses of the sworderi group demonstrated that the karst-associated species C. gunungsenyumensis and C. guakanthanensis and the habitat generalists C. quadrivirgatus, C. sworderi, and C. tebuensis only overlapped slightly in the PCA and did not overlap in the DAPC ( Fig. 2A, B) in accordance with their morphology and not their phylogenetic relationships. Principal component (PC) 1 accounted for 31% of the variation and loaded most heavily for SNT, FLL, HDL, HL while PC2 accounted for 16% of the variation, loaded most heavily for AXG (Fig 2A; Table 2). Cyrtodactylus gunungsenyumensis and C. guakanthanensis show no morphological separation from one another nor is there any significant separation among C. quadrivirgatus, C. sworderi, and C. tebuensis in the PCA. The morphological separation between the karst-associated species and the habitat generalists is even more pronounced in the DAPC where the 66% confidence ellipses between C. gunungsenyumensis and C. guakanthanensis nearly completely overlap, as well as those of C. quadrivirgatus and C. sworderi. That of C. tebuensis is separated from all others but in close approximation to those of C. quadrivirgatus and C. sworderi (Fig. 2B). The ANOVA analyses recovered no significant difference between C. gunungsenyumensis and C. guakanthanensis in any character. However, each species bears significant morphological differences in all characters  among varying combinations of the habitat generalists C. quadrivirgatus, C. sworderi, and C. tebuensis, as a result of the karst-associated species having longer snouts (SNT), narrower pelvises (PW), and longer limbs (FLL and HDL) ( Table 1 Similar to the karst-adapted species in the sworderi group, the PCA and DAPC of the swamp clade species indicate that the karst-associated Cyrtodactylus metropolis is well-differentiated from all other members of that clade (Fig. 4A, B). PC1 accounted for 35% of the variation and loaded most heavily for PH, SNT, HL, HDL, ED, HW while PC2 accounted for 16% of the variation and loaded most heavily for AXG and FLL ( Fig. 4A; Table 3). The ANOVA analyses indicate that C. metropolis bears a number of significant differences between its sister species C. payacola as well as C. pantiensis but there are no significant differences between it and C. majulah (Table 1). However, this is most likely the result of the small sample size of C. majulah (n = 3), as the means of PW, PH, AXG, HL, SNT, HDL, FLW, and FLL between the two species are widely separated (Fig. 5).
Cyrtodactylus jarakensis and C. rosichonarieforum were unavailable for study. As with C. gunungsenyumensis and C. guakanthanensis, C. metropolis bears a number of significant morphological differences among other members of the swamp clade owing to its longer forelimbs (FLL) and flatter and narrower pelvis (PH and PW: Table 1). Although not significant, the boxplot demonstrates it also has a relatively longer snout (Fig. 5).
To test whether or not the karst-adapted species-Cyrtodactylus metropolis of the swamp clade and C. gu nungsenyumensis and C. guakanthanensis of the sworderi group converged on each other morphologically, a PCA and DAPC were performed on a concatenated data set composed of all the species of each group. The analyses showed that the karst-adapted C. gunungsenyumensis and C. guakanthanensis grouped together as before and that the swamp-adapted species (C. majulah, C. payaco la, and C. pantiensis) and habitat generalist (C. quadrivir ga tus, C. sworderi, and C. tebuensis) also grouped together but separately from the karst-adapted species (Fig. 6A, B). In both analyses, however, C. metropolis plotted intermediately between the two groups, showing only slight overlap among both in the PCA but no overlap of its 66% confidence ellipse with either     group in its intermediate position in the DAPC. PC1 accounted for 30% of the variation, loading most heavily for HDL, SNT, HL, and FLL whereas PC2 loaded most heavily for FLW and PW and accounted for 16% of the variation ( Fig. 6A; Table 4). These analyses indicate that species with a swamp and general habitat preference have the same body shape that differs significantly in a number of morphological characters from that of the karst-adapted species (Fig. 7). The intermediate position of C. metropolis is more clearly visualized by coding the individuals in a concatenated data set for habitat preference (i.e. karst, general, swamp) separate from C. metropolis (Figs. 8A, B, 9). K and λ statistics for the mean values of each trait mirrored one another and indicated there was no phylogenetic signal in PW, PH, AXG, SNT, ED, HDL, and FLL but both statistics indicated that HL showed phylogenetic signal in both groups and FLW and HW were phyloge-    (Table 5).
These data indicate that although head length (HL) and head width (HW) are under the influence of phylogeny, an increase in snout length (SNT) is not. A similar relationship among head shape was noted in the cave-adapted species of the condorensis group (sec. Grismer and Grismer 2017) where head length remained constant but relative snout length increased as a result in the more posterior displacement of the orbit. Overall, the K and λ statistics indicate that the morphological diversity within these lineages is not greatly influenced by phylogeny. This result is corroborated by the SCM that demonstrated that the three karst-associated species converged on habitat preference independently.

Discussion
The analyses demonstrated that the evolution of a karstadapted ecomorphology evolved independently in three species of Cyrtodactylus from Peninsular Malaysia from two unrelated lineages, and that it co-evolved with a karst habitat preference. Grismer et al. (2020b) demonstrated that throughout Indochina, karstic habitats have been a unique environment that has driven the independent evolution of some of the most diverse radiations in the genus. Studying the evolution of karstic habitat preference in non-sister species among these much smaller Sundaic lineages, provides an opportunity to assess, and test, its direct correlation with morphology in relatively short evolutionary time scales. Although C. gunungsenyumensis, C. gua kan tha nensis, and C. metropolis evolved a karst ic habitat preference independently, they converged mor pho lo gi cally in much the same way, indicating that selection pressures in a karstic habitat are strong enough to force convergence among unrelated species. Other distantly related karst-adapted species bear many of the same charac teristic as those seen in C. gunungsenyumensis, C. gua kan tha nensis, and C. metropolis-namely, longer limbs and snout, and narrower and flatter pelvises (Neilson and Oliver 2017; Grismer et al. 2020a). Interestingly however, in the concatenated analysis, Cyrtodactylus metropolis does not unequivocally group with the karst ecomorphs of the sworderi group nor with the swamp-dwellers or generalist of the swamp clade but instead, is somewhat intermediate among them-although distinctively closer to the karst-adapted species (Fig. 6A, B). Grismer et al. (2014c) noted that C. metropolis was very common on the karst vegetation surrounding the exterior walls of the limestone karst tower it inhabits (the Batu Caves in Selangor, Peninsular Malaysia) but has only rarely been observed within the cave itself. This prompted them to posit that C. metropolis is not a cave-adapted species. Instead, we hypothesize here, that prior to the urbanization of swamp habitat that once surrounded Batu Caves, C. metropolis was simply a swamp species inhabiting this portion of Peninsular Malaysia. Urbanization then transformed Batu Caves and its associ-ated vegetation into a habitat island that became the only place left wherein C. metropolis could survive. Prior to urbanization, C. metropolis may have resembled the other swamp clade species in morphology and color pattern and likely frequented the karst vegetation surrounding Batu Caves as does the habitat generalist, C. quadrivirgatus, that frequents various karstic habitats throughout its wide distribution in Peninsular Malaysia (Grismer 2011). However, we hypothesize, that subsequent to being refuged to this habitat-island, it began undergoing anagenetic change from having a swamp-adapted and generalist ecomorphology and color pattern to acquiring a more karst-adapted morphology and color pattern (i.e., wide, somewhat faded, regularly shaped dorsal interspaces) seen in C. gunungsenyumensis and C. guakanthanensis (compare Figs 10 and 11).    Appendix 2 Table S2. Summary statistics of the scaled morphometric data of swamp clade and sworderi group species. SD = standard deviation, n = sample size.