A new granite cave-dwelling Bent-toed Gecko from Viet- nam of the Cyrtodactylus irregularis group (Squamata; Gekkonidae) and a discussion on cave ecomorphology

An integrative analysis of genetic, morphological, and ecological data recovered a new granite cave-adapted species, Cyrtodactylus raglai sp. nov., from the Song Giang River Valley, Khanh Hoa Province, Vietnam. Cyrtodactylus raglai sp. nov. is nested within one of two major clades within the irregularis species group where it forms a monophyletic group with C. cryptus and its sister species C. kingsadai. It differs from its sister species by an uncorrected pairwise sequence divergence of 16.5% and 16.8% based on the mitochondrial genes NADH dehydrogenase subunit 2 and its flanking tRNAs and the first subunit of cytochrome oxidase c (CO1), respectively. Cyrtodactylus raglai sp. nov. is a narrow-range endemic restricted to a riparian, granite cave microhabitat and its overall morphology bears that of other granite cave-dwelling ecomorphs in the genus. The Song Giang River Valley and its associated waterways are currently threatened by the construction of a hydropower station which will negatively impact the surrounding ecosystem. Urgent field surveys in this region are necessary in order to obtain critical data on its biodiversity and importance to conservation efforts in southern Vietnam.


Introduction
The Cyrtodactylus irregularis group is the largest, most taxonomically diverse species group within the genus. It contains 22 nominal species and at least 11 undescribed species (Grismer et al 2021) that range throughout cen-tral and southern Vietnam, eastern and southern Laos, and is represented by isolated populations in northern and eastern Cambodia (Fig. 1). The diversity of the irregularis group is centered in Vietnam wherein 18 of the 20 Vertebrate Zoology 71, 2021, 155-174 | DOI 10.3897/vz.71.e60225 Copyright Anh The Nguyen. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
nominal species occur, most of which are narrow-range endemics (Orlov et al 2007;Nazarov et al 2008Nazarov et al , 2012Ngo and Bauer 2008;Geissler et  . The taxonomic diversity of the irregularis group is mirrored by its wide range in habitat preferences, having species adapted to karstic and granitic forests, tree trunks, and generalist that range through many different habitats (see above references). Grismer et al (2020) demonstrat- Figure 1. Distribution of the species of the Cyrtodactylus irregularis group. Dot in the center of an icon indicates the type locality of a species. The small inset map shows Vietnam and the location of the area of studies detailed in a larger map (red rectangle). For details on species localities see Pauwels et al. (2018); Neang et al. (2020); and Ostrowski et al. (2020). ed that species within the irregularis group bearing the same habitat preference did not form monophyletic subgroups, but were nested in different places throughout the group's phylogeny. Nonetheless, the statistically inferred ancestral habitat preference for the group was that of a habitat generalist from which all other habitat preferences ultimately evolved (Grismer et al 2020). Reported here is a newly discovered population of the irregularis group from Khanh Hoa Province, Vietnam. Genetic data from the mitochondrial gene NADH dehydrogenase subunit 2 and its flanking tRNAs (ND2) and the first subunit of cytochrome oxidase c (CO1), indicate it is not nested within any other species of the irregularis group. Morphological analyses demonstrate it bears several statistically significant mensural and meristic differences from its closest relatives, and an ecological assessment of its lifestyle indicates it is the first member of the irregularis group with a granite cave habitat preference. Integrating the results of these different analyses leads us to the robust hypothesis that this population is a new species and as such, is described below.

Morphological data and analyses
Morphological data included both meristic and mensural characters. Measurements were taken on the left side of the body when possible to the nearest 0.1 mm by NAP using Mitutoyo dial calipers under a Nikon SMZ 1500 dissecting microscope and follow Grismer and Grismer (2017) and Grismer et al (2018). Measurements taken were: snout-vent length (SVL), taken from the tip of the snout to the vent; tail length (TL), taken from the vent to the tip of the tail; tail width (TW), taken at the base of the tail immediately posterior to the postcloacal swelling; sternal width (SW), distance between the ventral edges of the opposing glenoid fossae; forearm length (FL), taken on the ventral surface from the posterior margin of the elbow while flexed 90° to the inflection of the flexed wrist; forelimb width (FLW), measured from the anterior to the posterior margins of the brachium immediately adjacent to their insertion points on the body; tibia length (TBL), taken on the ventral surface from the posterior surface of the knee while flexed 90° to the base of the heel; hind limb width (HLW), measured from a point equidistant between its anterior and posterior insertion points on the body to the tip of the fourth toe; axilla to groin length (AG), 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; pelvic width (PW), distance between the dorsal tips of the ilia; pelvic height (PH), distance from the dorsal tip of the ilium to the ventral surface of the pubis; 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; head depth (HD), the maximum height of head measured from the occiput to base of the lower jaw; eye diameter (ED), the greatest horizontal diameter of the eye-ball; eye to ear distance (EE), measured from the anterior edge of the ear opening to the posterior edge of the bony orbit; snout length (SNT), measured from anterior-most margin of the bony orbit to the tip of snout; eye to nostril distance (EN), measured from the anterior margin of the bony orbit to the posterior margin of the external nares; interorbital distance (IO), measured between the anterior-most edges of the bony orbits; ear length (EL), measured as the greatest vertical distance of the ear opening; and internarial distance (IN), measured between the nares across the rostrum.
Meristic characters evaluated were the number of supralabial scales (SL) counted from the largest scale immediately below the middle of the eyeball to the rostral scale; infralabial scales (IL) counted from the mental to the termination of enlarged scales just after the upturn of the mouth; the number of paravertebral tubercles (PVT) between limb insertions counted in a straight line immediately left of the vertebral column; the number of longitudinal rows of body tubercles (LRT) counted transversely across the center of the dorsum from one ventrolateral fold to the other; the number of longitudinal rows of ventral scales (VS) counted transversely across the center of the abdomen from one ventrolateral fold to the other; the number of expanded subdigital lamellae on the fourth toe (4TL) counted from the base of the first phalanx where it contacts the body of the foot to the claw including the claw sheath (see Murdoch et al 2019: Fig. 2); and the total number of enlarged femoral scales (FS) from each thigh combined as a single metric. In some species, only the distalmost FS are greatly enlarged and the proximal scales are smaller whereas in others, the enlarged scales are continuous with the enlarged precloacal scales and the separation of the two scales rows was determined to be at a point even with the lateral body margin (see Murdoch et al 2019: Fig. 3). The number of enlarged precloacal scales (PS); the number of precloacal pores in (PP) in males; the number of rows of enlarged post-precloacal scales (PPS) on the midline between the enlarged precloacal scales and the vent; the number of postcloacal tubercles (PCT) on either side of the base of the tail; and the number of dark body bands (BB) between the nuchal loop (the dark band running from eye to eye across the nape) and the hind limb insertions. Other morphological characters evaluated were whether or not the enlarged femoral scales and precloacal scales were separated by a diastema at the base of the femora or contiguous and if so, if the proximal femoral scales were less than one-half the size of the distal femorals (see Murdoch et al 2019: Fig. 3); and the presence or absence of a pocket in the skin webbing between the digits of the hind and forefeet (see Murdoch et al 2019: Fig. 2).
Color pattern characters evaluated were the nuchal loop being continuous from orbit to orbit across the nape, or separated medially into paravertebral halves, bearing an anterior azygous notch or not, and the shape of its posterior border; presence or absence of a band on the nape; dorsal body bands bearing paired, having paravertebral elements or not; dark dorsal body bands wider than light interspaces, with or without lightened centers, edged with light tubercles or not, irregularly shaped or more regular (straight-edged or even); dark markings present or absent in the dorsal interspaces; top of head bearing combinations of dark, diffuse mottling or dark, distinct blotches overlain with a light-colored reticulating network or not; anterodorsal margin of thighs and brachia whitish due to a lack of dark pigment; light caudal bands bearing dark markings or immaculate; light caudal bands encircling the tail or not; and dark caudal bands wider than light caudal bands. Majority-rule consensus tree from 1000 ML bootstrap pseudoreplicates of Cyrtodactylus. Phylogeny is based on 1469 bp of ND2 Black circles represent nodes the UB and BPP support > 90 and 0.90, respectively. Grey circles represent nodes with UB support > 90 only. The white circle represents a node with BPP support > 0.90 only. Nodes lacking circles are unsupported. Scale bar denotes uncorrected pairwise sequence divergence. Species group designations follow Grismer et al. (2021).Photograph by Anh The Nguyen.
Analyses of variance (ANOVA) were conducted on selected meristic and mensural characters (see below) with 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 a subset of the irregularis group. Character means bearing statistical differences were subjected to a Tukey-HSD 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 among the sampled individuals from the selected subset of species were visu-alized using principal component analysis (PCA) from the ADEGENET package in R (Jombart et al 2010) implemented by the prcomp() command. PCA is a dimension reducing algorithm that decreases the complexity of a data set by finding a subset of input variables that contain the most relevant information (i.e. the greatest variance in the data) while de-emphasizing those characters that do not, thus increasing the overall accuracy of the results by eliminating noise and the potential of overfitting (Agarwal et al 2007). PCA is an unsupervised analysis that recovers morphospatial relationships among the sampled individuals (i.e. data points) and how well they form uncoerced clusters that may or may not align themselves with putative species boundaries delimited by phylogenetic analyses and defined by univariate anal- yses (see below). It is important to note that clusters of conspecific individuals were not pre-defined in the analysis, but simply color-coded in the scatter plot in order to observe their positions and morphospatial relationships with respect to one another. Given that the meristic data ranged from 8-50, these data were log-transformed prior to analysis in order to normalize their distribution so as to ensure characters with very large or very low values could not over-leverage the results owing to intervariable nonlinearity. To remove potential effects of allometry in the mensural characters, size was normalized using the following equation: X adj =log(X)-β[log(SV-L)-log(SVL mean )], where X adj =adjusted value; X=measured value; β=unstandardized regression coefficient for each population; and SVL mean =overall average SVL of all populations (Thorpe 1975(Thorpe , 1983Turan 1999;Lleonart et al 2000). The morphometrics of each species were normalized separately and then concatenated so as not to conflate intra-with interspecific variation (Reists 1986). All data were scaled to their standard deviation to insure they were analyzed on the basis of correlation and not covariance.
A discriminant analysis of principal components (DAPC) was performed on both data sets. The DAPC places individuals from each predefined population into separate clusters (i.e., plots of points) bearing the smallest within-group variance that produce linear combinations of centroids having the greatest between-group variance (i.e. linear distance; Jombart et al 2010). DAPC relies on scaled data calculated from its own internal 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-95% of the variation in the data set (Jombart and Collins 2015) as determined from a scree plot generated as part of the analysis. Retaining too many PCs forces false structure to appear in the data while retaining too few, runs the risk of missing true structure (Cangelosi and Goriely 2007).
The type material is deposited in the herpetological collection of the Department of Zoology, Southern Institute of Ecology (SIEZC) in Ho Chi Minh City, Vietnam, and the Zoological Museum of Lomonosov Moscow State University, Moscow, Russia (ZMMU).

Genetic data and analyses
A data set composed of the 310 species of Cyrtodactylus from Grismer et al (2021) plus one specimen of Cyrtodactylus sp. from Song Giang River Valley, Khanh Hoa Province, Vietnam (ZMMU R16688), was assembled and aligned in Geneious. A phylogeny was constructed based on 1469 base pairs of ND2. Mediodactylus russowii and Hemidactylus frenatus were used as outgroups to root the tree based on Gamble et al (2012).
Additionally, to compare the newly discovered population of Cyrtodactylus sp. from Khanh Hoa Province with other members of the irregularis species group, we generated a 672 base pair segment of the 5'-end of (CO1) which included 29 sequences of 26 species and compared it with the CO1-barcoding dataset from Pauwels et al (2018) and Neang et al (2020). Cyrtodactylus intermedius was used as an outgroup to root the tree based on Wood et al (2012. Since at present, the CO1 dataset for members of the irregularis species group has better taxon coverage than the ND2 dataset, we herein applied the CO1 fragment primarily as a DNA-barcoding marker (Murphy et al 2013; Brenan et al 2017), and did not rely on it for reconstructing phylogenetic relationships within the genus. Cyrtodactylus intermedius CO1 sequence was used as an outgroup to root the tree based on Grismer et al (2021). Due to a significant lack of overlap of specimens sequenced for both genes, we elected not to concatenate the gene fragments from the different data sets.
Genomic DNA was isolated from liver tissue stored in 95% ethanol using the following DNA extraction protocol. A small amount of liver tissue (about 1-3 mg) was dried and transferred to 500µl of extraction buffer including 1% SDS, 0.4 mg/ml proteinase K, 10mM Tris HCl (pH 8.0), 1mM EDTA, 100mM NaCl, and subsequently incubated at 56 °C for 3 hours, then centrifuged at 12,000 rpm for one minute. The 450 µl of supernatant were transferred into a new 1.5 ml tube, where 50 µl of 3M sodium acetate solution and 500µl of ice-cold isopropanol were added, mixed and refrigerated at -20 °C for 15 minutes, again centrifuged again at 12,000 rpm for one minute, the supernate was discarded, and the DNA sediment was washed with 500 µl of 70% ethanol. After the ethanol was removed, the total DNA was dried at 50 °C for 8 minutes and subsequently resuspended with 50 µl TE.
We also amplified the 5'-fragment of the CO1 mtDNA gene; the PCR parameters and the amplification protocol and conditions were the same as for the ND2 gene, but the annealing temperature was at 45 °C. Primers used both for PCR and sequencing were the VF1-d (5'-TTCT-CAACCAACCACAARGAYATYGG-3') and the VR1-d (5'-TAGACTTCTGGGTGGCCRAARAAYCA-3') following Ivanova et al (2006).
All PCR products were visualized using 1.0 % agarose gel electrophoresis. Successful PCR products were sent to National Key Laboratory of Gene Technology (Institute of Biotechnology, VAST, Hanoi, Vietnam) for PCR purification, cycle sequencing, sequencing purification, and sequencing using the same primers as in the amplification step. Sequences were analyzed from both the 3' and the 5' ends separately to confirm congruence between reads. Forward and reverse sequences were uploaded and edited in Geneious 2019.0.4 (https://www.geneious.com). Following sequence editing we aligned the protein-coding region of ND2 and the flanking tRNAs using the MAFTT v7.017 (Katoh and Kuma, 2002) plugin under the default settings in Geneious 2019.0.4 (https://www.geneious. com). Mesquite v3.04 (Maddison and Maddison, 2015) was used to calculate the correct amino-acid reading frame and to confirm the lack of premature stop codons in the ND2 portion of the DNA fragment. The 5'-end partial CO1 sequences were aligned in MAFTT v7.017, and similarly checked for the absence of stop-codons in Mesquite v3.04.
The sequence of ND2 gene and flanking tRNAs fragment and the 5'-end partial sequence of CO1 gene of the newly collected specimen ZMMU R16688 from Khanh Hoa Province, Vietnam, were submitted to GenBank with the accession numbers MW675652 and MW675653, respectively [available upon manuscript acceptance]. Gen-Bank accession numbers for the ND2 sequences of the remaining Cyrtodactylus species in the dataset and the outgroups are listed in Grismer et al (2021).

Phylogenetic analyses
A Maximum likelihood (ML) analysis for the ND2 data set was implemented using the IQ-TREE webserver (Nguyen et al 2015;Trifinopoulos et al 2016) preceded by the selection of substitution models using the Bayesian Information Criterion (BIC) in Model-Finder (Kalyaanamoorthy et al 2017) which selected TVM+F+I+G4 for codon position 1, TIM3+F+I+G4 for codon position 2, GTR+F+ASC+G4 for codon positon 3, and GTR+F+I+G4 for the tRNAs. The best-fitting model selected for CO1 dataset was GTR+G for the first and the third codon positions and HKY+I for the second codon position as suggested by the Akaike Information Criterion (BIC). One-thousand bootstrap pseudoreplicates via the ultrafast bootstrap (UB; Hoang et al 2018) approximation algorithm were employed, and nodes having UB values of 95 and above were considered strongly supported (Minh et al 2013). We considered nodes with values of 90-94 as weakly supported. MEGA7 (Kumar et al. 2016) was used to calculate uncorrected pairwise sequence divergence among the ingroup species using the complete deletion option which removes gaps and missing data from the alignment prior to the anlysis.
A Bayesian phylogenetic tree (BI) was estimated using Bayesian Evolutionary Analysis by Sampling Trees ( . 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 300,000,000 generations and logged every 30,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 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 (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 (Huelsenbeck et al 2001;Wilcox et al 2002). We considered nodes with values of 0.90-0.94 as weakly supported.

Species delimitation
The general lineage concept (GLC: de Queiroz 2007) adopted herein proposes that a species constitutes a population of organisms evolving independently from other such populations owing to a lack of gene flow. By "independently," it is meant that new mutations arising in one species cannot spread readily into another species (Barraclough et al 2003;de Queiroz 2007). Under the GLC implemented herein, molecular phylogenies were used to recover monophyletic mitochondrial lineages of individual(s) (i.e. populations) in order to develop initial species-level hypotheses-the grouping stage of Hillis (2019). Discrete color pattern data and morphological data were then used to search for unique characters and patterns consistent with the previous designations of the species-level hypotheses-the construction of boundaries representing the hypothesis-testing step of Hillis (2019)-thus providing independent diagnoses to complement the molecular analyses.

Phylogenetic data
For the ND2 dataset, the ML and BI analyses recovered trees with very similar topologies, and the ML topology used here (Fig. 2), is the same as that in Grismer et al (2021). The Cyrtodactylus irregularis group is composed of two clades (clade 1 and 2) which are strongly supported in the ML analysis (UB 99 and 96, respectively) but only clade 1 is strongly supported in the BI analysis (BPP 0.96) . The new population from Song Giang River Valley is nested within clade 1 and forms a well-supported monophyletic group in both analyses with C. cryptus and  its sister species C. kingsadai (100/1.00). Together, they form the strongly supported (99/1.00) sister lineage to the remaining species of clade1. The Song Giang River Valley population bears a 13.3% uncorrected pairwise sequence divergence from its sister species C. kingsadai and a 16.5% divergence from C. cryptus. From the remaining species in clade 1, its difference ranges from 8.1-17.0% (Table 1). The CO1 barcoding data set recovered a polytomy with no backbone support and little resemblance to the ND2 tree (Fig. 3). In both trees however, the new population from the Song Giang River Valley formed a monophyletic group with C. cryptus, an undescribed population from Ba Na, Danang Province, Vietnam, and its sister species C. kingsadai, although none of these relationships were supported in the CO1 tree. The uncorrected pairwise sequence divergence using CO1 between the Song Giang River Valley population and its sister species C. kingsadai was 16.8%, clearly indicating these two populations are not conspecific.

Meristic data
The PCA, based on the meristic data set, demonstrated that the new population from Song Giang River Valley clustered separately from its closest relatives and overlapped C. kingsadai, along PC2 (Fig. 4A). PC1 accounts for 43.8% of the variation in the data and loads most heavily for longitudinal rows of tubercles, fourth toe lamellae, and enlarged femoral scales (LPT, 4TL, and FS, respectively; Fig. 4C; Table 2). PC2 accounts for an additional 26.1% of variation in the data set and loads most heavily for labial scales (SL and IL) and ventral scales (VS). Similarly, the Song Giang River Valley population is well-separated from all other species in the meristic DAPC with no overlap of its 67% inertia ellipses (Fig.  4B). The ANOVAs and subsequent TukeyHDS tests of the meristic data demonstrated that the Song Giang River Valley population bears statistically different mean values between it and C. cryptus and C. kingsadai across various combinations of characters ( Fig. 4D; Table 3).

Mensural data
The PCA, based on the mensural data set, demonstrated that the Song Giang River Valley population clustered separately from its closest relatives and overlapped with them along PC2 but not PC1 (Fig. 5A). PC1 accounted for 72.9% of the variation in the data set and loaded most heavily for head length, head width, and eye diameter (HL, HW, ED, respectively; Fig. 5C; Table 4). PC2 accounted for an additional 15.8% of the variation and loaded most heavily for head depth and snout length (HD and SNT, respectively). Similarly, the Song Giang River Valley population is well-separated from all species in the mensural DAPC with no overlap of its 67% inertia ellipses (Fig. 5B). The ANOVAs and subsequent TukeyHDS tests of the mensural data demonstrated that the Song Giang River Valley population bears statistically different mean values between it and C. cryptus and C. kingsadai across various combinations of characters ( Fig. 5D; Table 3). Based on the phylogenetic, PCA, DAPC, and ANO-VA analyses and the high degree of uncorrected pairwise sequence divergence between the Song Giang River Val-ley population and the remaining species in clade 1 (Figs 1, 2, 3), we hypothesize that the Song Giang River Valley population is a discretely diagnosable lineage which shows no evidence of reticulation with any other lineage and as such, should be accorded species status.  Paratypes. Two adult females, SIEZC 2.0243 and ZMMU R16688 bearing the same data as the holotype except that the latter was collected on 16 August 2020.

Diagnosis.
Cyrtodactylus raglai sp. nov. can be separated from all other species of clade1 of the C. irregularis group by having 10 or 11 supralabials; nine or 10 infralabials; 44-47 paravertebral tubercles; 14 or 15 rows of longitudinally arranged tubercles; 36-39 ventrals; 8-10 expanded subdigital lamellae, 12 or 13 unexpanded subdigital lamellae, and 21-22 total subdigital lamellae on the fourth toe; 18 or 19 enlarged femorals; 12 enlarged precloacals; four rows of enlarged post-precloacals; three postcloacal tubercles in males; five precloacal pores in the male; no pitted precloacal scales in females; enlarged femorals and enlarged precloacals not continuous; proximal femorals less than one-half size of distal femorals; enlarged subcaudals; maximum SVL 111.7 mm; small, irregularly shaped dark blotches on top of head; and four irregularly shaped body bands edged with white tubercles wider than the interspaces (Tables 5, 6). (Figs. 6, 7). Adult male SVL 95.0 mm; head moderate in length (HL/SVL 0.28) and width (HW/HL 0.66), flattened (HD/HL 0.38), distinct from neck, triangular in dorsal profile; lores concave anteriorly, weakly inflated posteriorly; prefrontal region slightly concave; canthus rostralis rounded; snout elongate (ES/HL 0.42), flat, rounded in dorsal profile; eye large (ED/HL 0.25); ear opening narrow, elliptical, obliquely oriented, moderate in size; eye to ear distance slightly greater than diameter of eye; rostral rectangular, partially divided dorsally by inverted Y-shaped furrow, bordered posteriorly by large left and right supranasals and one small azygous internasal, bordered laterally by first supralabials; external nares bordered anteriorly by rostral, dorsally by large supranasal, posteriorly by two moderately sized postnasals, bordered ventrally by first supralabial; 10(R,L) rectangular supralabials extending to below midpoint of eye, second supralabial slightly larger than first; 9(R,L) infralabials tapering smoothly to slightly past the termination of enlarged supralabials; scales of rostrum and lores flat, same size as granular scales on top of head and occiput; scales of occiput intermixed with   small, distinct, tubercles; superciliaries elongate, largest anteriorly; mental triangular, bordered laterally by first infralabials and posteriorly by large left and right trapezoidal postmentals contacting medially for ~60% of their length posterior to mental; one row of slightly enlarged, elongate sublabials extending posteriorly to fourth infralabial; gular and throat scales small, granular, grading posteriorly into slightly larger, flatter, smooth, imbricate, pectoral and ventral scales. Body relatively short (AG/SVL 0.43) with well-defined ventrolateral folds; dorsal scales small, granular, interspersed with small, rounded, semi-regularly arranged, smooth tubercles; tubercles extend from occiput onto base of tail forming transverse rows; approximately 14 longitudinal rows of tubercles at midbody; approximately 46 paravertebral tubercles; 36 flat, imbricate, ventral scales much larger than dorsal scales; five large, pore-bearing, precloacal scales; no deep precloacal groove or depression; and four rows of large post-precloacal scales on midline.

Description of holotype
Forelimbs thin, relatively long (FL/SVL 0.16); lacking tubercles, granular scales slightly larger than those on body; palmar scales rounded, slightly raised; digits well-developed, inflected at basal interphalangeal joints; digits slightly narrower distal to inflections; subdigital lamellae transversely expanded, those proximal to joint inflections wider than those distal to inflection; claws well-developed, sheathed by a dorsal and ventral scale; hind limbs thin, more robust than forelimbs, long (TBL/ SVL=0.20), covered dorsally by granular scales interspersed with slightly larger, weakly keeled tubercles and anteriorly by flat, slightly larger scales; ventral scales of thigh flat, imbricate, larger than dorsals; subtibial scales large, flat, imbricate; one row of 10(R)9(L) enlarged femoral scales terminating distally before knee, not continuous with enlarged precloacal scales; proximal femoral scales much smaller than distal femorals, the latter forming an abrupt union with smaller, rounded, ventral scales of posteroventral margin of thigh; femoral pores absent; plantar scales flat; digits well-developed, inflected at basal interphalangeal joints, slightly narrower distal to inflections; subdigital lamellae transversely expanded, those proximal to joint inflections wider than those distal to inflection, 10(R,L) transversely expanded subdigital lamellae on fourth toe proximal to joint inflection that extends onto the sole; 12(R,L) narrower lamellae distal to inflection; 22 total subdigital lamellae; and claws well-developed, sheathed by a dorsal and ventral scale at base.
Tail long (TL/SVL 1.25), original, 119 mm in length, 9.0 mm in width at base, tapering to a point; dorsal caudals small, generally square; median row of subcaudals transversely expanded, significantly larger than dorsal caudals, not extending up onto lateral side of tail; transverse rows of 2-4 keeled tubercles on anterior one-half of tail, paravertebral tubercles largest; tubercle rows separated by 7-9 rows of dorsal caudals; base of tail bearing small hemipenal swellings with three large postcloacal tubercles on either side; and postcloacal scales flat, imbricate. (Figs. 6, 7). Ground color of top of head, limbs, and dorsum gray; top of head heavily mottled with irregularly shaped, small, dark-brown blotches, light-colored reticulate pattern absent; dark-brown, nuchal  loop extends from posterior margin of one orbit to posterior margin of other orbit, anterior and posterior margins on nape jagged, no anterior azygous notch; no dark banding on nape; four dark-brown, irregularly shaped, dorsal body bands reaching ventrolateral folds extending from shoulders to pelvis, wider than light-colored interspaces, bearing lightened vertebral areas, edged in light-colored tubercles; interspaces bearing dark markings; limbs dark- brown with irregularly shaped, light-colored markings; eight, wide, dark-brown caudal bands much wider than seven; light-colored caudal bands do not encircle tail and have darkened centers; iris gold bearing thin, black reticulations; edges of pupils orange; venter beige with faint, dark mottling on lateral edges of belly, undersides of limbs; and subcaudal region dark-brown.

Coloration in life
Variation (Fig. 7). The paratypes resemble the holotype very closely in all aspects of coloration and pattern. Spec-  imens ZMMU R16688 and SIEZC 2.0243 have eight and nine light-colored caudal bands, respectively, as opposed to seven in the holotype and SIEZC 2.0243 has 10 darkbrown caudal bands as opposed to eight in the holotype. Meristic and mensural variation is presented in Table 6.
Distribution (Fig. 1). Cyrtodactylus raglai sp. nov. is known only from the type locality in the Song Giang River valley, Khanh Trung Commune, Khanh Vinh District, Khanh Hoa Province, Vietnam.
Etymology. The new species name "raglai" is given in a reference to the Raglai people, an ethnic group living in the forested mountain areas of Khanh Hoa Province of Vietnam, including the Song Giang River Valley where the new species was found. In Raglai language, the self-designating word "raglai" also means "forest", stressing the importance of the tropical forest ecosystem for this people. To reflect this polysemy, the new species name is given as a noun in apposition and hence is invariable.
Natural history. Cyrtodactylus raglai sp. nov. was recorded in the forested valley of Song Giang River in the northwestern part of Khanh Hoa Province (Fig. 1), in the vicinity of the Song Giang Hydropower Station. The surrounding mountains form the northeastern slopes of Langbian Plateau and are covered with polydominant montane evergreen tropical forest and are dissected by a rich network of small streams and rivulets that feed into the Song Giang River (Fig. 8A). Scattered outcroppings of large granite boulders (ca. 3-5 m in height, 6-8 m in width) occur throughout the river valley and many of these have small permanent streams running through and beneath them. Here, soil erosion has formed small narrow caves among the boulders varying from 3-5 m in depth (Fig. 8B, C) and C. raglai sp. nov. was only observed inside such caves or within crevices between the boulders.
Lizards were elusive and difficult to find and only four specimens were recorded during three excursions to the area, and all were restricted only to specific outcroppings. Lizards were active from 19:00 h to 21:00 h, and were usually observed while foraging on the vertical surfaces of the boulders (Fig. 8C)

Discussion
Evolution of cave ecomorphology Grismer and Grismer (2017) demonstrated that in the Cyrtodactylus condorensis group, the granite cave-adapted species C. eisenmanae and C. grismeri, had a distinctive, gracile overall morphology adapted for climbing on vertical surfaces in low levels of illumination that was very different from the stout, robust morphology of the closely related sympatric habitat generalists C. condorensis and C. leegrismeri. The cave ecomorphs, C. eisenmanae and C. grismeri, have shorter trunks; longer limbs; a longer head and snout; larger eyes; and a flatter pelvis and head.
In a genus-wide analysis involving 243 species, Grismer et al (2020) statistically demonstrated that a granite cave habitat preference evolved independently from a general habitat preference at least twice in the Cyrtodactylus badenensis -condorensis group lineage, once in C. hontreensis of the intermedius group (sec. Grismer et al 2021), and at least one more time from a karst-adapted ancestor in C. nigriocularis of the angularis group (sec. Grismer et al 2021). Although no mensural data have yet been taken to bear out the cave-adapted morphology of C. badenensis, C. hontreensis, or C. nigriocularis, it is clear from their overall body stature, they bear many of the characters listed above for C. eisenmanae and C. grismeri (Fig. 9). Cyrtodactylus raglai sp. nov. is the sixth species of Cyrtodactylus to have a cave habitat preference and the first species to do so in the irregularis group. Morphometric comparisons of C. raglai sp. nov. to other closely related members of the irregularis group could not be evaluated using the complete data set of Grismer and Grismer (2017) because those data for the other irregularis group species do not exist. Therefore, only morphometrics involving head shape and eye size were evaluated. Despite small samples sizes, C. raglai sp. nov. differed statistically from the habitat generalist C. cryptus in having a longer and wider head, longer snout, and larger eyes (Fig. 5D)-characters observed in the cave ecomorphs, C. eisenmanae and C. grismeri. Its granite-adapted sister species C. kingsadai, has an intermediate morphology between C. raglai sp. nov. and C. cryptus and differs statistically from C. raglai sp. nov. only by having a shorter snout. However, given its intermediate position in the PCA, DAPC, and the boxplots (Fig. 5), a larger sample size would likely increase their statistical differences. Nonetheless, based on the limited morphometric data and ecological observations made at the time of collection, we hypothesize that C. raglai sp. nov. is a cave ecomorph similar in morphology to all the other cave ecomorphs (Fig. 9). A study using a complete morphometric data set comparing all the cave ecomorphs to each other and to their closest relatives is currently under way (Grismer et al in progress).

Conservation
The description of Cyrtodactylus raglai sp. nov. continues to underscore the fact that the mountainous areas of southern Vietnam harbor more species of Cyrtodactylus than any other area in the country. This recent discovery of yet another narrow-range endemic, is a clear indication that the true diversity of this area is not yet known. This makes it all the more important the urgent need for continued field work in this region, and in particular, forested areas with granitic outcroppings where more narrow-range endemics are likely to be discovered. To date, Cyrtodactylus raglai sp. nov. is known only from a narrow area within the Song Giang River Valley. In 2014, a significant part of this valley was partially flooded after the construction of the Song Giang Hydropower Station. A second hydropower station on the Song Giang River is currently under construction and is expected to be completed in 2022. Economic growth in Khanh Hoa Province requires additional electricity, and the plans for further development of hydropower stations on the Song Giang River will likely have a strong impact on the hydrological regime of the river and the surrounding forest ecosystems. Currently, this forested area has no legal protection. Therefore, an urgent need for additional herpetological surveys in the Song Giang River Valley is crucial for estimating its biodiversity and importance for nature conservation in southern Vietnam. Given the relatively small estimated range of Cyrtodactylus raglai sp. nov. and the increasing threats to the ecosystem where it was found, further research is urgently required to clarify the extent of its distribution, population trends, and conservation status. Because C. raglai sp. nov. is known only from three specimens, we tentatively suggest it should be categorized as Data Deficient (DD) according to the IUCN's Red List categories (IUCN Standards and Petitions Committee, 2017).