Integrative taxonomy of the rock-dwelling gecko Cnemaspis siamensis complex (Squamata, Gekkonidae) reveals a new species from Nakhon Si Thammarat Province, southern Thailand

Abstract The rock-dwelling gecko genus Cnemaspis is one of the most species-diverse genera of gekkonid in Thailand. Earlier studies relied on morphological data to identify species, but cryptic morphology often obscured species diversity in Cnemaspis. In this study, an integrative taxonomic approach based on morphological characters and sequences of the mitochondrial NADH dehydrogenase subunit 2 (ND2) gene were used to clarify current taxonomy of the Cnemaspis siamensis complex and delimit a new species from Lan Saka District, Nakhon Si Thammarat Province, southern Thailand. Cnemaspis lineatubercularissp. nov. is distinguished from other congeneric species by the combination of morphological characters: (1) maximum snout-vent length (SVL) of 40.6 mm (mean 38.8 ± SD 1.4, N = 12) in adult males and maximum SVL of 41.8 mm (mean 39.5 ± SD 1.9, N = 7) in adult females; (2) 8–9 supralabial and infralabial scales; (3) gular, pectoral, abdominal, and subcaudal scales keeled; (4) rostral, interorbitals, supercilium, palmar scales, and ventral scales of brachia smooth; (5) 5–6 small, subconical spine-like tubercles present on flanks; (6) 19–21 paravertebral tubercles linearly arranged; (7) 27–29 subdigital lamellae under the fourth toe; (8) 4–7 pore-bearing precloacal scales, pores rounded arranged in chevron shape and separated only in males; (9) one postcloacal tubercles each side in males; (10) ventrolateral caudal tubercles present anteriorly; (11) caudal tubercles restricted to a single paravertebral row on each side; (12) single median row of subcaudal scales keeled and lacking enlarged median row; and (13) gular region, abdomen, limbs and subcaudal region yellowish only in males. Genetically, the uncorrected pairwise divergences between the new species and their congeners in the C. siamensis group were between 15.53–28.09%. The new species is currently known only from granitic rocky streams at Wang Mai Pak Waterfall in the Nakhon Si Thammarat mountain range. Its discovery suggests that additional unrecognized species of Cnemaspis may still occur in unexplored areas of southern Thailand.


Introduction
The rock-dwelling gecko genus Cnemaspis Strauch, 1887 is one of the most speciose genera in the family Gekkonidae. The genus is geographically widespread from tropical Africa eastward through South Asia, southward to Southeast Asia (Bauer et al. 2007;Grismer et al. 2014). However, recent molecular phylogenetic analyses of Cnemaspis suggest the genus may be polyphyletic, with three separate, unrelated clades consisting of African, South Asian, and Southeast Asian clades (Gamble et al. 2012;Pyron et al. 2013).
Southeast Asian Cnemaspis is a monophyletic group (Gamble et al. 2012;Pyron et al. 2013) that contains 59 species distributed from Laos, southern Vietnam, Cambodia, Thailand, southward through the Thai-Malay Peninsula to Borneo, Java, and Sumatra (Bauer and Das 1998;Das 2005;Bauer et al. 2007;Grismer and Ngo 2007;Grismer et al. 2009Grismer et al. , 2014Grismer and Chan 2010;Kurita et al. 2017;Riyanto et al. 2017;Uetz et al. 2019). In Thailand, there are currently 17 recognized species of Cnemaspis (Grismer et al. , 2014Wood et al. 2017;Ampai et al. 2019;Uetz et al. 2019), ranging from Kanchanaburi Province, western Thailand  to Chanthaburi Province, eastern Thailand (Bauer and Das 1998), southward through southern Thailand and its offshore islands (Grismer et al 2014;Wood et al. 2017;Ampai et al. 2019). Based on the combination of morphological characters and molecular data, Grismer et al. (2014) indicated that Cnemaspis species from Thailand belong to four species groups, consisting of the affinis group, the chanthaburiensis group, the kumpoli group (= Pattani clade of Grismer et al. 2014) and the siamensis group.
More than half of Southeast Asian Cnemaspis species have been described primarily or solely on the basis of morphological characteristics (Smith 1925;Taylor 1963;Bauer and Das 1998;Das and Leong 2004;Das 2005;Bauer et al. 2007;. However, morphological data alone has been insufficient for resolving some taxonomic issues that are confounded by morphological crypsis. During the past decade, an integrative taxonomic approach that uses multiple sources of data (e.g., morphology, DNA sequencing, ecology, biogeography, behavior) to delimit species and describe taxa (Dayrat 2005) has been shown to be very effective in revealing cryptic species diversity and microhabitat specialization in Southeast Asian Cnemaspis (e.g., Grismer et al. 2013Grismer et al. , 2014Wood et al. 2013Wood et al. , 2017Kurita et al. 2017;Ampai et al. 2019).
During fieldwork in October 2016 and January 2019, we collected specimens of the C. siamensis group at Wang Mai Pak Waterfall, Lan Saka District, Nakhon Si Thammarat Province, southern Thailand that could not be referred to any named species. We examine qualitative and quantitative (univariate and multivariate analyses) variation in morphology and mitochondrial DNA sequence data and show that the Lan Saka specimens differ from all other species of Cnemaspis. On the basis of this integrative approach, we described the Lan Saka population as a new species.

Sampling
Specimens of Cnemaspis were collected by hand during the day (1000-1800 h) and at night (1900-2200 h) between October 2016 and January 2019 from Wang Mai Pak Waterfall, Lan Saka District, Nakhon Si Thammarat Province, Thailand. Liver samples for genetic analysis were taken from euthanized specimens and preserved in 95% ethanol. Specimens were then fixed in 10% formalin and later transferred to 70% ethanol for permanent storage. Specimens and tissue samples were deposited in the herpetological collection at the Zoological Museum of Kasetsart University, Bangkok, Thailand (ZMKU) and the Thailand Natural History Museum, Pathum Thani, Thailand (THNHM). Comparative material was also examined in the holdings of these institutions (Appendix 1), and comparative data were obtained from the original descriptions of other Cnemaspis species in Thailand Wood et al. 2017;Ampai et al. 2019).

Morphological measurements
The following morphometric measurements were taken by the first author on the left side of preserved specimens to the nearest 0.1 mm using digital calipers under a Nikon SMZ 745 dissecting microscope. Morphological measurements were taken only from adult individuals as determined by the presence of secondary sexual characteristics including the presence of hemipenes or pore-bearing precloacal scales in males, and the presence of calcium glands or eggs in females. Sixteen morphological measurements were taken following Grismer et al. (2014) and Wood et al. (2017): snout-vent length (SVL), taken from tip of snout to the anterior margin of vent; tail width (TW) at the base of the tail immediately posterior to the postcloacal swelling; tail length (TL), as distance from the vent to the tip of the tail, whether original or regenerated; forearm length (FL), taken on the dorsal surface from the posterior margin of the elbow while flexed 90° to the inflection of the flexed wrist; tibia length (TBL), taken on the ventral surface from the posterior surface of the knee while flexed 90° to the base of the heel; head length (HL), as distance from the posterior margin of the retroarticular process of the lower jaw to the tip of the snout; head width (HW) at the angle of the jaws; head depth (HD), as the maximum height of head from the occiput to the throat; axilla-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; eye diameter (ED), as the maximum horizontal diameter of the eyeball; eye-ear distance (EE), measured from the anterior margin of the ear opening to the posterior edge of the eyeball; ear length (EL), taken from the greatest vertical distance of the ear opening; eye-nostril distance (EN), measured from the anterior most margin of the eyeball to the posterior margin of the external nares; eye-snout distance (ES), measured from the anterior margin of the eyeball to the tip of snout; inner orbital distance (IO), as the width of the frontal bone at the level of the anterior edges of the orbit; and internarial distance (IN), measured between the medial margins of the nares across the rostrum.
Meristic characters of scales and qualitative observations of other structures were made through a Nikon SMZ 745 dissecting microscope. The external observations of meristic characters were taken following Grismer et al. (2014) and Wood et al. (2017): number of supralabial (SUP) and infralabial (INF) scales, counted from below the middle of the orbit to the rostral and mental scales, respectively; texture of scales on the anterior margin of the forearm; number of paravertebral tubercles (PVT) between limb insertions, counted in a straight line immediately left of the vertebral column; general size (i.e., strong, moderate, weak) and arrangement (i.e., random or linear) of dorsal body tubercles; number of subdigital lamellae beneath the fourth toe (= 4 th toe lamellae), counted from the base of the first phalanx to the claw; presence or absence of a row of enlarged, widely spaced, tubercles along the ventrolateral edge of the body flank between limb insertions; number, orientation and shape of pore-bearing precloacal scales; relative size of subcaudal and subtibial scales; and number of postcloacal tubercles on each side of tail base.

Morphological analysis
Statistical analyses were used to compare differences in size and shape in the siamensis group, including the Lan Saka samples (N = 19) and four congeners in the siamensis group: C. adangrawi (N = 8), C. chanardi (N = 7), C. omari (N = 5) and C. siamensis (N = 8). Other species in the siamensis group (C. huaseesom, C. kamolnorranathi, C. phangngaensis, C. punctatonuchalis, C. roticanai, C. thachanaensis, and C. vandeventeri) were not included in the morphometric analyses due to lack of specimens. Five putative operational taxonomic units (OTUs) were assigned on the basis of observed variation in morphometric analysis. Fifteen morphometric variables (SVL, TW, FL, TBL, HL, HW, HD, AG, EE, ED, EL, EN, ES, IO, and IN) were corrected for differences in ontogenetic composition by the following allometric equation: X adj = X β(SVL -SVL mean ), where X adj is the adjusted value of the morphometric variable; X is the original value; β is the within-clade coefficient of the linear regression of each original character value (X) against SVL; SVL = snout-vent length; SVL mean = overall average SVL length of OTUs (Thorpe 1975(Thorpe , 1983Turan 1999;Lleonart et al. 2000). Tail length (TL) was not included due to the differences in length between original and regenerated tails. Univariate analyses were implemented in the statistic software PAST 3.24 (Hammer et al. 2001) using an analysis of variance (ANOVA) to compare morphological differentiation in traits among the Lan Saka samples and the five congeners in the siamensis group. ANOVAs having p-value less than 0.05 were subjected to a Tukey's honestly significant difference (HSD) test to identify all pairwise comparisons among sample means for significant differences (p < 0.05).
Multivariate analyses were performed using the base statistical software in RStudio v. 1.2.1335 (RStudio Team 2018). A principal component analysis (PCA) using the built-in R functions: prcomp (R Core Team 2018) and ggplot2 (Wickham 2016) were performed to find the best low-dimensional space of morphological variation in data. Principal components (PCs) with eigenvalues greater than 1.0 were retained in accordance to the criterion of Kaiser (1960). A discriminant analysis of principal components (DAPC) was performed using the adegenet function (Jombart 2008) to characterize clustering and distance in morphospace. The DAPC was used for all congeners to find the linear combinations of morphological variables that have the greatest between-group variance and the smallest within-group variance. The DAPC relies on data transformation using PCA as a prior step to ensure that variables included in the discriminant analysis (DA) are uncorrected and number fewer than the sample size (Jombart et al. 2010).

Genetic analysis
Genomic DNA was extracted from liver tissue of five individuals of Cnemaspis (Table 1) using the Qiagen DNAeasy tissue kit (Valencia, CA, USA). A 1,251 bp fragment of mitochondrial (mt) DNA consisting of the NADH dehydrogenase subunit 2 (ND2) gene and the flanking tRNAs Trp, Ala, Asn, and Cys was amplified using polymerase chain reaction (PCR) under the following conditions: initial denaturation at 95 °C for 2 min, followed by a second denaturation at 95 °C for 35 sec, annealing at 52 °C for 35 sec, followed by a cycle extension at 72 °C for 35 sec, for 33 cycles using the light strand primer L4437b (5'-AAGCAGTTGGGCCCATACC-3'; Macey et al. 1997) and heavy strand primer H5934 (5' AGRGTGCCAATGTCTTTGTGRTT-3'; Macey et al. 1997). PCR products were purified using the AccuPrep PCR Purification Kit (Bioneer, Daejeon, Korea), and were sequenced using the amplifying primers on an ABI 3730 automatic sequencer (Applied Biosystems, CA, USA). Sequences were edited and aligned using Geneious R11 (Biomatters, Ltd, Auckland, New Zealand).
All new sequences were deposited in GenBank under accession numbers MT112890-MT112894 (Table 1).

Phylogenetic analysis
Homologous sequences of 69 Cnemaspis, and the outgroups Cyrtodactylus bokorensis and Hemidactylus garnotii based on Bauer et al. (2008) and Grismer et al. (2015a), were downloaded from GenBank and aligned to the five newly generated Cnemaspis sequences using Geneious R11 (Biomatters, Ltd, Auckland, New Zealand). The aligned dataset was partitioned into four partitions consisting of ND2 codon positions and tRNAs. Phylogenies were reconstructed with the maximum likelihood (ML) criterion using IQ-TREE 1.6.7 (Nguyen et al. 2015) on the IQ-TREE web server (Trifinopoulos et al. 2016). The best-fit model of substitution for each partition was estimated using IQ-TREE's ModelFinder function (Kalyaanamoorthy et al. 2017) under the Akaike Information Criterion (AIC). The selected models were TIM+F+I+G4 for first, second and third codon partitions, and HKY+F+G4 for the tRNA partition. Bootstrap analysis was performed using the ultrafast bootstrap approximation (Minh et al. 2013) with 1,000 replicates and 0.95 minimum correlation coefficient.
Phylogenies were also reconstructed with Bayesian Inference (BI) in MrBayes v3.2 on XSEDE on the Cyberinfrastructure for Phylogenetic Research (CIPRES; Miller et al. 2010) computer cluster. The best-fit model of substitution was estimated for each partition with jModelTest 2.1.10 (Posada 2008) under AIC. The selected models were GTR+ I+G4 for each ND2 codon partition, and HKY+ I+G4 for the tRNA partition. Two simultaneous runs, each with three heated and one cold chain, were performed using the default priors for 10,000,000 generations, with trees sampled every 1,000 generations from the Markov Chain Monte Carlo (MCMC). Runs were halted after the average standard deviation split frequency was below 0.01 and convergence was assumed. The first 25% of the trees were discarded as burnin using the sumt command. The convergence of the two simultaneous runs, and stationary state of each parameter, were assessed by examining Trace plots and histograms in Tracer v1.6 (Rambaut et al. 2014). Runs were terminated when the effective sample sizes (ESS) of all parameters ≥ 200.
The most likely tree in the ML analysis, and the 50% majority-rule consensus of the sampled trees from the BI analysis, were visualized using FigTree v1.4.3 (Rambaut 2009). Nodes having bootstrap support (BS) of ≥ 95 and posterior probabilities (PP) of ≥ 0.95 were considered to be well-supported (Huelsenbeck and Ronquist 2001;Wilcox et al. 2002). Uncorrected pairwise genetic distances were calculated using MEGA v7.0.26 (Kumar et al. 2016).

Morphological analyses
The ANOVA found statistically significant differences in morphometric characters of the Lan Saka samples and four congeners in the siamensis group (p < 0.05) for all fifteen variables, as did the Tukey's HSD pairwise (p < 0.05; Table 2).
The PCA of five species of Cnemaspis showed large morphometric differences on a scatter plot of the first four components with eigenvalues greater than 1.0 (Fig. 1A). These four components accounted for 85.40% of the total variance (Table 3). The first principal component (PC1) accounted for 33.88% of the most of variance and loaded heavily on the head proportions (interorbital distance, eye-nostril distance and eye-snout distance) and the shape of tail (tail width). The second principal component (PC2) accounted for 25.70% and mostly loaded for the body proportion (axillar-groin length) and the head proportions (internarial distance, head length, eye-ear distance and ear length). The third principal component (PC3) accounted for 17.10% and loaded heavily on the head proportions (head width and head depth) and forearm length whereas the fourth (PC4) accounted for 8.72% and loaded heavily on the head proportions (head width, ear length and head length) and the body proportions (axillagroin length and tibia length). Factor loadings for each component are provided in Table 3. The ordination of the first two components showed separation between the Lan Saka samples and four congeners in the siamensis group. The PC2 axis showed separation between C. adangrawi, C. omari, and C. siamensis from C. chanardi and the Lan Saka samples. The biplot analysis showed that the Lan Saka samples overlapped slightly with C. chanardi. The DAPC (97.70% of cumulative variance) discriminated among groups and supported distinct clusters that corresponded to five Cnemaspis species (Fig. 1B).

Molecular analyses
The aligned dataset contained 1,251 characters of 69 individuals of Cnemaspis and two individuals of the outgroup species. The standard deviation of split frequencies among the two simultaneous BI runs was 0.001646. The ESS values were greater than or equal to 2,944 for all parameters. The maximum likelihood value of the best ML tree was lnL = -54,716.041. The most likely ML tree and the 50% majority rule consensus tree from the BI analysis resulted in trees with similar topologies (Fig. 2). The Lan Saka samples represented a well-supported clade (100 BS, 1.0 PP) within the siamensis group and the sister taxon of a clade containing C. adangrawi, C. chanardi, C. phangngaensis, C. omari, and C. roticanai (Fig. 2), although relationships within that sister clade were not resolved (Fig. 2). Sequence divergences (uncorrected p-distance for ND2) ranged from 0.00-0.40% within the Lan Saka samples and 15.53-28.09% among the Lan Saka samples and other species in the siamensis group (Table 4).
Body slender, elongate (AG/SVL 0.43); small, keeled, dorsal scales equal in size throughout body intermixed with several large, keeled, multicarinate tubercles; 19 paravertebral tubercles linearly arranged; tubercles present on lower flanks; tubercles extend from occiput to tail; five small, subconical spine-like tubercles on flanks; dorsal scales raised and keeled; pectoral and abdominal scales keeled, round, flat to concave, slightly larger than dorsal and not larger posteriorly; ventral scales of brachia smooth, raised and juxtaposed; six separated pore-bearing precloacal scales with rounded pores; precloacal depression absent; femoral pores absent.
Coloration in life (Fig. 4). Dorsal ground color of head light brown, top of head bearing small, diffuse, faint black and yellowish markings; dark postorbital stripes faint, extending to nape; large, round, whitish marking on nape; single light-yellowish prescapular crescent on the shoulder, located dorsoanteriorly of forelimb insertion; dorsal ground color of body, limbs and tail light brown with black irregular blotches; ground color of ventral surfaces grayish-white intermixed with yellowish blotches; ventral pattern sexually dimorphic, anterior gular, abdominal, and caudal regions yellowish in males; two dark blotches on nape form a bipartite pattern; light sage vertebral blotches extending from the nape to tail; flanks with irregular, incomplete brown to yellowish blotches becoming smaller posteriorly; tubercles on the whole body white or yellow; widely separated, white or yellow tubercles occur on flanks; subconical spine-like yellowish tubercles on flanks; limbs beige with dark brown mottling; tail faintly marked with dark brown.
Coloration in preservative (Figs 5, 7, 8). Color pattern similar to that in life with some fading of markings. Dorsal ground color of head, body, limbs and tail darker brown than in life, with indistinct, irregular markings. Yellow coloration in gular, pectoral, abdominal regions, flanks, and tail faded to light-yellow and creamy-white.
Variation. Most paratypes approximate the holotype in general aspects of morphology (Figs 6-8), with most differences found in the degree of vertebral blotches. All adult female paratypes lack the yellowish coloration in the gular, abdominal, and caudal regions. ZMKU R 00821-00825, ZMKU R 00829-00831, THNHM 28694 (nine adult males) and ZMKU R 00826, ZMKU R 00832 and ZMKU R 00833 (three adult females) have regenerated tails of uniform tan coloration. ZMKU R 00821, 00824, 00827 (three adult males) and ZMKU R 00832, 00835 (two adult females) have lighter dorsal markings that appear more as transverse bands than as paravertebral blotches. THNHM 28696 (adult female) has a broken tail. Differences in meristic and morphometric characters within the type series are presented in Table 6.
The male holotype was found during the night (1943 h) perched head down on a vertical surface in a crevice of a granitic rock boulder near a stream. A female paratype (ZMKU R 00832) was found with the male holotype, separated by only a distance of approximately 10 cm.
Paratypes that were found during the day were in shaded areas, crevices of boulders, rock walls and on boulder outcrops near streams. Paratypes found at night were in shaded surfaces of the boulders, within deep crevices, or perched on vegetation near a rocky stream. Three gravid females (ZMKU R 00832-00834) contained one or two eggs during January 2019. Some juveniles (SVL < 30 mm; not collected) were found in rock cracks and perched on a rock near a stream on 25 January 2019.
Cnemaspis lineatubercularis sp. nov. appears to be a diurnal species in that observed specimens during daytime were active and fast-moving when disturbed, but those at night were inactive, slow-moving or asleep on dry granitic rocks and vegetations. At night, Cyrtodactylus lekaguli and Gehyra mutilata were found in syntopy with the new species on a rock wall and vegetation near a stream. A summary of ecological parameters of activity periods, elevation (lowland < 600 m), microhabitat preference and presence or absence of ocelli (eyespots) of Cnemaspis in Thailand is shown in Table 7.
Etymology. The specific epithet lineatubercularis is taken from linea (Lat. for line) and tubercularis (Lat. for having tubercles), in reference to the new species having paravertebral tubercles linearly arranged.

Discussion
The complex geological history of Thailand created a large number of granitic rocky outcrop ecosystems in southern Thailand (Charusiri 1993;Cobbing et al. 2011). This ecosystem supports high levels of species endemism and species diversity of gekkonid lizards, especially species in the genus Cnemaspis (see figure 5 in Grismer et al. 2014). The findings of this study provide new data from a poorly studied area in Nakhon Si Thammarat Province, southern Thailand. The results suggest that additional unexplored regions may still harbor unrecognized species of Cnemaspis in Thailand.
A decade ago, only four species of Cnemaspis were known from Thailand, including C. biocellata, C. chanthaburiensis, C. kumpoli, and C. siamensis (Smith 1925;Taylor 1963;Bauer and Das 1998;Grismer et al. 2008a, b).  described seven new species of Cnemaspis (C. chanardi, C. huaseesom, C. kamolnorranathi, C. narathiwatensis, C. niyomwanae, C. punctatonuchalis, and C. vandeventeri) from Thailand. Previously, Grismer et al. (2014) described a new species C. omari from Perlis, Malaysia, that is also distributed in adjacent Satun Province, Thailand. Wood et al. (2017) described three additional new species of Cnemaspis (C. lineogularis, C. phangngaensis, and C. thachanaensis) from southern Thailand. Most recently, two new insular species of Cnemaspis (C. adangrawi and C. tarutaoensis) were described from Tarutao, Adang and Rawi islands of southern Thailand (Ampai et al. 2019). The discovery and description of C. lineatubercularis sp. nov. brings the total number of Thai Cnemaspis species to 18, representing one-third (33%) of the 60 named species in Southeast Asia.
The new species is known only from the type locality and likely has a narrow geographic distribution. It is expected to be found in other nearby granitic rocky streams in Kam Lon Subdistrict, Lan Saka District, Nakhon Si Thammarat Province. However, additional surveys for this species are needed to clarify the geographic range of the new species. Our findings agree with those of Grismer et al. (2014) that most species of Cnemaspis in Thailand are diurnal, granite-associated, lowland species that lack ocelli ( Table 7). Further research and additional field surveys in unexplored regions of lowland forest in southern Thailand are needed to better understand the taxonomy, ecology, distribution, biogeography, and conservation of Cnemaspis in the region.