Karstic Landscapes Are Foci of Species Diversity in the World’s Third-Largest Vertebrate Genus Cyrtodactylus Gray, 1827 (Reptilia: Squamata; Gekkonidae)

: Karstic landscapes are immense reservoirs of biodiversity and range-restricted endemism. Nowhere is this more evident than in the world’s third-largest vertebrate genus Cyrtodactylus (Gekkonidae) which contains well over 300 species. A stochastic character mapping analysis of 10 different habitat preferences across a phylogeny containing 344 described and undescribed species recovered a karst habitat preference occurring in 25.0% of the species, whereas that of the other eight speciﬁc habitat preferences occurred in only 0.2–11.0% of the species. The tenth category—general habitat preference—occurred in 38.7% of the species and was the ancestral habitat preference for Cyrtodactylus and the ultimate origin of all other habitat preferences. This study echoes the results of a previous study illustrating that karstic landscapes are generators of species diversity within Cyrtodactylus and not simply “imperiled arks of biodiversity” serving as refugia for relics. Unfortu-nately, the immense ﬁnancial returns of mineral extraction to developing nations largely outweighs concerns for biodiversity conservation, leaving approximately 99% of karstic landscapes with no legal protection. This study continues to underscore the urgent need for their appropriate management and conservation. Additionally, this analysis supports the monophyly of the recently proposed 31 species groups and adds one additional species group.


Introduction
The dramatic topography of karstic landscapes composes some of the most surreal images of our world and has stirred the emotions of ancient artisans and natural historians for time on end. But not only are these crenulated, repeating layers of rugged terrain steeped in natural beauty (Figure 1), they are the only refuge for some of the most seriously endangered species on the planet [1]. Asia contains 8.35 million km 2 of karstic habitat with some of the most extensive concentrations ranging from China to western Melanesia ( Figure 2). These formations are notable for their fragmented, island-like nature, with hills, caves, and towers forming archipelagos of habitat-islands stretching across broad geographic areas. This, and their fractured and eroded surfaces-which provide a myriad of microhabitats in which many taxonomic groups have specialized-have contributed to their extraordinarily high degrees of range-restricted endemism [2][3][4][5]. Karst formations are often referred to as "imperiled arks of biodiversity" [5]. However, a stochastic character mapping analysis of habitat preference using 243 species of the gekkonid genus Cyrtodactylusthe third most speciose vertebrate genus on the planet-indicated just the opposite [6]. Grismer et al. [6] demonstrated that karstic landscapes not only harbor range-restricted endemics, but have been the foci of speciation for the largest independent gekkonid radiations across all of Indochina and Southeast Asia. They went on to show that even in this ecologically labile genus, karst-associated species outnumbered by threefold all other species bearing other specific habitat associations. As such, this has transformed our view of karstic landscapes from that of "limestone museums" harboring relictual endemics, to platforms of speciation and generators of biodiversity across a broad taxonomic landscape (e.g., [7][8][9][10]).  Cyrtodactylus is by far the most speciose and ecologically diverse gekkotan genus [6,11]. It currently contains 306 nominal species (as of 14 February 2021; [12]) ranging from South Asia to Melanesia ( Figure 3) where they occupy a vast diversity of habitats. As would be expected from a group this large and widely distributed, it bears a broad variety of ecotypes ranging from short robust terrestrial species to cryptically colored arboreal species to gracile cave-dwelling and karst-adapted specialists (e.g., [13][14][15][16][17][18][19][20][21]; Figure 4). The annual rate at which new species are being described is unprecedented and shows no signs of leveling off ( Figure 5) and the majority of the most recently described species are associated with karst formations. In some cases, multiple species from distantly related clades may be found throughout a single karstic archipelago [14,20], and even more remarkable, different species from distantly related clades may even occupy the same small karst formation [14,20]. The intent of this paper is to test, (1) whether or not the same clades bearing the same specific habitat preferences presented by Grismer et al. [6] are recoverable, (2) whether or not the relative frequencies of species in each habitat preference category are not significantly different than that reported by Grismer et al. [6], (3) and specifically, is the hypothesis that karstic landscapes are generators of biodiversity further supported. We test these hypotheses by augmenting Grismer et al.'s [6] original phylogeny of 243 species with an additional 101 species (a 44% increase in species coverage) and by adding a new category of habitat preference. Additionally, with this significant influx of species, we test the monophyly of the 31 different species groups recently designated by Grismer et al. [11] based on their phylogeny of 310 named and unnamed species (an 11% increase in species coverage). specifically, is the hypothesis that karstic landscapes are generators of biodiversity further supported. We test these hypotheses by augmenting Grismer et al.'s [6] original phylogeny of 243 species with an additional 107 species (a 44% increase in species coverage) and by adding a new category of habitat preference. Additionally, with this significant influx of species, we test the monophyly of the 31 different species groups recently designated by Grismer et al. [11] based on their phylogeny of 310 named and unnamed species (an 11% increase in species coverage).

Habitat Preferences and Ecotypes
Here we refine some of the criteria for designating habitat preference used by Grismer et al. [6] based on newly acquired data from recent publications and fieldwork. We also add an additional habitat preference (sandstone), bringing the total to 10 as opposed to nine categories (Table S1). Habitat preference for each species was coded as a discrete character state and ascertained by integrating data from the literature, firsthand experience of the authors, and personal communication with researchers familiar with particular species. Grismer et al. [6] acknowledged that some of these categories could be further subdivided (e.g., arboreal into branch, twig, and leaf), but those subdivisions become far less defensible owing to a lack of detailed microhabitat information. In this regard, many species can be considered data deficient, inasmuch as baseline information on their ecological requirements are often limited to anecdotal observations made at the time of their collection (e.g., [14]). The potential biases of using limited observations from a single locality at one point in time to ascertain the habitat preference of an entire species does not go unnoticed. However, in many cases, these are the only data available. Nonetheless, judiciously vetted, natural history observations summarized across the literature coupled with our own field observations and those of others, can provide a useful framework for supporting robust, testable, downstream hypotheses regarding habitat preference. The habitat preferences and their associated ecotypes bearing the same categorical names are described below. Obvious morphological correlates associated with some ecotypes are noted only for additional clarity.

1.
General ( Figure 4A). Species that use the majority of the microhabitats in their immediate surroundings in whatever environment they inhabit. The microhabitats may include rocks of all types (when present), logs, tree trunks (with or without holes and crevices), and all vegetative structures of various dimensions, the ground, and humanmade structures in many cases. No particular microhabitat is notably preferred over any other although some species may be most often observed in low vegetation. Trunk ( Figure 4B). These are species generally found on the trunks and large branches of large trees at varying heights and often take refuge in cracks, crevices, or holes in the trunks. They may occasionally occur on large granite rocks but only if the rocks are near the trees. These species are generally the largest and most robust species in the genus [22][23][24]. None have been reported to have prehensile tails although some species may coil the tail horizontally similar to that seen in arboreal species. 3.
Karst ( Figure 4C). These are generally more gracile species that are restricted to habitats where limestone rock (karst) is present. Individuals use this substrate (including cliff faces, small rocks, and boulders) as well as adjacent vegetation. If caves are present, they will enter only into the twilight zone and usually no deeper than 50 m from the entrance [14]. Despite what has been written about many karst-associated species being cave species or cave adapted (e.g., [25]), none truly are and most are more commonly found on the outside of caves (see below). These species do not occur in habitats lacking karstic substrates. 4.
Cave ( Figure 4D). These are species that occur exclusively in the cave-like environments formed by large granite boulders. Open spaces between the boulders can be quite extensive and contain areas where very little light penetrates. These species rarely occur on the out-facing (i.e., the forest-side) surfaces of the boulders and for the most part, are restricted to the spaces between the boulders at varying depths below the surface of the ground in extremely low levels of illumination. These are truly cave-adapted species with notably thin, gracile bodies, long limbs, flat heads, large eyes, and faded color patterns [13,26,27].

5.
Terrestrial ( Figure 4E). These are species that generally occur only on the ground and may take refuge beneath natural and human-made surface objects. They may occasionally be found on the tops of small rocks (when present) or on the bases of small trees and shrubs but never higher than 1 m above the ground. These species are relatively small and notably squat, with short fat tails, thick heads, and short digits [28,29]. 6.
Arboreal ( Figure 4F). These are cryptically colored species [30,31] generally restricted to small branches, leaves, trunks of varying sizes, and shrubs. Some may take refuge beneath exfoliating bark often as high or higher than three meters above the ground. These species are rarely observed on the ground or lower than 1.5 m above the ground. In such instances, it is usually during windy and/or rainy nights (perhaps forced down from higher up; [32]; authors pers. obs.) or during egg laying. All species have a prehensile tail used as a climbing aid [31][32][33] that is often carried in a coiled, elevated position. 7.
Swamp ( Figure 4G). These are species restricted to swampy habitats that use low, viny vegetation, the trunks of small trees and shrubs, or small logs often above, but always in close proximity to water. These species generally have large eyes with notably reddish-orange irises [34,35]. 8.
Granite ( Figure 4H). These are generally more robust, strongly tuberculated species found in forested habitats bearing large granite boulders (not just small, scattered, granite rocks or rocks of other types). Vegetation is often used, especially by hatchlings and juveniles, but individuals occur more commonly on the granite boulders in all planes of orientation. These species do not occur in forested areas lacking granite boulders. 9.
Intertidal ( Figure 4I). This category contains a single species that occurs exclusively in the rocky intertidal zones of small islands in the Seribuat Archipelago off the southeastern coast of Peninsular Malaysia and avoids nearby forested regions even if they lack other species of Cyrtodactylus [19,36]. 10. Sandstone ( Figure 4J). This category was not included in Grismer et al. [6]. It contains a single species endemic to a forested sandstone massif isolated in the lowlands of northwestern Cambodia [11]. This species is known to forage only on the surface or within crevices of sandstone rocks and was not observed on the nearby vegetation [37]. This species is similar in body shape to closely related granite-associated species (Grismer unpublished).

Phylogenetic Analyses
A Maximum likelihood (ML) analysis was implemented using the IQ-TREE webserver [42,43] preceded by the selection of substitution models using the Bayesian Information Criterion (BIC) in ModelFinder [44] which selected TVM+F+I+G4 for the tRNAs and codon position 1 and GTR+F+I+G4 for codon positions 2 and 3. One-thousand bootstrap pseudoreplicates via the ultrafast bootstrap (UFB; [45]) approximation algorithm were employed, and nodes having UFB values of 95 and above were considered strongly supported [46]. We considered nodes with values of 90-94 as well supported.
A Bayesian inference (BI) phylogeny was estimated using Bayesian Evolutionary Analysis by Sampling Trees (BEAST) version 2.4.6 [47] implemented in CIPRES (Cyberinfrastructure for Phylogenetic Research; [48]). Input files were constructed in Bayesian Evolutionary Analysis Utility (BEAUti) version 2.4.6 using a lognormal relaxed clock with unlinked site models, linked trees and clock models, and a Yule prior and run in BEAST version 2.4.6 [47] on CIPRES. bModelTest, implemented in BEAST, 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 [49] to ensure effective sample sizes (ESS) were well above 200 for all parameters. A maximum clade credibility tree using mean heights at the nodes was generated using TreeAnnotator v. 1.8.0 [50] with a burn-in of 1000 trees (10%). Nodes with Bayesian posterior probabilities (BPP) of 0.95 and above were considered strongly supported [51,52]. We considered nodes with values of 0.90-0.94 as well supported.
Grismer et al. [6] demonstrated that in their 243-species data set, the third codon position contributed significantly to the strongly supported topological resolution of the tree and showed no signs of codon saturation. In their 310-species tree, Grismer et al. [11] demonstrated that their mito-nuclear tree constructed from ND2 and three nuclear genes did not improve the resolution or the nodal support of the deep nodes in their ND2 tree. Therefore, only ND2 was used in this analysis.

Ancestral State Reconstruction
In order to estimate the probability of each habitat preference at each node in the tree, we employed a stochastic character mapping (SCM) analysis implemented in R [v3. 4.3] using the R package Phytools [53] on the BEAST tree converted to newick format. The transition-rate matrix that best fit the data was identified by comparing corrected Akaike Information Criterion (AICc) scores among alternate models using the R package ape 5.2 [54]. Three transition-rate models were considered: a 90-parameter model having different rates for every transition type (the ARD model); a 45-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, an MCMC approach was used to sample the most probable 1000-character histories from the posterior using make.simmap() and then summarized them using the summary() command.

Results
The ML analysis recovered essentially the same well to strongly supported tree ( Figure 6) recovered in Grismer et al. [11]. The same 31 monophyletic species groups designated in Grismer et al. [11] were recovered here even though sampling in was greatly expanded with additional species (Table S1). The ML analysis also recovered a new clade, designated here as the tibetanus group, that is composed of Cyrtodactylus tibetanus, C. cf. tibetanus, and C. zhaoermii. Cyrtodacylus cf. tibetanus and C. zhaoermii were unavailable for the analysis of Grismer et al. [11], where C. tibetanus was recovered as the earliest diverging member of the lawderanus group. Che et al. [55] recovered the same new clade in a less inclusive (i.e., fewer species) mito-nuclear phylogeny. Although Grismer et al. [11] recovered C. rubidus as the sister species of the lateralis group, it was not included in that group because this relationship was well supported only in the ML analysis and not the BI analysis. Here, it is placed in the lateralis group with high support in both analyses (90 UFB, 0.90 BPP), a grouping also supported by the fact that all members of this group have prehensile tails.  The BEAST analysis recovered a tree with generally strong nodal support throughout with a 94.4% topological consistency (recovering 322 of the same 347 nodes) as the ML tree (Figure 7). The AICc scores for the three transition-rate models were ARD = 1101.751; SYM = 1035.445; and ER = 890.9552. The results of the SCM analysis were consistent with those of Grismer et al. [6] in that the ER model recovered large and small clades that independently evolved the same habitat preferences throughout the geographic range of the genus ( Figure 7A). The SCM recovered a general habitat preference as being ancestral for not only the genus Cyrtodactylus but for all other major clades and ultimately all other habitat preferences as well. Notably for this study, however, the two largest independently evolved lineages of karst-associated species-the lineage composed of the sadansinensis, yathepyanensis, oldhami, sinyineensis, and chauquangensis groups and the angularis group-were also recovered, even with their expanded species contents. Their parapatric distributions across much of western and northern Indochina coincide with regions bearing the most extensive karstic landscapes (Figure 8). Other less diverse, independently evolved karstic lineages, such as the linnwayensis group from the Shan Plateau in Myanmar and a karst-associated subclade from the Thai-Malay Peninsula in the pulchellus group, were also recovered and associated with regions rich in karstic habitats ( Figures 7A and 8). Several isolated instances of the independent evolution of karst habitat preference are scattered across the tips of the tree, representing species from Borneo (C. cavernicolus, C. limajalur, C. muluensis), Cambodia (C. laangensis), China (Cyrtodactylus sp. SYS r1232), Indonesia (C. darmandvillei), Myanmar (C. aunglini, C. chrysopylos, C. myaleiktaung), Papua New Guinea (C. tanim), Peninsular Malaysia (C. evanquahi, C. guakanthanensis, C. gunungsenyumensis, C. metropolis, C. lenggongensis, C. sharkari), and Vietnam (C. sp. nov., C. yangbayensis) ( Figure 7A).  (13) cave (7) general (134) granite (32) intertidal (1) karst (86) sandstone (1) swamp (5) terrestrial (29) trunk (38) Node probability of ancestral habitat preference Percentage of species in each habitat    These data are consistent with those of Grismer et al. [6] in showing that the frequency of karst-associated species far out-numbers that of any other specific habitat preference and is nearly two and one-half times more prevalent than any other specific habitat preference in that it contains 25.0% of the species followed by trunk (11.0%), granite (9.2%), terrestrial (8.4%), arboreal (3.8%), cave (2.0%), swamp (1.4%), and intertidal and sandstone (0.2%; Figure 7B). In Grismer et al. [6], granite-associated species comprised the second highest habitat preference and trunk-associated species the third. That ranking has been reversed here. The percentage of species with a karst habitat preference was 29.6% in Grismer et al. [6] but dropped to 25.0% here. We posit that this drop of nearly 5% is a direct result of our inability to explore unsurveyed karstic regions on the Shan Plateau and in the Salween Basin of Myanmar during 2020 due to COVID-19.

Discussion
The analysis presented here is based on the most complete phylogeny of the genus Cyrtodactylus to date with an increase of 101 species from that of Grismer et al. [6] and 35 from that of Grismer et al. [11]. The hypotheses marshaled by Grismer et al. [6] concerning the evolution of habitat preference is supported here in that there was no notable change in the frequencies of species bearing different habitat preferences across the genus-even with the addition of 107 species. More specifically, however, a karst habitat preference retained a higher frequency than that of any other specific habitat preference (25.0% versus 0.2-11%), supporting the hypothesis that these landscapes are platforms for the generation of biodiversity. This pattern is particularly strong in Indochina and less so on islands throughout the Indo-Australian Archipelago, reflecting the sharp contrast in the extent of karstic landscapes between these regions ( Figure 8). These data clearly underscore the importance of karstic habitats to this hyper-diverse genus and continue to amplify the work of many other authors indicating that the high levels of biodiversity and rangerestricted endemism in karstic habitats rivals that of most other habitats throughout the tropics (see discussions in [1,4,5,10,[56][57][58][59][60][61]). The sad irony is that, although these are some of the most imperiled ecosystems on the planet due to unregulated and unsustainable quarrying practices, only 1% of these terrains throughout Asia are afforded any form of legal protection. Therefore, the diversity of the karst-associated species in general-and Cyrtodactylus in particular-are, for the most part, without legal protection. Unfortunately, the immense financial returns from cement manufacturing makes the challenge of karst conservation difficult and many governments from developing nations that are willing to overlook sustainable quarrying policies in order to expand their economy [1]. Continued exploitation of karstic habitats for limestone shows no signs of abating.

Conclusions
This study echoes the results of Grismer et al. [6] in that karstic landscapes are exceedingly important for maintaining Cyrtodactylus diversity and serve as foci for their speciation and maintenance of their diversity. Referring to them as "imperiled arks of biodiversity" is somewhat misleading as these are ecological platforms for speciation that not only continue to generate the most speciose, independent, radiations of the Gekkota, but do so across a broad range of other taxonomic groups (e.g., [7][8][9][10]62]). Referring to them as "imperiled arks of biodiversity" instead of centers for speciation draws attention away from their importance as generators of biodiversity in an era of biodiversity crisis and could potentially lessen the urgency for legislative conservation measures.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/d13050183/s1, Table S1: Species, habitat preference with supporting references, species group designations, and GenBank accession numbers for specimens used in the SCM analysis. Species can be cross-referenced to Figure 6 by their GenBank no.