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, or limited 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 recovered monophyletic mitochondrial lineages of individuals (populations) used to develop initial species-level hypotheses—the grouping stage of Hillis (2019). Discrete color pattern data and univariate and multivariate analyses of morphological data were then used to search for characters and morphospatial patterns consistent with the tree-designated species-level hypotheses—the construction of boundaries representing the hypothesis-testing step of Hillis (2019)—thus providing independent diagnoses to complement the molecular analyses. It is important to note, that delimiting species (phylogeny) and diagnosing species (taxonomy) are independent but overlapping operations that should not be conflated (Frost and Hillis 1990; Frost and Kluge 1994; Hillis 2019).

These boundaries were cross-checked using a Generalized Mixed Yule Coalescent (GMYC) approach (Pons et al. 2006) employed using the splitSelect package version 1.0.3 in R (Christidis et al. 2021) and a Bayesian Poisson Tree Process for species delimitation (bPTP), thus providing an independent framework to complement the empirically based thresholds of the morphological analyses. Both approaches are a method for delimiting species from single-locus gene trees with low population samples (Lin et al. 2018) by detecting genetic clustering beyond the expected levels of a null hypothesis which infers that all individuals of a population form a genetically, interacting nexus. In clades where effective population sizes are relatively low and divergence times among the populations are relatively high, the single-threshold version of the model (such as that used herein) for the GMYC outperforms the multi-threshold version (Fujisawa and Barrenclough 2013). For the bPTP, Markov Chain Monte Carlo (MCMC) was run for 10,000 generations on the bPTP web server (Zhang et al. 2013) and checked for convergence. Both models rely on the prediction that independent evolution leads to the appearance of distinct genetic clusters, separated by relatively longer internal branches (Barraclough et al. 2003; Acinas et al. 2004). Such groups therefore, diverge into discrete units of genetic variation that are recovered with surveys of higher clades.

The morphological data set of Grismer et al. (2021c) was augmented with the paratypes of C. interdigitalis (THNHM 20226–29), seven specimens of C. cf. rukhadeva from Khaeng Krachan National Park, Phetchaburi Province, Thailand (THNHM 01807, 03251–54, 24838), two specimens of sp. 13 from Thung Yai Naresuan Wildlife Sanctuary, Tak Province (THNHM 00104) and Ban Saphan Lao, Kanchanaburi Province (THNHM 27821), one specimen of sp. 14 from Khlong Nakha Wildlife Sanctuary, Ranong Province (THNHM 01667), and two specimens of C. brevipalmatus (THNHM 10670, 14112) from Khao Nan National Park and Khao Luang National Park, Nakhon Si Thammarat Province (Fig. 1). Grismer et al.’s (2021c) molecular data set was augmented with one specimen of C. interdigitalis (YC00952) from the type locality and one specimen of C. cf. interdigitalis from Khammouane Province, Laos (VNUF R.2014.50). Methods for DNA extraction, sequencing and editing followed Grismer et al. (2021c) and resulted in a 1399 base pair segment of ND2 and adjacent tRNAs. All material examined is listed in Table 1 along with GenBank accession numbers for the new and published genetic material.

Morphological data included both meristic and morphometric characters. To reduce the degree of researcher bias, data were taken using the protocol of Le et al. (2021) and where possible, double checked by LLG using high resolution digital photographs and/or the actual specimens. All data were taken on the left side of the body (when possible) and measured to the nearest 0.1 mm using dial calipers under a dissecting microscope.

Morphometric data 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 proximal 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 proximal 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 margin 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 posterior to the eyes; 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 dorsomedial-most edges of the bony orbits; internarial distance (IN), measured between the external 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 at the corner of the mouth or posterior to the eye, to the rostral scale; infralabial scales (IL), counted from termination of enlarged scales at the corner of the mouth to the mental scale; number of paravertebral tubercles (PVT) between the limb insertions counted in a straight line immediately left of the vertebral column; number of longitudinal rows of body tubercles (LRT) counted transversely across the body midway between the limb insertions from one ventrolateral body fold to the other; number of longitudinal rows of ventral scales (VS) counted transversely across the abdomen midway between limb insertions from one ventrolateral fold to the other; number of transverse rows of ventral scales (VSM) counted along the midline of the body from the postmentals to just anterior to the cloacal opening, stopping where the scales become granular; 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; 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; total number of subdigital lamellae (TL4T) beneath the fourth toe (i.e. TL4E + TL4U = TL4T); number of expanded subdigital lamellae on the fourth finger proximal to the digital inflection (FL4E) counted the same way as with TL4E; number of small generally unmodified subdigital lamellae distal to the digital inflection on the fourth finger (FL4U) counted the same way as with TL4U; total number of subdigital lamellae (FL4T) beneath the fourth toe (i.e. FL4E + FL4U = FL4T); total number of enlarged femoral scales (FS) from each thigh combined as a single metric; number of enlarged precloacal scales (PCS); number of precloacal pores (PP) in males; the number of femoral pores (FP) in males; and the number of dark body bands (BB) between the dark band on the nape and the hind limb insertions on the body. Categorical characters evaluated were the presence or absence of tubercles on the flanks (FKT; Fig. 2), single, enlarged, unmodified, medial subcaudal scales (SC2; Fig. 3A), enlarged, posteriorly emarginated, medial subcaudals bearing a median furrow (SC3; Fig. 3B), or no enlarged medial subcaudals (SC1; Fig. 3C); large or small dorsolateral caudal tubercles (DCT) forming a continuous or discontinuous fringe (VL1; Fig. 4), ventrolateral caudal fringe narrow or wide (VL2; Figs 3 and 4), and the cross-section of the tail round or square (TLcross; Fig. 4). The raw morphological data for all characters and specimens are presented in Table 2.

Figure 2. 

Tubercles on the left flank of Cyrtodactylus sp. 13 (THNHM 00104) from Thung Yai Naresuan Wildlife Sanctuary, Tak Province, Thailand.

Figure 3. 

Subcaudal scale and ventrolateral caudal fringe morphology. A Cyrtodactylus cf. rukhadeva (THNHM 03251) from Kaeng Krachan National Park, Phetchburi Province, Thailand. B Cyrtodactylus interdigitalis (THNHM 20226) from Tham Yai Nam Nao, Nam Nao National Park, Phetchabun Province, Thailand. C Cyrtodactylus sp. 13 (THNHM 00104) from Thung Yai Naresuan Wildlife Sanctuary, Tak Province, Thailand.

Figure 4. 

Dorsolateral caudal tubercle and ventrolateral caudal fringe morphology. A Cyrtodactylus cf. brevipalmatus (LSUHC 11788) from Langkawi Island, Kedah State, Peninsular Malaysia. B Cyrtodactylus elok (LSUHC 8238) from Negeri Sembilan State, Peninsular Malaysia.

Ingroup samples consisted of 16 individuals of the brevipalmatus group representing five nominal species (Table 1). All species of the pulchellus group were used to root the tree following Grismer et al. (2021b). The protein-coding region and the flanking tRNAs were aligned using the MAFTT v7.017 (Katoh et al. 2002) plugin under the default settings in Geneious 2019.0.4 (https://www.geneious.com). Maximum likelihood (ML) and Bayesian Inference (BI) analyses were used to estimate phylogenetic trees. Best-fit models of molecular evolution determined in IQ-TREE (Nguyen et al. 2015) using the Bayesian information criterion (BIC) implemented in ModelFinder (Kalyaanamoorthy et al. 2017) indicated that K2P+I was the best-fit model of evolution for the tRNAs and TN+F was the best model of evolution for codon positions 1 and 2 and HKY+F for codon position 3. 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 1,500 generations and the first 10% of each run 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 using Tracer v1.6 (Rambaut et al. 2018) to ensure effective sample sizes (ESS) for all parameters were above 200. Bayesian posterior probabilities (BPP) of 0.95 and above and ultrafast bootstrap support values (UFB) of 95 and above were considered an indication of strong nodal support (Huelsenbeck et al. 2001; Minh et al. 2013). Uncorrected pairwise sequence divergences 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.

All statistical analyses were conducted using R Core Team (2018). Morphometric characters used in statistical analyses were SVL, AG, HumL, ForL, FemL, TibL, HL, HW, HD, ED, EE, ES, EN, IO, EL, and IN. Tail metrics were not used due to a high degree incomplete sampling (i.e. regenerated, broken, or missing). To remove potential effects of allometry (sec. Chan and Grismer 2022), 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 potential intra- with interspecific variation (Reist 1986; McCoy et al. 2006). The juvenile C. ngati (HNUE-00112) was removed from the data so as not to skew the normalization results. All data were scaled to their standard deviation to ensure they were analyzed on the basis of correlation and not covariance. Meristic characters analyzed were SL, IL, PVT, LRT, VS, VSM, TL4E, TL4U, TL4T, FL4E, FL4U, FL4T, FS, PCS, and BB. Precloacal and femoral pores were omitted from the multivariate analyses due to their absence in females. Categorical characters analyzed were DCT, VFL1, VFL2, TLcross, SC1, SC2, and SC3.

Small sample sizes (n=1 or 2) for some of the species/ populations precluded standard statistical analyses for relevant groups closely related to C. interdigitalis. However, morphospatial clustering and positioning among the species/populations was analyzed using multiple factor analysis (MFA) on a concatenated data set comprised of 15 meristic characters, 16 normalized morphometric characters, and seven categorical characters (Table 2). The MFA was implemented using the mfa function in the R package FactorMineR (Husson et al. 2017) and visualized using the Factoextra package (Kassambara and Mundt 2017). MFA is a global, unsupervised multivariate analysis that incorporates qualitative and quantitative data (Pagès 2015), making it possible to analyze different data types simultaneously in a nearly total evidence environment. Grismer and Chan (in progress) have empirically demonstrated that this approach outperforms less data-rich multivariate methods in differentiating monophyletic lineages in multivariate space. In an MFA, each individual is described by a different set of variables (i.e. characters) that are structured into different data groups in the data frame—in this case, quantitative data (i.e. scale counts and normalized morphometrics) and categorical data (i.e. scale, tubercle, and caudal morphology). In the first phase of the analysis, separate multivariate analyses are carried out for each data group—principal component analysis (PCA) for quantitative data and multiple correspondence analysis (MCA) for categorical data. The first eigenvalue from each of these analyses is retained and used to weight the data groups in the second phase of the analysis—a PCA of the weighted data. Standardizing the data in this manner prevents one data type from overleveraging another. In other words, the standardization of the data in the first phase prevents data types with the most number of characters from outweighing data types with fewer characters in the second phase. This way, the contribution of each data type to the overall variation in the data set is scaled to define the morphospatial distance between individuals as well as calculating each data type’s contribution to the variation in the overall analysis (Pagès 2015; Kassambara and Mundt 2017).

In order to further examine the morphometric dissimilarity among C. ngati and closely related specimens from Laos and Thailand that have been referred to as C. interdigitalis (see below), a PCA and discriminant analysis of principal components (DAPC) of the 16 normalized morphometric characters was employed. 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 undesignated clusters that may or may not align with the putative species boundaries delimited by phylogenetic analyses and when possible, defined by univariate analyses. It is important to understand that clusters of conspecific individuals delimited a priori in the phylogeny are not pre-defined in the analysis but simply color-coded in the scatter plot in order to observe their positions and morphospatial relationships. DAPC from the adegent package 2.1.5 in R (Jombart 2021) is a supervised analysis (i.e. groups are specified a priori in the analysis) that relies on scaled data calculated from a 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 and then choosing an appropriate number of linear discriminants on which to display those data on a scatterplot (Jombart and Collins 2015).

A two-sample Student t-test on data meeting the assumptions of normality and homogeneity of variance, and a Welch’s two-sample t-test on data violating those assumptions, were run for each morphometric character to test for the presence or absence of significantly different mean values (p<0.05). The ranges, means, medians, and 50% quartiles of each character were visualized using violin plots with embedded boxplots. All analyses were performed in R [v3.4.3].

Following the MFA and PCA, a non-parametric permutation multivariate analysis of variance (PERMANOVA) from the vegan package 2.5-3 in R (Oksanen et al. 2020) was used to determine if the centroid locations of each species/population were statistically different from one another (Skalski et al. 2018). The analyses were based on the calculation of a Euclidean (dis)similarity matrix using 5,000 permutations. A pairwise post hoc test calculates the differences between all combinations of population pairs, generating a Bonferroni-adjusted p value and a pseudo-F ratio (F statistic). A p < 0.05 is considered significant and larger F statistic values indicate more pronounced group separation. A rejection of the null hypothesis (i.e. centroid positions and/or the spread of the data points (i.e. clusters) are no different from random) signifies a statistically significant difference between species/populations.

Mismatch distributions using the PopGenome package version 2.7.5 in R (Pfeifer et al. (2020) based on the infinite-sites model (Kimura 1969, 1971) were performed on the genetic data in order to compare observed base pair differences to those of a simulated distribution under the sudden population expansion model. Populations at demographic equilibrium or in decline should present a multimodal distribution of pairwise differences, while populations that have experienced a sudden demographic expansion should display a unimodal distribution (Slatkin and Hudson 1991; Rogers and Harpending 1992). To complement the mismatch distributions, an independent neutrality test (Tajima 1989) was employed using the pegas package version 1.1 in R (Paradis et al. 2021) to test for possible selection and population expansion. Tajima’s D statistic was calculated under the infinite sites model using 1000 simulated samples, which assumes average heterozygosity for a pair of randomly chosen alleles and is compared with the expected number of sites segregating in each sample.

A Mantel randomization test from the adegent package 2.1.5 in R (Jombart 2021) and a partial distance-based redundancy analysis (dbRDA) from the vegan package 2.5-3 in R (Oksanen et al. 2020) were used to ascertain if there is a correlation between (dis)similarity matrices of genetic and geographic distances at nodes within the interdigitalis clade and C. ngati (see below) that could be a function of isolation by distance (IBD; Slatkin 1993)—just one of a number of possibilities. It should be noted, however, that Legendre et al. (2015) demonstrated the potentially inappropriate use of the Mantel test to analyze spatial data (but see Diniz-Filho et al. 2013). Uncorrected pair-wise genetic distances calculated in MEGA 7 were converted into Euclidean distances using the dist() function. A principal coordinate analysis (PCoA) was performed on the GPS coordinates of all samples in order to transform them into a Gower pairwise (dis)similarity matrix. To test for IBD in both analyses, genetic distances were treated as response variables and the transformed GPS data were the independent variables. The mantel.randtest() function was implemented and run for 10,000 permutations and the observed r statistic was visualized on a histogram of permutations. Observing the location of the r statistic on the histogram, allows one to assess the likelihood of the observed correlation arising by chance (Jombart 2021). An r statistic falling in the middle of the histogram indicates no correlation, whereas an r statistic skewed to right would indicate a correlation. A p_adj<0.05 is considered evidence of significant correlation.

To complement the Mantel test, a dbRDA analysis was performed on the two matrices using the capescale() function from the vegan package. Statistical significance using 999 permutations was assessed by calculating and R², R²_adj, and F test p-values. The data matrices were regressed in a linear regression analysis using the lm() function and plotted on a heat map using a 2-dimensional kernal density estimation (kde2d function) from the MASS package 7.3-54 (Ripley et al. (2021). The R² and p_adj-values were compared to the results of the Mantel and dbRDA tests.

The ML and BI analyses recovered trees with well-supported (UFB 100/BPP 1.00) identical topologies (Fig. 5). Cyrtodactylus brevipalmatus and C. cf. brevipalmatus of the Thai-Malay Peninsula were recovered as the well-supported sister lineage to the remainder of the brevipalmatus group unlike in Grismer et al. (2021c), where it was recovered as the poorly supported sister species of C. elok. Here, C. elok is recovered as the poorly supported (61/--) sister species to the remainder of the group. Other than these rearrangements, the phylogenies recovered here are essentially the same as those in Grismer et al. (2021c). The relevant difference being the addition of C. interdigitalis (YC00952) from the type locality being recovered as the well-supported sister species of ZMMU R-16492 (sp. 11 in Grismer et al. 2021c) from Phu Hin Rong Kla National Park, Phitsanulok Province, Thailand. Given that ZMMU R-16492 and YC00952 are sister species, the former is referred to here as sp. 11 also given that they are separated by a straight-line distance of approximately 70 km that encompasses a 22 km wide uninhabitable lowland river basin. Additionally, they share an uncorrected pairwise sequence divergence between them of 4.6%. A complete morphological comparison between them is not possible given that the preserved specimen of ZMMU R-16492 lacks a tail that contains a number of diagnostic characters (see Table 2). Together, these specimens and the remainder of the group are referred to here as the interdigitalis clade (Fig. 5).

Grismer et al. (2021c) referred to specimens FMNH 265806 from Thailand and FMNH 255454 and 270492 and NCSM 79472 and 80100 from Laos as C. cf. interdigitalis owing to a lack any sort of data from C. interdigitalis from the type locality. The molecular data clearly show that FMNH 265806, 255454, and 270492 do not form an exclusive clade with C. interdigitalis but do so with two of the paratypes of C. ngati with an uncorrected pairwise sequence divergence among them of only 0.9–1.4%. We therefore refer to these specimens as C. ngati. NCSM 79472 from Ban Pha Liep, Houay Liep Stream, Xaignabouli Province and NCSM 80100 from Houay Wan Stream, a tributary of Nam Pha River, Laos, however, fall outside C. ngati but are sequentially related to it, and are referred to here as C. cf. ngati 1 and 2, respectively (Fig. 5). The taxonomy of these populations awaits the acquisition of additional material.

The GMYC species delimitation independently recovered FMNH 265806, 255454, 270492 and VNUF R.2014.50 and C. ngati (IBER 4829 and VNUF R.2020.12) as conspecific. However, the likelihood ratio test was insignificant (p=0.10000) which is not surprising given the high percentage of singletons in the data set (Talavera et al. 2013). The bPTP analysis also recovered these specimens as conspecific with a 0.84 posterior probability.

Figure 5. 

A Maximum likelihood consensus tree based on 1399 base pairs of ND2 and flanking tRNAs with UFBS and BBP nodal support. B Histogram of uncorrected pairwise sequence divergence of specimens in the interdigitalis clade from Cyrtodactylus interdigitalis YC00952 from the type locality of Tham Yai Nam Nao, Nam Nao National Park, Phetchabun Province, Thailand.

The MFA of the concatenated data sets corroborate the phylogenetic analyses, in part, in that C. brevipalmatus (combined with C. cf. brevipalmatus from here on out), C. elok, C. interdigitalis, and C. rukhadeva (including the type series and specimens from Kaeng Krachan National Park; see below), and sp. 9. were recovered as distinct well-separated species along the combined axes of Dim-1 and Dim-2 (Fig. 6). The extreme morphological distinction of C. elok forced the close proximity of the other species/populations along a longer Dim-1, thus obscuring the magnitude of the relative differences among the others (Fig. 6A). Although removing C. elok from the data set did not alter their positional relationships of the others, it recovered a clearer separation among them (Fig. 6B). The MFAs recovered close morphological similarity among C. cf. ngati 1 and 2 and the Thai and Lao specimens of C. ngati along with a wide separation of those populations from C. ngati from the type locality in Vietnam (Figs 1 and 6A, B). In the data set lacking C. elok (Fig. 6B), meristic data contributed to 40% of the variation along Dim-1 followed by the categorical and morphometric data (Fig. 6C). For Dim-2 two, morphometric data contributed 40% of the variation followed by meristic and categorical data. So as not to inflate the p-adjusted values of the data set lacking C. elok, populations whose species status was not in question (based on their phylogenetic relationships) and represented by only one or two specimens, were deleted prior to the implementation of the PERMANOVA test. These included C. cf. ngati 1 and 2, sp. 9, and sp. 12. The PERMANOVA test returned a significant difference in the centroid placement among all possible combinations of population pairs (p<0.05) and a more highly significant difference between C. rukhadeva and all other populations (p-adjusted<0.05) (Table 3).

Figure 6. 

A. MFA of all species/populations of the brevipalmatus group. B MFA of all species/populations of the brevipalmatus group minus Cyrtodactylus elok. C Histograms of the percent contribution of each data type to the variation in dimensions 1 (Dim-1) and dimensions 2 (Dim-2) based on 6B. The dashed red line represents the average of percent contribution.

The MFA recovered the holotype of C. rukhadeva from Khao Laem Mountain, Suang Phueng District, Ratchaburi Province, Thailand and the paratype from Hoop Phai Tong, Suang Phueng District 7.7 km to the east, as overlapping along Dim-2 and in close morphospatial proximity to each other and a series of seven specimens (THNHM 01807, 03251–54, 24622, 24838) of C. cf. rukhadeva (Grismer et al. 2021c) from Kaeng Krachan National Park, Phetchaburi Province, ~ 83 km to the south with which they overlapped along Dim-1 (Figs 1, 6B). The two taxa are identical in all seven categorical characters and overlap in all 19 meristic characters except for the number of femoral pores in males (Table 2). The holotype of C. rukhadeva (ZMMU R-16851) has 17 femoral pores whereas the number of pores in five males from Kaeng Krachan (THNHM 01807, 03251–52, 03254, 24622) ranges from 11–14. However, given C. rukhadeva is known from only a single male and there is a range of at least four pores in the Kaeng Krachan population, it is likely that additional specimens from both populations will result in an overlap in this character as with the 26 other characters (Table 2). Based on these data and the MFA, we hypothesize that these three populations are likely conspecific and we recognize the Kaeng Krachan population as C. cf. rukhadeva. Additional evidence will be required to test this hypothesis.

The MFA recovered the specimen from the Khlong Nakha Wildlife Sanctuary, Ranong Province, Thailand (sp. 12, THNHM 01667) as widely separated from C. brevipalmatus (Fig. 6) that occurs ~130 km to the southeast across the Isthmus of Kra—a well-known biogeographical barrier (Fig. 1). THNHM 01667 differs from the five specimens of C. brevipalmatus examined here (AUP-00573, LSUHC 11788, THNHM 10670, 14112, USMHC 2555) in having large as opposed to small dorsolateral caudal tubercles and ventrolateral caudal fringes—characters that are invariant in the other species (Table 2). The next geographically closest species is C. cf. rukhadeva ~412 km to the north and C. elok ~524 km to the southeast (Fig. 1). Based on these data, THNHM 01667 is hypothesized to be a distinct species. The acquisition of additional specimens and genetic data (in progress) will test this hypothesis.

The MFA recovered significant morphological distinction between the type series of C. ngati from Vietnam and the geographically distant Thai and Lao populations of C. ngati plus C. cf. ngati 1 and 2 (Figs 1 and 6). Because body proportions of karst-adapted Cyrtodactylus vary greatly relative to their closely related non-karst adapted species (see Grismer et al. 2016, Nielson and Oliver 2017, Quah et al. 2019, Wood et al. 2020, Kaatz et al. 2021), a PCA and DAPC examining only morphometric data from C. cf. ngati 1 and 2 and C. ngati returned a wide separation of the types of C. ngati from the others (Fig. 7A). Principal component (PC) 1 accounted for 52.7% of the variation and loaded heavily for AG, FemL, HL, HW, HD, ED, and ES (Table 4). PC2 accounted for an additional 12.9% of the variation and loaded heavily for HumL and IO. The DAPC also returned a wide separation of the two morphological groups based on the retention of the first five (PCs), which accounted for 97.3% of the variation (Fig. 7B). Student and Welch’s two sample t-tests demonstrated that C. ngati from the type locality is significantly smaller in having a shorter SVL and AG; has a smaller and flatter head (HL and HW, and HD, respectively) with a smaller eye (ED), a shorter snout (ES), a smaller ear opening (EL); and shorter hind limbs (FemurL and TibL) but not forelimbs (Figs 7C; Table 5). The PERMANOVA test returned a highly significant p-adjusted value (0.0087) for the separation of their centroids.

Figure 7. 

A PCA of specimens in the interdigitalis clade. B DAPC of specimens in the interdigitalis clade. C Violin plots overlain with box plots showing the range, frequency, mean (white dot), and 50% quartile (black rectangle) of the size-adjusted morphometric characters. * denotes characters with significantly different mean values.

The mismatch analyses returned multimodal distributions for the interdigitalis clade and C. ngati, indicating they are not undergoing range expansion (Fig. 8). This was supported by significantly negative Tajima’s D statistics: D = –2.1128; p_norm= 6.75e-05 and p_beta = 0.00 for the interdigitalis clade and D = –6.2531; p_norm = 4.071e-10 and p_beta = 0.00 for C. ngati.

Based on eigenvalue decomposition, the loadings of the first two principal components from the PCoA for the geographic data were retained as spatial variables. The very low calculated Mantel r statistic and insignificant p-value (r=-0.1578; p=0.7438) between the genetic and transformed geographic distance (dis)similarity matrices of individuals in the interdigitalis clade fell midway within the range of the observed permutations, indicating there is no correlation between genetic and geographic distances (i.e. potentially no IBD; Fig. 9A). This was consistent with the regression analysis which returned an R² of -0.0038 (i.e. no correlation; p=0.3581) and a negatively sloping regression line (Fig. 9B). Additionally, no correlation was recovered in the more exclusive analysis within C. ngati (Mantel r=0.2420; p=0.2599; Fig. 9C). The regression analysis here was also consistent with the Mantel test (R²=0.0591; p=0.5005; Fig. 9D) indicating that only 5.9% of the genetic variation may be due to geographic location. Results of the dbRDA analyses mirrored those of the Mantel tests in that no correlations were recovered between the different data sets: p=0.396 and R²_adj=0.190 for the interdigitalis clade and p=0.233 and R²_adj=0.608 for C. ngati.

Figure 8. 

Mismatch distribution analyses and calculated Tajima’s D statistics. A The Cyrtodactylus interdigitalis clade. B Cyrtodactylus ngati.

Figure 9. 

A Observed simulations of the Mantel test of the individuals in the interdigitalis clade based on 10,000 permutations. B Regression analysis of the geographic and genetic distance matrices and heat map of clusters of the individuals in the interdigitalis clade. C Observed simulations of the Mantel test of the individuals of Cyrtodactylus ngati based on 10,000 permutations. D Regression analysis of the geographic and genetic distance matrices and heat map of clusters of individuals of C. ngati. Dashed red line is the regression line. Vertical line with diamond tip refers to the value of the r statistic.