Field work was conducted between 2009 and 2018 in Laos and Vietnam. Specimens were euthanized with ethyl acetate, fixed in approximately 85% ethanol, then transferred to 70% ethanol for permanent storage. Specimens were subsequently deposited in the collections of the Institute of Ecology and Biological Resources (IEBR), Vietnam Academy of Science and Technology, Hanoi, Vietnam; the Vietnam National Museum of Nature (VNMN), Hanoi, Vietnam; the Vietnam National University of Forestry (VNUF), Hanoi, Vietnam; the National University of Laos (NUOL), Laos; and the Zoological Research Museum Alexander Koenig (ZFMK), Bonn, Germany.

Main morphological characters were rechecked: Measurements were taken with a digital caliper to the nearest 0.1 mm. Abbreviations are as follows: snout-vent length (SVL, from tip of snout to anterior margin of cloaca); tail length (TaL, from posterior margin of cloaca to tip of tail).

Scale counts were taken using stereo microscopes (Leica S6E, Keyence VHX-500F): ventral scales in longitudinal rows at midbody (V) counted transversely across the center of the abdomen from one ventrolateral fold to the other; dorsal tubercle rows (DTR) counted transversely across the center of the dorsum from one ventrolateral fold to the other; supralabials (SL) and infralabials (IL) counted from the first labial scale to the corner of mouth; enlarged femoral scales (EFS); femoral pores (FP); precloacal pores (PP) or the total number of femoral pores and precloacal pores (i.e. the contiguous rows of femoral and precloacal scales bearing pores combined as a single meristic character referred to as the femoroprecloacal pores); number of subdigital lamellae on the fourth finger (LD4) and number of subdigital lamellae on the fourth toe (LT4) counted from the base of the first phalanx to the claw.

Most described taxa of the genus Cyrtodactylus in Laos and Vietnam, except for C. buchardi, C. raglai, and C. septimontium were included in the study. In addition, samples of the species from different localities were sequenced to determine their distribution range. In total, 84 new samples from 26 provinces were incorporated (Suppl. material 1: Table S1). The tissue samples of muscle, liver or tail tissue was preserved separately in 70% ethanol. The mitochondrial DNA cytochrome c oxidase subunit I (COI) gene was selected as the markers have been widely used in previous studies and for some geographic populations, only COI data were available (Nguyen et al. 2015, 2017; Luu et al. 2016a; Brennan et al. 2017; Ngo et al. 2017). In several cases, where comparative data for COI were not available, the mitochondrial gene NADH dehydrogenase subunit 2 (ND2) was generated for specimens under consideration. In addition, we obtained 90 sequences of the mitochondrial COI for the ingroup taxa and one outgroup species, Hemidactylus frenatus, from GenBank (Wood et al. 2012).

Total genomic DNA was extracted using DNeasy Blood and Tissue Kit (Qiagen, Germany), following protocols by the manufacturer’s instructions. PCR was performed using HotStar Taq Mastermix (Qiagen, Germany) to amplify approximately 657 bp fragment of the mitochondrial gene COI and approximately 1200 bp fragment of the mitochondrial gene ND2. We used two primer pairs for PCR with VF1d (5’–TTCTCAACCAACCACAARGAYATYGG-3’), VR1d (5’– TAGACTTCTGGGTGGCCRAARAAYCA–3’) (Ivanova et al. 2006) for generating the COI and MetF1 (5’– AAGCTTTCGGGCCCATACC–3’), COIR1 (5’– AGRGTGCCAATGTCTTTGTGRTT–3’) (Macey et al. 1997 for the ND2. The PCR volume consisted of 21 µl (2 µl each primer, 5 µl water, 10 µl of Taq mastermix and 1 – 2 µl of DNA depending on the quality of DNA in the final extraction solution). The reactions were performed at 95 °C for 15 min, followed by 35 cycles of 30 s at 95 °C, 45 s at 48 °C, and 60 s at 72 °C with a final elongation step of 6 min at 72 °C. A negative and positive control was used for every DNA extraction and PCR reactions. PCR products were visualized using electrophoresis through a 1% agarose gel, marker 1 kb, 1X TBE and stained with ethidium bromide and photographed under UV light. Successful amplifications were purified using GeneJet PCR Purification Kit (ThermoFisher Scientific, Lithuania). Cleaned PCR products was sent to 1^st Base (Malaysia) for sequencing.

Newly generated sequences were checked by eye using Sequencher v5.4 (Gene Codes Corp, Ann Arbor, MI, USA), aligned by ClustalX v2.1 (Thompson et al. 1997) with default setting. The data were then analyzed using Bayesian inference (BI) as implemented in MrBayes v3.1.2 (Ronquist et al. 2012) and maximum likelihood (ML) analysis using IQ-TREE v.1.6.7.1 (Nguyen et al. 2015) with a single molecular model and 10,000 ultrafast bootstrap (UFBP) replications. For BI, the analysis was conducted with a random starting tree and run for 1 × 10⁷ generations. Four Markov chains, one cold and three heated (utilizing default heating values), were sampled every 1000 generations. Log likelihood scores of sample points were plotted against generation time to detect stationarity of the Markov chains. The burn-in value was set to 59 in the BI analysis, as –lnL scores reached stationarity after 59,000 generations in both runs. The optimal model, GTR+I+G, for nucleotide evolution was set to BI and ML analysis as selected by jModeltest v2.1.4 (Darriba et al. 2012). Nodal support was evaluated using UFBP as estimated in IQ-TREE v1.6.7.1 and posterior probability (PP) in MrBayes v3.2. UFBP and PP ≥ 95% are regarded as strong support for a clade (Ronquist et al. 2012; Nguyen et al. 2015). Uncorrected pairwise distance (p) was calculated in PAUP*v4.0b10 (Swofford 2001).

We successfully sequenced a fragment of the COI gene for 90 samples and ND2 gene for four samples. The final concatenated matrix consisted of 216 terminals, including 90 from this study, 126 from previous studies, including one outgroup, Hemidactylus frenatus, following Wood et al. (2012). Both BI and ML analyses based on a total of 657 aligned characters with no gaps and using a single model of molecular evolution produced very similar topologies (Fig. 1). All species groups were generally recovered with strong support values from both analyses, except for the C. irregularis group, which was corroborated only by the BI analysis (Fig. 1). The results show that Cyrtodactylus species in Laos and Vietnam fall into six species groups, namely C. angularis, C. brevipalmatus, C. chauquangensis, C. condorensis, C. intermedius, and C. irregularis. While Vietnam harbors members of all groups, Lao species mostly belong to the two karst-dwelling C. angularis and C. chauquangensis groups (Grismer et al. 2021a, b). Two other groups contain one taxon each from Laos, including C. cf. cryptus of the C. irregularis group (samples of C. buchardi and C. pseudoquadrivirgatus from Laos not included in the analysis) and C. cf. ngati of the C. brevipalmatus group. Two remaining mostly insular groups, the C. condorensis and C. intermedius groups, do not have any species from Laos (Fig. 1, Suppl. material 1: Table S1).

Figure 1. 

Bayesian cladogram based on 657 bp of the partial COI gene. Numbers above branches are Bayesian posterior probabilities and ultrafast bootstrap values of ML analysis, respectively. ★ = new sequences used in the phylogenetic analyses.

The main difference between this and previous studies is that the Cyrtodactylus angularis group was rendered paraphyletic. Nguyen et al. (2015) and Ngo et al. (2017) used the same COI region and confirmed that the monophyly of the group is significantly supported in the BI analysis, while its monophyletic relationship based on 1474 base pairs of the ND2 gene received perfect statistical values from both BI and ML phylogenetic estimates (Grismer et al. 2021a). All other species groups were recovered with strong support values from both analyses, except for the C. irregularis group, which was corroborated only by the BI analysis but received perfect statistical support in both analysis in Grismer et al. (2021a).

In the Cyrtodactylus angularis group, almost all species are well defined and supported by both analyses. According to our tree, the group contains 19 known species and one undescribed taxon (Cyrtodactylus sp. 1) in both countries. There are only three species that have notable genetic sub-structuring, i.e., C. darevskii, C. multiporus, and C. pageli, and samples from genetically distinct populations are labeled as cf. (Fig. 1). Another population from Khammouane Province, Laos with three representative samples, KM2012.52, KM2012.54–1, and KM2012.54–2, is clearly differentiated from other species and likely forms a new species. The former two samples were incorrectly assigned to C. lomyenensis in previous studies (Nazarov et al. 2014, Schneider et al. 2014). It is marked as Cyrtodactylus sp. 1 on the tree (Fig. 1). The highest pairwise genetic divergence between species of the Cyrtodactylus angularis group is 24.04% based on the fragment of COI gene (Suppl. material 2: Table S2). The lowest pairwise genetic divergence between species that are not morphologically well differentiated (SVL, LT4, LD4, Infralabials, EFS, color pattern of dorsum, and enlarged subcaudals specifications, see Tables 2, 4) is 2.44% (C. rufford and C. lomyenensis). The lowest divergence between species that are morphologically distinct is 4.57% (C. jarujini and C. thathomensis).

Members of the Cyrtodactylus brevipalmatus group, recently discovered in Vietnam (Le et al. 2021), are present in both countries. The genetic divergence between the Vietnamese and Lao populations of C. ngati is 2.13% and they are ~ 3.81–4.41% separated from C. cf. interdigitalis from Thailand. This group contains the lowest number of taxa (Fig. 1, Suppl. material 1: Table S1, Suppl. material 6: Table S6). The C. chauquangensis group is also distributed in both countries with 16 described species and one undescribed form, Cyrtodactylus sp. 7, from Vietnam. Samples of purported C. bichnganae are all labeled as cf. because genetic sequences of the true C. bichnganae have not been made available (see Pham et al. 2019). The second and third smallest groups, the C. intermedius and C. condorensis groups, have three and four species, respectively. All of the species occur in Vietnam and mostly inhabit offshore islands in the southern part of the country.

The largest group, the Cyrtodactylus irregularis group, consists of more than 30 species with at least five undescribed forms, Cyrtodactylus sp. 2 – Cyrtodactylus sp. 6. The highest and lowest pairwise genetic divergence that exists between species of the C. irregularis group are 21.41% and 0.74% (Suppl. material 3: Table S3). In addition, C. cf. pseudoquadrivirgatus is recovered in three distinct places of the tree, one close to C. taynguyenensis, one embedded within a clade comprising Cyrtodactylus sp. 4, Cyrtodactylus sp. 5, C. culaochamensis, and the other was recovered as a sister taxon to C. cryptus. Finally, C. dati and C. thuongae are not highly divergent in terms of pairwise genetic distance and only separated by 0.74% (Table 3, Suppl. material 3: Table S3).

The Cyrtodactylus condorensis group is composed of four well defined species with pairwise divergences of 5.48 – 18.05% (Suppl. material 4: Table S4). The Cyrtodactylus chauquangensis group is composed of 17 species which have remained stable in many analyses of phylogenetic relationships. The genetic divergences between the members of the group are ~ 3.81–19.54% (Suppl. material 5: Table S5).

Our results based on a fragment of the mitochondrial gene COI show that the lineage containing Cyrtodactylus thuongae with the holotype and paratype from Tay Ninh Province and C. dati from the Lam Dong, Binh Phuoc populations is divided into two sub-lineages. However, the PP value is insignificant (Fig. 1). These populations are separated by no greater than 1% in sequence divergence using COI data (Table 3). In addition, genetic distances between the holotype of C. dati (UNS 0543) and the holotype and paratype of C. thuongae (IEBR A.2013.23, IEBR A.2013.25) are 2.21% and 1.85%, respectively, based on 951 bp of a fragment of the mitochondrial ND2 gene. At the time of the latter species’ description, molecular data for both species were not available for comparison. Instead, C. thuongae was described on the basis of its highly variable morphology (Table 1). Morphological examination of specimens of C. thuongae and published data of C. dati show that morphologically they do not have a high level of distinction (size, number of ventral scales, infralabials, EFS, FP, LT4, enlarged subcaudals, color pattern of dorsum. The only differences between the two species are: DTR: 17 or 18 in C. thuongae vs. 20–22 in C. dati, SL: 7–9 in C. thuongae vs. 10–12 in C. dati and 0 or 1 pitted scales in males in C. thuongae vs. five or six in C. dati.

Much the same is true for C. lomyenensis from Khammouan (paratype – UNS0527) and C. rufford also from Khammouan (holotype – VFU R.2015.14). Genetic divergence between the two species is less than 2.44% based on a fragment of mitochondrial gene COI. According to our morphological examinations between specimens of C. rufford and published data of C. lomyenensis, C. rufford differs from C. lomyenensis by having fewer ventral scale rows (27–29 vs. 35 or 36), fewer supralabials (10–12 vs. 13 or 14), and more femoral and precloacal pores in males (42 or 43 vs. 39 or 40) (Table 2). The differences are quite small, except for the ventral scale row. In addition, C. rufford is similar to C. lomyenensis in size and coloration: head dorsum yellowish with irregular brown blotches, dorsal pattern with transverse bands, rings on original tail with dark brown transversal bands wider than light brown spaces, and median row of enlarged subcaudal scales. In addition, genetic distance between the paratype of C. lomyenensis (UNS0527) and holotype of C. rufford (VFU R.2015.14) is 0.21% based on 413 bp of a fragment of the mitochondrial ND2 gene.