Specimens were collected in Suan Phueng District, Ratchaburi Province, Thailand by Parinya Pawangkhanant, Platon V. Yushchenko, Mali Naiduangchan, Chatmongkon Suwannapoom, and Nikolay A. Poyarkov during several field surveys in 2019 (Fig. 2). Geographic coordinates and elevations were obtained using a Garmin GPSMAP 60CSx and recorded in WGS 84 datum. Specimens were collected by hand, anaesthetized and euthanized with 20% benzocaine solution, femoral muscle or liver tissue samples removed, and then fixed in 10% buffered formalin before preserving in 70% ethanol. The tissue samples were stored in 95% ethanol. Specimens and tissues were subsequently deposited in the herpetological collections of the School of Agriculture and Natural Resources, University of Phayao (AUP, Phayao, Thailand) and of the Zoological Museum of Moscow University (ZMMU, Moscow, Russia).

Specimen collection and animal use protocols were approved by the Institutional Ethical Committee of Animal Experimentation of the University of Phayao, Phayao, Thailand (certificate number UP-AE61-01-04-022, issued to Chatmongkon Suwannapoom) and were strictly complacent with the ethical conditions of the Thailand Animal Welfare Act. Field work, including collection of animals in the field and specimen exportation, was authorized by the Institute of Animals for Scientific Purpose Development (IAD), Bangkok, Thailand (permit numbers U1-01205-2558 and UP-AE59-01-04-0022, issued to Chatmongkon Suwannapoom).

Morphological data included both meristic and morphometric characters. Measurements were taken on the left side of the body when possible to the nearest 0.1 mm using dial calipers under dissecting microscope following Murdoch et al. (2019). 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, original or partially regenerated; tail width (TW), taken at the base of the tail immediately posterior to the postcloacal swelling; humeral length (HumL), taken from the medial end of the humerus at its insertion point in the glenoid fossa to the distal margin of the elbow while flexed 90°; forearm length (ForL), taken on the ventral surface from the posterior margin of the elbow while flexed 90° to the inflection of the flexed wrist; femur length (FemL), taken from the medial end of the femur at its insertion point in the acetabulum to the distal margin of the knee while flexed 90°; tibia length (TibL), taken on the ventral surface from the posterior surface of the knee while flexed 90° to the base of the heel; 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; 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; eye to snout distance or snout length (ES), measured from anteriormost 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 dorsal-most edges of the bony orbits; internarial distance (IN), measured between the nares across the rostrum; and ear length (EL), greatest oblique length across the auditory meatus.

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 transverse rows of ventral scales (VSM), counted along the midline of the body from the postmentals to the cloacal opening; the number of expanded subdigital lamellae on the fourth toe proximal to the digital inflection (TL4E), counted from the base of the first phalanx where it contacts the body of the foot to the largest scale on the digital inflection—the large contiguous scales on the palmar and plantar surfaces were not counted; the number of small, generally unmodified subdigital lamellae distal to the digital inflection on the fourth toe (TL4U), counted from the digital inflection to the claw including the claw sheath (see Murdoch et al. 2019: Fig. 2); the total number of subdigital lamellae (TL4T) beneath the fourth toe (i.e. TL4E + TL4U = TL4T); the number of expanded subdigital lamellae on the fourth finger proximal to the digital inflection (FL4E) counted from the base of the first phalanx where it contacts the body of the foot to the largest scale on the digital inflection—the large contiguous scales on the palmar and plantar surfaces were not counted; the number of small, generally unmodified subdigital lamellae distal to the digital inflection on the fourth finger (FL4U) counted from the digital inflection to the claw including the claw sheath (see Murdoch et al. 2019: Fig. 2); the total number of subdigital lamellae (FL4T) beneath the fourth toe (i.e. FL4E + FL4U = FL4T); the total number of enlarged femoral scales (FS) from each thigh combined as a single metric; the number of enlarged precloacal scales (PCS); the number of precloacal pores in (PP) in males; the number of femoral pores in (FP) in males; the number of rows of post-precloacal scales (PPS) on the midline between the enlarged precloacal scales and the vent; 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 the presence or absence of paravertebral tubercles, enlarged femoral scales, small tubercles on the forelimbs and flanks, paried or single enlarged subcaudal scales; large or small dorsolateral and ventrolateral caudal fringe scales (i.e. fringes), and the cross-section of the tail round or more square in shape.

Small sample sizes (n=1–3) from some of the populations/species precluded meaningful statistical analyses. However, the morphospatial clustering among the species for the meristic, morphometric, and concatenated datasets was visualized using principal component analysis (PCA) and discriminant analysis of principal components (DAPC) from the ADEGENET package in R (Jombart et al. 2010). We analyzed each data type (i.e. meristic and morphometric) separately in order to visualize the contribution of each to the variation in the group and then we concatenated the data types to visualize their overall contribution to the variation of the group. Previous studies have demonstrated that concatenated data sets generally outperform single data sets in capturing the overall morphological variation among closely related species (Grismer et al. 2018, 2020a, 2021).

Morphometric characters analyzed were SVL, AG, HumL, ForL, FemL, TibL, HL, HW, HD, ED, EE, EN, IO, EL, and IN. To remove potential effects of allometry in the morphometric dataset, size was normalized using the following equation: X_adj=log(X)-β[log(SVL)-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, 1983; Turan 1999; Lleonart et al. 2000, accessible in the R package GroupStruct (available at https://github.com/chankinonn/GroupStruct). The morphometrics of each species were normalized separately and then concatenated so as not to conflate intra- with interspecific variation (Reist 1986). All data were scaled to their standard deviation to insure they were analyzed on the basis of correlation and not covariance.

Meristic characters analyzed were SL, IL, LRT, VS, VSM, TL4E, TL4U. TL4T, FL4E, FL4U, FL4T, PCS, PPS, and BB. Paravertebral tubercles (PVT), enlarged femoral scales (FS) were not included due to their absence in Cyrtodactylus elok. Precloacal and femoral pores were also omitted due to their absence in females. For corroboration of the PCA, a discriminant analysis of principal components (DAPC) was performed on the same data sets. DAPC is a supervised analysis that 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 90–99% of the variation in the data set (Jombart and Collins 2015) as determined from a scree plot generated as part of the analysis. The ranges, means, medians, and 50% quartiles were visualized using boxplots for the meristic data and violin plots with embedded boxplots for the mensual data. All analyses were performed in R [v3.4.3].

For the molecular phylogenetic analyses, we extracted the total genomic DNA from ethanol-preserved femoral muscle tissue of six specimens of the new Thai populations using standard phenol-chloroform-proteinase K extraction procedures with consequent isopropanol precipitation, for a final concentration of about 1 mg/ml (protocols followed Hillis et al. (1996) and Sambrook and David (2001)). We visualized the isolated total genomic DNA in agarose electrophoresis in the presence of ethidium bromide. We measured the concentration of total DNA in 1 μl using NanoDrop 2000 (Thermo Scientific), and consequently adjusted to ca. 100 ng DNA/μL.

We amplified mtDNA fragments of ND2 and its flanking tRNAs using a double-stranded Polymerase Chain Reaction (PCR) to obtain a 1397 base pairs under the following conditions: 1.0 μl genomic DNA (10–30 μg), 1.0 μl light strand primer (concentration 10 μM), 1.0 μl heavy strand primer (concentration 10 μM), 15 μl Master Mix 2x (CWBIO, China), and 12 μl ultra-pure H2O. PCR reactions were executed on Bio-Rad T100™ gradient thermocycler under the following conditions: initial denaturation at 94 °C for 5 min, followed by a second denaturation at 94 °C for 60 s, annealing at 58 °C for 60 s, followed by a cycle extension at 72 °C for 60 s, for 35 cycles with the final extension step 72 °C for 10 min. For amplification of the full ND2 gene with parts of adjacent tRNAs we used the NADH2-Metf6 (5’-AAGCTTTCGGGCCCATACC-3’) and CO1H (5’-AGRGTGCCAATGTCTTTGTGRTT-3’) primers following Macey et al (1997).

PCR products were loaded onto 1.5% agarose gels in the presence of ethidium bromide and visualized in agarose electrophoresis. When distinct bands were produced, we purified the PCR products using 2 μl of a 1:4 dilution of ExoSapIt (Amersham) per 5 μl of PCR product prior to cycle sequencing. A 10 μl sequencing reaction included 2 μL of template, 2.5 μl of sequencing buffer, 0.8 μl of 10 pMol primer, 0.4 μl of BigDye Terminator version 3.1 Sequencing Standard (Applied Biosystems) and 4.2 μl of water. The cyclesequencing used 35 cycles of 10 sec at 96°C, 10 sec at 50°C and 4 min at 60°C. We purified the cycle sequencing products by ethanol precipitation. We carried out sequence data collection and visualization on an ABI 3730xl Automated Sequencer (Applied Biosystems).

Ingroup samples consisted of 14 individuals of the brevipalmatus group representing three nominal species including a sample from one of the new populations in Suan Phueng District plus five individuals of Cyrtodactylus cf. interdigitalis from Thailand and Laos. Four species of the linnwayensis group were used to root the tree following Grismer et al. (2021b). Sequence data for some specimens were acquired from GenBank and newly generated sequences were deposited in GenBank (Table 1).

Maximum Likelihood (ML) and Bayesian Inference (BI) were used to estimate phylogenetic trees. Best-fit models of evolution determined in IQ-TREE (Nguyen et al. 2015) using the Bayesian information criterion (BIC) implemented in ModelFinder (Kalyaanamoorthy et al. 2017) indicated that TNe+I+G4 was the best-fit model of evolution for the tRNAs and HKY+F+G4, TIM3+F+G4, and TPM3u+F+I+G4 were the best models of evolution for codon positions 1, 2, and 3, respectively. The ML analysis was performed using the IQ-TREE webserver (Trifinopoulos et al. 2016) with 1000 bootstrap pseudoreplicates using the ultrafast bootstrap (UFB) analysis (Minh et al. 2013; Hoang et al. 2018). The BI analysis was performed on CIPRES Science Gateway (Miller et al. 2010) using MrBayes v3.2.4 (Ronquist et al. 2012). Two independent runs were performed using Metropolis-coupled Markov Chain Monte Carlo (MCMCMC), each with four chains: three hot and one cold. The MCMCMC chains were run for 15,000,000 generations with the cold chain sampled every 1500 generations and the first 10% of each run being discarded as burn-in. The posterior distribution of trees from each run was summarized using the sumt function in MrBayes v3.2.4 (Ronquist et al. 2012). Stationarity was checked with Tracer v1.7 (Rambaut et al. 2018) to be sure effective sample sizes (ESS) for all parameters were well above 200. We considered Bayesian posterior probabilities (BPP) of 0.95 and above and ultrafast bootstrap support values (UFB) of 95 and above as an indication of strong nodal support (Huelsenbeck et al. 2001; Minh et al. 2013). Uncorrected pairwise sequence divergences (p-distance) were calculated in MEGA 7 (Kumar et al. 2016) using the complete deletion option to remove gaps and missing data from the alignment prior to analysis.

The ML and BI analyses recovered trees with identical topologies (Fig. 2A). Cyrtodactylus elok and C. brevipalmatus of the Thai-Malay Peninsula were recovered as poorly supported sister species (BPP 0.57/UFB 59) and the strongly supported (1.00/100) sister lineage to a strongly supported clade (1.00/100) containing the remaining taxa. Within that clade, the individual of the new population from Suan Phueng District was recovered as the sister lineage to a monophyletic group bearing the phylogenetic sequence of sp. 9, sp. 10, sp. 11, and terminating with the strongly supported (1.00/100) sister species C. ngati and C. cf. interdigitalis. Support values for these nodes in are shown in Figure 2A. The phylogenetic analyses strongly supported the recognition of the Suan Phueng population as a new species, which bears a 15.4–22.1% uncorrected pairwise sequence divergence from the other species in the brevipalmatus group (Table 2). The phylogenetic analyses also supported the separate species status of sp. 9 and sp. 10 as previously recovered by Chomdej et al. (2021).

The phylogeny indicates species recognition for ZMMU R-16492 from Phu Hin Rong Kla National Park, Phitsanulok Province, Thailand (Fig. 1; referred to here as sp. 11) in that it is recovered as the sister species to the strongly supported (1.00/99) lineage comprised of C. ngati and C. cf. interdigitalis (Fig. 2). Although the locality of ZMMU R-16492 at 1147 m in elevation is only 45 km to the northwest of the type locality of C. interdigitalis at 700 m in Nam Nao National Park, Petchabun Province (Ulber 1993), these upland sites are separated by a 20 km wide low-lying river basin of 136 m in elevation at its lowest point with no intervening hilly terrain (Fig. 1). Similarly, C. cf. interdigitalis FMNH 265806 from the Phu Luang Wildlife Sanctuary occurs at 850 m in elevation and is 65 km north of the type locality of C. interdigitalis but separated from it by an approximately15 km wide low-lying river basin at approximately 200 m at its lowest point. ZMMU R-16492 lies 50 km to the southwest of FMNH 265806 and it too is separated by a low-lying river basins. Given the close geographic proximity of these three topographically separated upland localities and the fact that ZMMU R-16492 and FMNH 265806 are not each other’s closest relatives, makes it premature to consider either as conspecific with C. interdigitalis in the absence of genetic data from the type locality or at least morphological data from the type series. Both options are currently being explored.

Although the monophyly of Cyrtodactylus cf. interdigitalis and C. ngati is strongly supported, the monophyly of each of them is not. Although it is difficult to conceive that C. ngati would not be monophyletic given that all specimens were found on the same karst formation approximately 200 km north of the nearest population of C. cf. interdigitalis (Fig. 1), nodal support for its monophyly is only 0.62/82. That of C. cf. interdigitalis ranging from northern Laos to central Thailand is even less—0.51/52. If these nodes were collapsed on the basis of low support values, it would be phylogenetically conceivable to consider C. ngati and C. cf. interdigitalis as conspecific. However, these two species are trenchantly different in morphology (see below) and again, the absence of any data from the type locality of C. interdigitalis, makes any taxonomic proposal premature.

The PCA of the concatenated meristic and adjusted morphometric data sets corroborated the phylogenetic analyses in that the Suan Phueng individuals are not embedded within the cluster of any other species along the combined ordination of principal component (PC) 1 and PC2 (Fig. 2B). PC1 accounted for 36.8% of the variation and loaded most heavily for AG, FemL, HW, HD, EL and FS. PC2 accounted for 15.3% of the variation and loaded most heavily for FL4T and BB (Table 3). These results were mirrored in the DAPC of the concatenated data (Fig. 2C). The PCA of the morphometric data also recovered separation of the Suan Phueng individuals along both axes but very close to C. cf. interdigitalis. In the DAPC, however, the Suan Phueng individuals were within the 66% confidence ellipsoid of C. cf. interdigitalis (Fig. 3). PC1 accounted for 50.2% of the variation and loaded most heavily for ES, HD, AG, FemL, and HL (Fig. 4) and PC2 accounted for an additional 16.2% of the variation, loading most heavily for IO, IN, and EN. The meristic data also recovered separation of the Suan Phueng populations along the combined axes in the PCA and DAPC but again, closest to C. cf. interdigitslis (Fig. 5). PC1 accounted for 35.0% of the variation and loaded most heavily for AG, FemL, HW, HD, EL and FS (Fig. 6) and PC2 accounted for 20.1% of the variation, loading most heavily for FL4E, TL4T, and BB.

Based on the phylogenetic relationships, PCA, DAPC, and several discrete morphological differences between the Suan Phueng individuals and all other species of the brevipalmatus group (see comparisons below and Table 4), we hypothesize that they comprise a discretely diagnosable lineage that is not reticulating with any other lineage and as such, should be accorded species status.