During October 2017 and June 2019, field surveys were conducted at five localities in southern Thailand, including the type locality of C. phuketensis (Fig. 1; Table 1). Specimens were investigated and captured by hand during the night (1900–2300). Liver or muscle tissues were individually preserved in 95% ethyl alcohol and stored at -20 °C for molecular analysis. Specimens were fixed in 10% formalin and later transferred to 70% ethyl alcohol. Preserved specimens were deposited in the herpetological collections of the Zoological Museum, Kasetsart University, Thailand (ZMKU). Additional specimens were also examined from the Princess Maha Chakri Sirindhorn Natural History Museum (PSU), Prince of Songkhla University, Thailand; Thailand Natural History Museum (THNHM), Thailand; La Sierra University Herpetological Collection (LSUHC), La Sierra University, Riverside, California, USA; and the Zoological Reference Collection (ZRC) of Lee Kong Chian Natural History Museum at National University of Singapore, Singapore.

Figure 1. 

Map illustrating the known geographic distribution of Cyrtodactylus macrotuberculatus and C. phuketensis. Yellow star: the type locality of C. macrotuberculatus at Gunung Raya, Pulau Langkawi, Kedah, Malaysia. Green star: the type locality of C. phuketensis at Thalang District, Phuket Island, Phuket Province. Yellow circles: C. macrotuberculatus samples used in this study. Green circle: C. phuketensis samples used in this study. Yellow squares: the distribution of C. macrotuberculatus taken from Grismer et al. (2012), and Quah et al. (2019). The samples used correspond to those in Table 1.

Total genomic DNA was extracted from 95% ethanol-preserved muscle or liver tissue using a NucleoSpin Tissue Kit (Macherey-Nagel GmbH & Co. KG, Germany). Mitochondrial NADH dehydrogenase subunit 2 (ND2) gene and flanking tRNAs were amplified via double-stand Polymerase Chain Reaction (PCR) using primers L4437a (tRNAmet: 5’ AAGCTTTCGGGCCCATACC 3’) and H5934 (COI: 5’ AGRGTGCCAATGTCTTTGTGRTT 3’) (Macey et al. 1997). PCR amplification occurred with an initial denaturation at 94 °C for 4 min, followed by 35 cycles of denaturation at 94 °C for 30 sec, annealing at 48–52 °C for 30 sec, extension at 72 °C for 1 min 30 sec, and final extension at 72 °C for 7 min. Amplification products were purified using NucleoSpin Gel and PCR Clean-Up kit (Macherey-Nagel GmbH & Co. KG, Germany) and visualized on 1.0% agarose gel electrophoresis. Purified PCR products were sequenced in both directions using amplifying primers on an ABI 3730XL DNA Sequencer (Applied Biosystems, CA, USA). Sequences were manually edited and aligned in Geneious R11 (Biomatters, Ltd, Auckland, New Zealand). The ND2 nucleotide sequences were translated to amino acid for the protein-coding region and to ensure the lack of stop codons. All sequences were deposited in GenBank under the accession numbers MW809294 to MW809309 (Table 1).

Phylogenetic relationships were inferred using two model based approaches, Bayesian Inference (BI) and Maximum Likelihood (ML). Outgroup species used to root the tree were Hemidactylus frenatus, Agamura persica, Tropiocolotes steudneri, C. elok, C. intermedius, C. interdigitalis, and Cyrtodactylus sp. based on Wood et al. (2012). The best-fit nucleotide substitution model for each of the three codon partitions and tRNAs was selected under the Bayesian Information Criterion (BIC) in PartitionFinder2 on XSEDE (Lanfear et al. 2016) using CIPRES (Cyberinfrastructure for Phylogenetic Research; Miller et al. 2010). BI analysis was executed in MrBayes 3.2.6 on XSEDE (Ronquist et al. 2012) using CIPRES with the TRN+I+G for the 1^st and 2^nd codon position, and TIM+I+G for the 3^rd codon position and tRNAs. Four chains (three hot and one cold) were run for 10,000,000 generations and sampled every 1,000 generations using Markov chain Monte Carlo (MCMC). To build a consensus tree, we discarded the first 25% of each run as burn-in and assessed stationarity by plotting log-likelihood score in Tracer ver. 1.7.1. (Rambaut et al. 2018). The ML analysis was performed on the web server W-IQ-TREE (Trifinopoulos et al. 2016) with 1,000 bootstrap replicates using ultrafast bootstrap approximation (Minh et al. 2013). Nodes having Bayesian posterior probabilities (BPP) of ≥ 0.95 and ultrafast bootstrap support (UFB) of ≥ 95 were considered to be strongly supported (Huelsenbeck and Ronquist 2001; Wilcox et al. 2002; Minh et al. 2013). Uncorrected pairwise sequence divergence was calculated to assess within and among species differences using the default settings in MEGA X 10.0.5 (Kumar et al. 2018).

Morphological and meristic characters were modified from the previous studies of Grismer and Ahmad (2008) and Grismer et al. (2012). Measurements were taken with digital calipers to the nearest 0.1 mm for the following sixteen characters:

SVL snout-vent length, taken from the tip of snout to the vent;

TW tail width, taken at the base of the tail immediately posterior to the postcloacal swelling;

TL tail length, taken from vent to the tip of the tail, original or regenerated;

FL forearm length, taken from the posterior margin of the elbow while flexed 90° to the inflection of the flexed wrist;

TBL tibia length, taken from the posterior surface of the knee while flexed 90° to the base of the heel;

AG axilla to groin length, 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;

HL head length, the distance from the posterior margin of the retroarticular process of the lower jaw to the tip of the snout;

HW head width, measured at the angle of the jaws;

HD head depth, the maximum height of head from the occiput to the throat;

ED eye diameter, the greatest horizontal diameter of the eye-ball;

EE eye to ear distance, measured from the anterior edge of the ear opening to the posterior edge of the eye-ball;

ES eye to snout distance, measured from anterior most margin of the eye-ball to the tip of snout;

EN eye to nostril distance, measured from the anterior margin of the eye-ball to the posterior margin of the external nares;

IO inter orbital distance, measured between the anterior edges of the orbit;

EL ear length, the greatest vertical distance of the ear opening;

IN internarial distance, measured between the nares across the rostrum.

Additional scale counts and non-meristic characters evaluated were the number of supralabial and infralabial scales counted from the largest scale immediately posterior to the dorsal inflection of the posterior portion of the upper jaw to the rostral and mental scales, respectively; the number of paravertebral tubercles between limb insertions counted in a straight line immediately left of the vertebral column; the number of longitudinal rows of body tubercles counted transversely across the center of the dorsum from one ventrolateral fold to the other; the number of longitudinal rows of ventral scales counted transversely across the center of the abdomen from one ventrolateral fold to the other; the presence or absence of tubercles on the ventral surface of the forearm; the presence or absence of tubercles in the gular region, throat, and lateral margins of the abdomen; the number of subdigital lamellae beneath the fourth toe counted from the base of the first phalanx to the claw; the total number of precloacal and femoral pores (i.e., the contiguous rows of femoral and precloacal scales bearing pores combined as a single meristic referred to as the femoroprecloacal pores); the presence or absence of a precloacal depression or groove; the degree of body tuberculation, weak tuberculation referring to dorsal body tubercles that are low and rounded whereas prominent tuberculation refers to tubercles that are raised and keeled; the width of the dark body bands relative to the width of the interspace between the bands; number of dark caudal bands on the original tail; the white caudal bands of adults immaculate or infused with dark pigment; and whether or not the posterior portion of the original tail in hatchlings and juveniles less than 50 mm SVL was white or whitish and faintly banded or boldly banded.

All statistical analyses were performed using the base statistical software in RStudio v. 1.2.1335 (RStudio Team 2018). To remove potential effects of allometry, mensural characters were scaled to SVL using the following allometric equation: X_adj = X-β(SVL-SVL_mean), where X_adj = adjusted value; X = measured value; β = unstandardized regression coefficient for each OTU; SVL = measured snout-vent length; SVL_mean = overall average SVL of each OTU (Thorpe 1975, 1983; Turan 1999; Lleonart et al. 2000). Male and female measurements were analyzed separately to remove potential effects of sexual dimorphism. For morphological analyses, TL (tail length) was excluded due to their different conditions (e.g., original, broken, and regenerated). Importantly, the following type material and topotypic specimens were included in the analysis: C. macrotuberculatus (holotype and three paratypes) and C. phuketensis (holotype, paratype and three topotypes). Prior to the morphological analyses, individuals were assigned on the basis of molecular data except C. phuketensis based on its distribution into three groups (= species): C. macrotuberculatus, C. phuketensis, and C. pulchellus.

Principal component analysis (PCA) was implemented in the R package FactoMineR (Lê et al. 2008) to discover or reduce dimensionality of the original character variables in order to recover characters bearing the highest degree of variation among groups. Fifteen scaled morphometric and seven meristic characters (scalations) were concatenated and used for the PCA analyses separately by sex. For females, femoroprecloacal pore counts were excluded from the PCA due to their presence in only males.

For univariate analyses, all transformed mensural characters were tested for normality using the Shapiro-Wilk Test. Equality of variances was tested using F-tests. Morphological differences of both males and females between C. macrotuberculatus and C. phuketensis were examined using a t-test (for normally distributed and equal variance data), Welch’s t-test (for unequal variance data) and Mann-Whitney U test (for non-normally distributed data) at a significant level of 95%.

The aligned matrix contained 1,453 bp of ND2 gene and its flanking tRNAs for 101 samples of the C. pulchellus complex including outgroups (Table 1). The standard deviation of split frequencies among the two Bayesian runs was 0.003263, and the Estimated Sample Size (ESS) of all parameters were ≥ 200. The BI and ML analyses generated similar topologies and strong nodal support for most clades, and only the ML tree is shown (Fig. 2). According to phylogenetic analyses, C. phuketensis is nested within C. macrotuberculatus with strong support (1.00 BPP, 100 UFB), thus rendering C. macrotuberculatus paraphyletic. Cyrtodactylus macrotuberculatus (including C. phuketensis) was recovered as sister lineage to a clade containing C. pulchellus and C. evanquahi. Uncorrected pairwise sequence divergence (p-distance) between C. phuketensis and this sister lineage was higher than 8.45% and within the C. phuketensis and C. macrotuberculatus clade, it ranged from 1.45–4.20% (mean 2.63%; Table 2). The p-distance within species ranged from 0.00–0.36% (mean 0.14%) for C. phuketensis and 0.00–4.38% (mean 2.48%) for C. macrotuberculatus.

Figure 2. 

Reconstructed phylogenetic relationships of the Cyrtodactylus pulchellus complex based on 1,453 bp of ND2 and flanking tRNAs. The phylogenetic tree is from the Maximum Likelihood analysis with Bayesian posterior probabilities (BPP) and ultrafast bootstrap support values (UFB), respectively. Black circles represent nodes supported by BPP and UFB of 1.0 and 100. Samples in bold are new sequence from this study.

A total of 45 preserved specimens from three species groups (C. macrotuberculatus = 29, C. phuketensis = 10, and C. pulchellus = 6) were used for principal component analysis (Table 3). The PCA of males showed complete overlap between C. macrotuberculatus and C. phuketensis, and they were completely separated from C. pulchellus along the first two principal components (Fig. 3A). The first three principal components of males accounted for 53.17% of the variation. The first principal component (PC1) accounted for 25.78% of the variation and was most heavily loaded on HL_adj, ES_adj, EN_adj and ventral scales; the PC2 accounted for 17.56% of the variation and was most heavily loaded on TBL_adj, IO_adj, supralabial and infralabial scales; and the PC3 accounted for 9.83% of the variation and was loaded most heavily on longitudinal tubercles (Table 3).

Figure 3. 

Plots for the first two principal components of morphological characters from A males, and B females resulting from the principal component analyses of Cyrtodactylus macrotuberculatus (yellow circles), C. phuketensis (green squares) and C. pulchellus (red diamonds). The letters in the scatter plots refer to holotype (= H), paratype (= P) and topotype (= T).

Along the first two PC plots, the PCA of females revealed complete overlap between C. macrotuberculatus and C. phuketensis, which were distinctly separated from C. pulchellus (Fig. 3B). The first three principal components of females accounted for 52.99% of the variation. The first principal component (PC1) accounted for 23.14% of the variation and was most heavily loaded on TW_adj, AG_adj, ES_adj, EN_adj, ventral scales and number of the 4^th toe lamellae; the PC2 accounted for 16.59% of the variation and was most heavily loaded on HL_adj, ED_adj, supralabial and infralabial scales; the PC3 accounted for 13.26% of the variation and was loaded most heavily on TBL_adj and EL_adj (Table 3).

Summary univariate statistics of morphological characters of adult males and females are shown in Table 4. In adult males, C. macrotuberculatus (N = 18) and C. phuketensis (N = 6) were not significantly different in most morphological characters (t-tests and Mann-Whitney U tests, p = 0.1765–0.9523) except only IO_adj (t-test, p = 0.0256). In adult females, C. macrotuberculatus (N = 11) and C. phuketensis (N = 4) were not significantly different in twelve morphological characters (t-tests and Welch’s t-test, p = 0.2325–0.9626) whereas only three characters were significantly different which are TBL_adj (t-test, p = 0.0495), AG_adj (Mann-Whitney U tests, p = 0.0176) and IN_adj (Welch’s t-test, p = 0.0129).

Sumontha et al. (2012) distinguished C. phuketensis from C. macrotuberculatus by using the number of subdigital lamellae on the fourth toe, number of dark body bands, and the presence of a precloacal groove in females. Based on examination of the type material and newly collected specimens of C. phuketensis from Phuket Island, these diagnostic characters overlap with those characters of C. macrotuberculatus. In this study, C. phuketensis has 19–21 total subdigital lamellae on the fourth toe (vs. 19–23 in C. macrotuberculatus); three or four dark body bands (vs. three or four in C. macrotuberculatus), and no precloacal groove in females (also absent in C. macrotuberculatus; Table 5).

The phylogenetic analyses recovered C. phuketensis as being nested within the C. macrotuberculatus and bearing a low genetic divergence (mean 2.63%) which was similar to that within C. macrotuberculatus populations (mean 2.48%). In concordance, the statistical analyses of meristic and mensural characters of C. phuketensis widely overlap with those of C. macrotuberculatus. Based on these data, we propose that C. phuketensis from Phuket Island, Phuket Province is a junior synonym of C. macrotuberculatus which can be recognized as follows.