Research Article
Print
Research Article
The taxonomy and phylogeny of the Cyrtodactylus brevipalmatus group (Squamata: Gekkonidae) with emphasis on C. interdigitalis and C. ngati
expand article infoL. Lee Grismer§, Attapol Rujirawan|, Siriporn Yodthong, Bryan L. Stuart#, Minh Duc Le¤«, Dzung Trung Le», Yodchaiy Chuaynkern˄, Perry L. Wood, Jr.˅, Anchalee Aowphol|
‡ La Sierra University, Riverside, United States of America
§ San Diego Natural History Museum, San Diego, United States of America
| Kasetsart University, Bangkok, Thailand
¶ Kasestart Universoty, Bangkok, Thailand
# North Carolina Museum of Natural Sciences, Raleigh, United States of America
¤ Vietnam National University, Hanoi, Vietnam
« American Museum of Natural History, New York, United States of America
» Ministry of Education and Training, Hanoi, Vietnam
˄ Khon Kaen University, Khon Kaen, Thailand
˅ University of Michigan, Michigan, United States of America
Open Access

Abstract

Convergent morphological specializations for an arboreal lifestyle in most species of the Cyrtodactylus brevipalmatus group have been a confounding factor for establishing a stable taxonomy among its species. Recent references to C. interdigitalis from throughout Thailand and Laos were made without comparisons to the type material from Tham Yai Nam Nao, Nam Nao National Park, Phetchabun Province, Thailand, but instead, were based on general morphological similarity and distribution. The taxonomy of C. interdigitalis is stabilized here by comparing the paratypes to other specimens from Thailand and Laos and recovering their phylogenetic relationships based on newly acquired genetic data, including those from the type locality. The phylogeny recovered all specimens outside the type locality to be either C. ngati from Vietnam or new species closely related to C. ngati. Cyrtodactylus interdigitalis is shown here to be a range-restricted upland endemic on the Phetchabun massif of northern Thailand. The phylogeny also indicates that C. ngati extends hundreds of kilometers farther south into northern Thailand and central Laos. We hypothesize that the significant morphological divergence in body shape of the types of C. ngati, compared to that of the Lao and Thai populations, may be due to local adaptions for utilizing karst (C. ngati) rather than vegetation (Lao and Thai populations). Additionally, phylogenetic and multivariate analyses identified a potentially new species from Phu Hin Rong Kla National Park, Phitsanulok Province, in northern Thailand and another from the Khlong Naka Wildlife Sanctuary, Ranong Province, in southern Thailand. A series of newly examined specimens from Kaeng Krachan National Park, Phetchaburi Province, Thailand represents a possible ~82 km range extension to the southeast of C. rukhadeva. This research continues to underscore the high diversity of range-restricted upland endemics in Thailand and the importance of examining type material (if possible) in the context of a phylogeny so as to construct proper taxonomies that reveal, rather than obscure, diversity.

Keywords

Convergent evolution, Gekkota, integrative taxonomy, multivariate analysis, phylogenetics, Thailand, Vietnam

Introduction

The bent-toed gecko genus Cyrtodactylus is an ecologically and morphologically diverse gekkotan lineage (Grismer et al. 2020a) that ranges from the Himalayan uplands and peninsular India and Sri Lanka eastward to western Melanesia (Grismer et al. 2021a). This wide distribution across such a varied ecological landscape has seeded the evolution of 32 geographically localized monophyletic species groups (Grismer et al. 2021b), and of these, the brevipalmatus group of Indochina and Southeast Asia, is one of the most ecomorphologically specialized (Grismer et al. 2020a). It currently contains five described and at least three undescribed species (Grismer et al. 2021c) bearing a suite of unique morphological adaptations suited for an arboreal lifestyle, including a prehensile tail and cryptic coloration (Grismer et al. 2020a, 2021b). High degrees of ecomorphological specialization within Cyrtodactylus are often accompanied with high degrees of convergent evolution (Grismer 2020a, 2021b; Kaatz et al. 2021) and for a long time, this had been a confounding factor for establishing a stable taxonomy among species within the brevipalmatus group (Smith 1935; Welch et al. 1990; Ulber 1993; Manthey and Grossmann 1997; Stuart 1999; Nabhitabhata et al. 2004; Nabhitabhata and Chan-ard 2005; Pauwels and Chan-ard 2006; Grismer 2008; Ellis and Pauwels 2012). An integrative taxonomic analysis of Grismer et al. (2021c) was a recent attempt to disentangle this taxonomy but fell short due to a lack of specimens from the type locality of the putatively widespread C. interdigitalis Ulber, 1993 (e.g. Manthey and Grossmann 1997; Chan-ard et al. 1999, 2015; David et al. 2004, 2011; Cox et al. 1998). As such, Grismer et al. (2021c) referred to all Lao and Thai populations outside the type locality as C. cf. interdigitalis. Here, new morphological data from the paratypes of C. interdigitalis are compared to populations of C. cf. interdigitalis and to all other species of the brevipalmatus group in the context of a well-supported molecular phylogeny that includes genetic sequence data from the type locality of C. interdigitalis from Tham Yai Nam Nao, Nam Nao National Park, Phetchabun Province, Thailand (Fig. 1). With these new data, the taxonomy of C. interdigitalis and other populations is evaluated and the morphological and genetic discordance within C. ngati is discussed.

Figure 1. 

Distribution of nominal species and unnamed populations and specimens of the Cyrtodactylus brevipalmatus group. A single asterisk denotes specimens in the phylogenetic tree (Fig. 5) but not examined, two asterisks denote specimens examined but not in phylogenetic tree, and three asterisks denote specimens in phylogenetic tree and examined. Colored squares are literature localities and species identification based on location. White squares are literature localities from which specimens were not examined and remain unidentified. Locality data for all material examined appears in Table 1.

Methods

Species delimitation

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.

Sampling

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.

Table 1.

Material examined in this study. Institutional abbreviations follow Sabaj (2020) except that YC = Yodchaiy Chuankern from Department of Biology, Faculty of Science, Khon Kaen University, Mueang, Khon Kaen, 40002 Thailand.

Species Catalog no. Location GenBank no.
C. brevipalmatus LSUHC 1899 Thailand, no data not in tree
C. brevipalmatus THNHM 10670 Thailand, Nakhon Si Thammarat Province, Nopphitam District, Khao Nan National Park, Huay Lak Protected Unit not in tree
C. brevipalmatus THNHM 14112 Thailand, Nakhon Si Thammarat Province, Lan Saka District, Khao Luang National Park not in tree
C. brevipalmatus AUP-00573 Thailand, Nakhon Si Thammarat Province, Khao Ram Mt. OK626313
C. brevipalmatus LSUHC 11788 Penisular Malaysia, Kedah State, Pulau Langkawi, Gunung Raya not in tree
C. cf. brevipalmatus USMHC 2555 Penisular Malaysia, Kedah State, Pulau Langkawi, Gunung Raya OK626314
C. elok ZRC 2.6091 Penisular Malaysia, Pahang State, Fraser‘s Hill, the Gap JQ889180
C. elok LSUHC 12180 Penisular Malaysia, Pahang State, near Cameron Highlands not in tree
C. elok LSUHC 12181 Penisular Malaysia, Pahang State, near Cameron Highlands not in tree
C. elok ZMMU R-16144 Malaysian pet trade, no data not in tree
C. interdigitalis THNHM 20226 paratype Thailand, Phetchabun Province, Nam Nao National Park, Tham Yai Nam Nao not in tree
C. interdigitalis THNHM 20228 paratype Thailand, Phetchabun Province, Nam Nao National Park, Tham Yai Nam Nao not in tree
C. interdigitalis THNHM 20229 paratype Thailand, Phetchabun Province, Nam Nao National Park, Tham Yai Nam Nao not in tree
C. interdigitalis THNMH 20227 paratype Thailand, Phetchabun Province, Nam Nao National Park, Tham Yai Nam Nao not in tree
C. interdigitalis YC000952 Thailand, Phetchabun Province, Nam Nao National Park, Tham Yai Nam Nao ON055281
C. cf. interdigitalis ZMMU R-16492 Thailand, Phitsanulok Province, Phu Hin Rong Kla National Park MW792061
C. ngati FMNH 255454 Laos, Khammouane Province, Phou Hin Poun National Biodiversity Conservation Area JQ889181
C. ngati FMNH 270493 Laos, Khammouane Province, Phou Hin Poun National Biodiversity Conservation Area not in tree
C. ngati FMNH 270492 Laos, Khammouane Province, Phou Hin Poun National Biodiversity Conservation Area OK626315
C. ngati FMNH 265806 Thailand, Loei Province, Nam San Noi River, Phu Luang Wildlife Sanctuary JX51947
C. ngati HNUE-R00111 holotype Vietnam, Dien Bien Province, Dien Bien District, Pa Thom Commune, Pa Xa Lao Village, Karst forest near Pa Thom Cave ON411220
C. ngati IEBR 4829 paratype Vietnam, Dien Bien Province, Dien Bien District, Pa Thom Commune, Pa Xa Lao Village, Karst forest near Pa Thom Cave OK626318
C. ngati VNUF R.2020.12 paratype Vietnam, Dien Bien Province, Dien Bien District, Pa Thom Commune, Pa Xa Lao Village, Karst forest near Pa Thom Cave OK626319
C. ngati HNUE-R00112 paratype Vietnam, Dien Bien Province, Dien Bien District, Pa Thom Commune, Pa Xa Lao Village, Karst forest near Pa Thom Cave not in tree
C. ngati VNUF R.2014.50 Laos, Khammoue Province, Hin Nam No National Protected Area ON411221
C. cf. ngati 1 NCSM 79472 Laos, Xaignabouli Province, Ban Pha Liep, Houay Liep Stream OK626316
C. cf. ngati 2 ZMMU R-14917 Vientiane Province, Laos not in tree
C. cf. ngati 2 NCSM 80100 Laos, Vientiane Province, tributary of Nam Pha River, Houay Wan Stream OK626317
C. cf. rukhadeva THNHM 24622 Thailand, Phetchaburi Province, Kaeng Krachan National Park not in tree
C. cf. rukhadeva THNHM 24838 Thailand, Phetchaburi Province, Kaeng Krachan National Park not in tree
C. cf. rukhadeva THNHM 03251 Thailand, Phetchaburi Province, Kaeng Krachan National Park not in tree
C. cf. rukhadeva THNHM 03252 Thailand, Phetchaburi Province, Kaeng Krachan National Park not in tree
C. cf. rukhadeva THNHM 03253 Thailand, Phetchaburi Province, Kaeng Krachan National Park not in tree
C. cf. rukhadeva THNHM 03254 Thailand, Phetchaburi Province, Kaeng Krachan National Park not in tree
C. cf. rukhadeva THNHM 01807 Thailand, Phetchaburi Province, Kaeng Krachan National Park not in tree
C. rukhadeva ZMMU R-16851 holotype Thailand, Ratchaburi Province, Suan Phueng District, Khao Laem Mountain OK626320
C. rukhadeva ZMMU R-16852 paratype Thailand, Ratchaburi Province, Suan Phueng District, Hoop Phai Tong not in tree
sp. 9 AUP-01715 Thailand, Kanchanaburi Province, Thong Pha Phum District, Thong Pha Phum National Park MT468909
sp. 13 THNHM 00104 Thailand, Tak Province, Umphang District, Thung Yai Naresuan Wildlife Sanctuary not in tree
sp. 13 THNHM 27821 Thailand, Kanchanaburi Province, Thong Pha Phum District, Ban Saphan Lao not in tree
sp. 10 AUP-00680 Thailand, Tak Province, Tha Song Yang District, Mae Moei National Park, Chao Doi Waterfall MT468902
sp. 14 THNHM 01667 Thailand, Ranong Province, Khlong Nakha Wildlife Sanctuary not in tree
sp. 12 ZMMU R-16492 Phu Hin Rong Kla National Park, Phitsanulok, Province, Thailand ON411222

Morphological data

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.

Table 2.

part a. Sex, meristic, categorical, and raw morphometric data used in the analyses from specimens in the Cyrtodactylus brevipalmatus group. Abbreviations for morphometric characters are in the Materials and methods. / = data unavailable; m = male; f = female.

Species brevipalmatus brevipalmatus brevipalmatus brevipalmatus brevipalmatus elok elok elok elok interdigitalis interdigitalis interdigitalis interdigitalis cf. interdigitalis
Institutional catalog number LSUHC 1899 LSUHC 15076 LSUHC 11788 THNHM 10670 THNHM 14112 LSUHC 8238 LSUHC 12180 LSUHC 12181 ZMMU R-16144 paratype THNHM 20226 paratype THNHM 20228 paratype THNHM 20229 paratype THNHM 20227 ZMMU R-16492
sex m f f f f f m m f f f f f m
Meristic data
supralabials (SL) 11 12 10 14 12 11 8 13 9 14 12 11 12 8
infralabials (IL) 8 10 9 11 11 11 8 11 9 9 8 8 7 9
paravertebral tubercles (PVT) 39 37 38 29 23 0 0 0 0 27 27 24 23 30
longitudinal rows of tubercles (LRT) 15 16 17 16 14 6 7 4 4 19 20 19 19 18
ventral scales (VS) 38 38 38 36 39 45 45 47 36 42 40 42 43 34
ventral scales along middle of body (VSM) 176 170 182 154 160 190 225 234 192 187 170 187 178 160
expanded subdigital lamellae on 4th toe (TL4E) 7 8 9 8 8 10 9 9 9 12 10 10 11 9
unmodified subdigital lamellae on 4th toe (TL4U) 13 11 11 11 12 11 10 11 9 14 13 12 14 10
total subdigital lamellae 4th toe (TL4T) 20 19 20 19 20 21 19 20 18 26 23 22 23 19
expanded subdigital lamellae on 4th finger (FL4E) 8 8 8 7 8 8 9 9 9 9 8 9 9 10
unmodified subdigital lamellae on 4th finger (FL4U) 9 11 10 10 10 12 13 9 8 12 11 12 13 9
total subdigital lamellae 4th finger (FL4T) 17 19 18 17 18 20 22 18 17 21 21 21 22 19
enlarged femoral scales (R/L) 0 0 0 8R/ 8L 7R/ 7L 0 0 0 0 11R/8L 10R/9L 8R/ 8L 9R/10L 9R/8L
total enlarged femoral scales (FS) 16 10 11 16 14 0 0 0 0 14 19 19 19 17
total femoral pores (FP) 7 0 0 0 0 0 0 0 0 0 0 0 0 17
enlarged precloacal scales (PCS) 7 7 7 8 7 8 8 8 7 14 15 13 19 13
precloacal pores (PP) 7 0 0 0 0 0 8 8 0 0 0 0 0 13
postcloacal tubercles (PCT) 3 3 2 3 3 3 2 3 3 3 2 3 3 3
body bands (BB) 4 6 3 5 5 5 5 3 3 5 5 5 5 3
Categorical data
small tubercles on flank (FK-tub) present present present present present absent absent absent absent present present present present present
dorsolteral caudal tubercles (DCT) small small small / small large large large large small / small small large
ventrolateral caudal fringe (VLF1) small small small / small large large large large small / small small large
ventrolateral caudal fringe scales generally homogenous (VFL2) no no no / no no no no no yes yes yes yes yes
tail crossection (TLcross) circular circular circular / circular square square square square circular / circular circular square
slightly enlarged medial subcaudals (SC1) present present present / absent absent absent absent absent absent / absent absent /
single enlarged medial subcaudal (SC2) absent absent absent / absent absent absent absent absent absent / absent absent /
single enlarged medial subcaudals intermittent, medially furrowed, posteriorly emarginate (SC3) no no no / no no no no no yes / yes yes /
Morphomeric data
SVL 68.8 70.8 64.1 66.0 63.8 80.2 78.2 84.8 78.6 81.2 74.8 78.6 59.7 68.1
AG 35.7 33.4 30.1 30.0 26.5 39.7 37.8 41.5 36.2 34.5 33.7 32.7 24.6 34.6
HumL 9.7 9.3 8.0 9.6 9.5 10.2 9.1 10.1 1.7 9.8 10.2 11.2 7.4 10.3
ForL 9.9 9.8 8.9 8.2 8.7 11.5 11.7 11.8 10.2 10.6 10.5 11.1 8.4 8.5
FemL 12.0 12.6 11.5 11.7 9.8 12.9 14.2 14.6 13.1 14.7 13.2 12.7 10.2 12.6
TibL 11.6 12.2 10.5 9.7 8.2 13.5 14.0 13.8 12.3 13.1 11.9 12.9 10.2 11.4
HL 19.3 19.3 19.0 17.9 18.2 21.8 21.6 21.9 21.7 20.8 19.9 21.7 16.7 18.4
HW 13.2 13.8 12.3 12.3 12.0 15.6 16.1 15.9 15.1 14.0 13.4 14.2 11.4 13.1
HD 8.0 7.6 7.6 7.3 7.0 9.6 9.8 10.4 9.8 3.4 8.6 8.7 6.6 8.3
ED 5.2 4.5 4.3 5.3 4.4 4.8 5.0 5.7 5.0 5.3 5.5 5.9 4.4 4.4
EE 5.7 5.9 4.9 5.7 5.7 6.4 7.1 7.0 6.8 5.8 6.2 6.4 4.8 6.2
ES 7.4 7.6 7.0 7.0 7.2 8.6 8.7 9.5 8.6 8.3 7.8 9.1 6.8 7.7
EN 5.7 5.4 4.9 5.3 5.4 6.0 6.2 6.5 6.2 6.0 5.5 6.8 5.1 5.5
IO 5.4 4.7 4.7 4.2 5.2 5.7 5.4 5.4 3.9 4.8 4.7 5.5 4.3 2.9
EL 1.0 1.4 1.1 1.3 1.0 1.9 1.4 1.5 1.4 1.3 1.3 1.6 1.2 0.9
IN 1.7 2.1 2.3 2.1 2.2 2.7 2.6 2.5 3.1 2.1 2.2 2.5 1.8 2.3
Table 2.

part b. Sex, meristic, categorical, and raw morphometric data used in the analyses from specimens in the Cyrtodactylus brevipalmatus group. Abbreviations for morphometric characters are in the Materials and methods. / = data unavailable; m = male; f = female.

Species ngati ngati ngati ngati ngati ngati ngati ngati cf. ngati 1 cf. ngati 2 cf. ngati 2
Institutional catalog number holotype HNUE-R00111 paratype IEBR 4829 paratype VNUF R.2020.12 paratype HNUE-R00112 FMNH 255454 FMNH 270493 FMNH 270492 FMNH 265806 NCSM 79472 ZMMU R-14917 NCSM 80100
sex m f f f f m m m f f f
Meristic data
supralabials (SL) 10 10 10 10 13 13 13 10 14 9 12
infralabials (IL) 9 9 9 9 10 9 11 8 11 10 12
paravertebral tubercles (PVT) 39 40 38 40 28 27 26 27 28 32 29
longitudinal rows of tubercles (LRT) 18 18 17 22 19 18 17 19 18 24 19
ventral scales (VS) 38 36 35 32 37 36 36 33 33 36 35
ventral scales along middle of the body (VSM) 168 164 178 158 159 166 156 158 164 166 165
expanded subdigital lamellae on 4th toe (TL4E) 8 10 9 9 10 10 8 10 9 8 10
unmodified subdigital lamellae on 4th toe (TL4U) 11 10 11 10 11 11 11 11 12 10 10
total subdigital lamellae 4th toe (TL4T) 13 16 17 16 21 21 19 21 21 18 20
expanded subdigital lamellae on 4th finger (FL4E) 16 16 17 16 8 8 8 8 9 7 9
unmodified subdigital lamellae on 4th finger (FL4U) 9 9 9 9 10 10 10 10 8 9 10
total subdigital lamellae 4th finger (FL4T) 15 15 18 15 18 18 18 18 17 16 19
enlarged femoral scales (R/L) 10R/10L 9R/8L 10R/9L 8R/9L 9R/7L 8R/9L 9R/9L 8R/8L 9R/8L 7R/8L 7R/8L
total enlarged femoral scales (FS) 20 17 19 17 16 17 18 16 17 15 15
total femoral pores (FP) 14 0 0 0 0 14 15 13 0 0 0
enlarged precloacal scales (PCS) 13 13 13 13 15 13 13 13 12 13 13
precloacal pores (PP) 0 0 0 0 13 13 13 13 0 0 0
postcloacal tubercles (PCT) 3 2 1 2 0 0 0 0 2 3 4
body bands (BB) 6 6 6 6 3 4 3 3 3 3 3
Categorical data
small tubercles on flank (FK-tub) present present present present present present present present present present present
dorsolteral caudal tubercles (DCT) small small small small small small small small small small small
ventrolateral caudal fringe (VLF1) small small small small small small small small small small small
ventrolateral caudal fringe scales generally homogenous (VFL2) no no no no yes yes yes yes yes yes yes
tail crossection (TLcross) circular circular circular circular circular circular circular circular circular circular circular
slightly enlarged medial subcaudals (SC1) present present present present / present present present present present present
single enlarged medial subcaudal (SC2) absent absent absent absent / absent absent absent absent absent absent
single enlarged medial subcaudals intermittent, medially furrowed, posteriorly emarginate (SC3) no no no no / no no no no no no
Morphometric data
SVL 66.5 68.1 69.3 46.6 83.56 70.24 74.13 73.76 77.95 87.1 77.66
AG 28.8 29.8 30.2 19.7 41.3 35.4 37.0 31.3 38.2 41.9 36.8
HumL 7.9 8.1 8.5 5.6 8.6 8.7 8.6 6.9 8.7 11.5 9.2
ForL 9.2 10.0 10.1 6.5 10.2 9.3 10.4 10.0 10.3 10.4 10.7
FemL 11.5 11.5 11.5 7.6 13.7 12.7 13.0 13.1 13.1 15.2 14.2
TibL 10.8 11.1 11.8 7.8 12.5 11.8 11.2 11.1 12.8 12.6 12.7
HL 20.1 20.4 20.7 16.1 21.7 20.6 20.3 20.7 21.2 22.1 21.4
HW 12.6 12.0 11.8 8.8 13.8 12.5 13.0 12.3 12.7 14.8 13.5
HD 7.4 7.2 6.6 5.1 9.2 8.4 9.1 7.6 8.3 8.7 9.2
ED 3.8 4.1 3.4 2.6 4.9 4.9 4.9 4.8 6.5 4.6 6.0
EE 5.8 5.5 5.9 4.4 6.9 6.1 6.2 5.7 5.3 6.5 6.2
ES 7.5 7.6 6.9 5.0 9.0 8.3 8.3 8.2 8.7 8.8 8.4
EN 6.7 6.3 6.2 4.5 6.5 6.2 6.1 6.2 6.2 6.6 6.0
IO 5.6 5.4 5.6 4.2 6.6 5.6 5.4 5.1 4.9 3.5 5.7
EL 0.8 0.8 0.7 0.3 1.3 1.1 1.2 1.0 1.5 1.2 0.9
IN 2.8 2.6 2.6 2.0 2.8 2.5 2.5 2.3 2.7 2.2 2.5
Table 2.

part c. Sex, meristic, categorical, and raw morphometric data used in the analyses from specimens in the Cyrtodactylus brevipalmatus group. Abbreviations for morphometric characters are in the Materials and methods. / = data unavailable; m = male; f = female.

Species rukhadeva rukhadeva cf. rukhadeva cf. rukhadeva cf. rukhadeva cf. rukhadeva cf. rukhadeva cf. rukhadeva cf. rukhadeva sp. 13 sp. 13 sp. 14
Institutional catalog number Holotype ZMMU R-16851 Paratype ZMMU R-16852 THNHM 24622 THNHM 24838 THNHM 03251 THNHM 03252 THNHM 03253 THNHM 03254 THNHM 01807 THNHM 00104 THNHM 27821 THNHM 01667
sex m f m f m m f m m f f m
Meristic Data
supralabials (SL) 11 9 11 13 13 11 12 13 12 12 15 12
infralabials (IL) 10 11 10 10 10 10 10 11 10 10 11 10
paravertebral tubercles (PVT) 27 30 26 28 27 27 30 30 26 33 29 29
longitudinal rows of tubercles (LRT) 19 20 18 19 18 18 19 19 19 18 20 19
ventral scales (VS) 34 43 38 36 37 37 39 34 35 37 36 34
ventral scales along middle of the body (VSM) 154 152 162 158 157 159 168 160 161 159 165 159
expanded subdigital lamellae on 4th toe (TL4E) 9 9 8 9 9 10 9 10 10 9 7 8
unmodified subdigital lamellae on 4th toe (TL4U) 11 11 11 13 12 12 15 13 13 12 12 13
total subdigital lamellae 4th toe (TL4T) 20 18 19 22 21 22 14 23 23 21 19 21
expanded subdigital lamellae on 4th finger (FL4E) 9 8 7 8 8 8 8 8 8 8 8 8
unmodified subdigital lamellae on 4th finger (FL4U) 10 9 10 11 10 10 12 12 12 11 10 12
total subdigital lamellae 4th finger (FL4T) 19 17 17 17 18 18 20 20 20 19 18 20
enlarged femoral scales (R/L) 9R/8L 8R/8L 9R/L 9R/ 9L 9R/ 7L 7R/ 7L 6R/ 7L 5R/ 8L 7R/ 7L 9R/ 9L 7R/10L 7R/ 7L
total enlarged femoral scales (FS) 17 16 18 18 16 14 13 13 14 18 17 14
total femoral pores (FP) 17 0 14 0 12 13 0 11 13 0 0 7
enlarged precloacal scales (PCS) 17 13 15 15 14 13 15 15 14 14 16 16
precloacal pores (PP) 17 13 15 0 14 13 0 15 14 0 0 16
postcloacal tubercles (PCT) 3 2 3 2 3 2 2 3 2 3 3 3
body bands (BB) 3 3 3 3 4 4 / / 5 3 / 5
Categorical data
small tubercles on flank (FK-tub) present present present present present present present present present present present present
dorsolteral caudal tubercles (DCT) small small small small small small small small / small small large
ventrolateral caudal fringe (VLF1) small small small small small small small small / small small large
ventrolateral caudal fringe scales generally homogenous (VFL2) yes yes yes yes yes yes yes yes / yes yes no
tail crossection (TLcross) square square square square square square square square / circular circular /
slightly enlarged medial subcaudals (SC1) absent absent absent absent absent absent absent absent / present present absent
single enlarged medial subcaudal (SC2) present present present present present present present present / absent absent absent
single enlarged medial subcaudals intermittent, medially furrowed, posteriorly emarginate (SC3) no no no no no no no no no no no no
Morphomeric data
SVL 74.9 71.7 68.3 71.8 73.6 75.3 74.7 73.2 61.5 63.7 72.9 70.2
AG 34.6 32.6 27.3 29.9 30.9 31.3 32.2 30.3 26.2 25.8 30.6 31.5
HumL 10.7 10.4 9.8 8.3 12.2 11.3 11.8 11.0 10.1 7.6 10.1 10.2
ForL 8.6 7.9 8.7 8.5 9.0 10.6 9.6 9.2 7.9 8.1 9.6 8.6
FemL 12.6 11.8 10.8 10.9 11.5 10.2 11.9 12.1 9.5 10.7 12.8 12.1
TibL 10.1 9.3 9.7 10.7 10.9 11.7 11.3 11.1 9.1 10.1 10.2 11.8
HL 20.2 19.2 19.7 19.9 20.8 21.3 20.8 21.5 17.9 17.6 19.9 18.3
HW 14.6 13.4 13.1 13.9 14.9 15.0 13.1 14.1 11.8 11.9 13.8 12.1
HD 9.2 8.5 7.3 8.9 8.2 8.2 8.1 8.9 7.5 7.7 8.4 7.8
ED 4.6 4.3 4.9 5.1 5.8 5.4 5.0 5.5 4.7 4.1 5.3 5.2
EE 6.2 6.2 5.1 6.2 5.6 5.7 5.4 6.2 4.3 4.9 6.3 4.9
ES 8.3 7.7 7.4 8.1 8.4 8.8 8.1 8.6 7.3 7.2 8.0 7.5
EN 6.3 5.7 5.4 6.0 6.2 6.4 5.8 6.2 5.3 5.6 5.9 5.5
IO 3.3 3.1 4.5 4.7 5.6 5.7 5.7 5.6 4.2 4.8 6.1 4.0
EL 1.2 1.0 1.6 1.5 1.2 1.3 1.2 1.2 0.9 1.4 1.4 1.3
IN 2.2 2.1 2.0 2.2 2.4 2.5 2.4 2.3 2.0 2.1 2.3 2.2

Phylogenetic analyses

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.

Statistical analyses

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: Xadj=log(X)-β[log(SVL)-log(SVLmean)], where Xadj=adjusted value; X=measured value; β=unstandardized regression coefficient for each population; and SVLmean=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 padj<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 R2, R2adj, 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 R2 and padj-values were compared to the results of the Mantel and dbRDA tests.

Results

Phylogenetic data

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.

Multivariate data

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.

Table 3.

Summary statistics of the PERMANOVA test. Shaded cells denote significant p values (< 0.05)

species pairs F Model R 2 p value p-adjusted
rukhadeva vs ngati 11.2567 0.5058 0.0014 0.0142
rukhadeva vs interdigitalis 4.6290 0.2962 0.0016 0.0158
rukhadeva vs ngati types 9.1022 0.4765 0.0046 0.0458
rukhadeva vs brevipalmatus 7.9317 0.3979 0.0006 0.0056
ngati vs interdigitalis 3.8323 0.3898 0.0286 0.2857
ngati vs ngati types 15.6712 0.7581 0.0286 0.2857
ngati vs brevipalmatus 18.0303 0.7203 0.0080 0.0804
interdigitalis vs ngati types 7.7098 0.6066 0.0286 0.2857
interdigitalis vs brevipalmatus 5.9252 0.4584 0.0082 0.0820
ngati types vs brevipalmatus 5.9327 0.4972 0.0179 0.1786

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.

Table 4.

PCA summary statistics for Cyrtodactylus ngati and C. cf. ngati 1 and 2. Shaded cells denote heavy loadings. Abbreviations are in the Materials and methods.

PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 3.02514 1.43855 1.31700 1.05543 0.94514 0.75287 0.64766 0.22607 0.00000 0.00000
Proportion of variance 0.57197 0.12934 0.10841 0.06962 0.05583 0.03543 0.02622 0.00319 0.00000 0.00000
Cumulative proportion 0.57197 0.70131 0.80971 0.87933 0.93516 0.97059 0.99681 1.00000 1.00000 1.00000
Eigenvalue 9.15147 2.06943 1.73449 1.11394 0.89329 0.56681 0.41947 0.05111 0.00000 0.00000
SVL –0.24017 0.02365 –0.05714 –0.04626 0.38558 –0.03327 0.88620 –0.03510 0.00000 0.00000
AG –0.31222 –0.03444 –0.16700 –0.13169 –0.04196 0.24410 –0.07073 0.08044 –0.09704 –0.44157
HumL –0.18395 0.30455 –0.04150 –0.10087 –0.70280 0.16579 0.23479 –0.28513 0.10144 –0.10135
ForL –0.20015 –0.17542 0.37514 0.50929 0.06400 0.22482 –0.00663 0.29195 0.33516 –0.05149
FemL –0.31397 0.15040 0.00333 0.10602 0.09091 –0.19909 –0.12166 0.22019 –0.43592 –0.21044
TibL –0.26011 –0.21195 0.05811 –0.33668 –0.28102 –0.36533 0.04617 0.41117 –0.28757 0.16460
HL –0.29184 0.00874 –0.02861 –0.13911 0.05482 –0.57390 –0.13660 –0.06990 0.65217 –0.23562
HW –0.29884 0.18474 –0.19528 0.18133 –0.09028 –0.00066 –0.04357 0.15712 0.11422 0.76179
HD –0.30883 –0.11021 –0.13985 0.18125 –0.07375 0.20441 –0.04838 –0.19786 –0.14403 –0.03920
ED –0.28495 –0.23283 0.25136 –0.15942 0.04516 –0.02978 –0.08982 –0.14963 –0.07552 0.08892
EE –0.13329 0.05995 –0.61810 0.38078 –0.05532 –0.05066 –0.02928 0.15765 0.06903 –0.15484
ES –0.30780 –0.06844 –0.06403 –0.13606 0.24823 0.07530 –0.22076 –0.57475 –0.01069 0.20829
EN 0.20879 0.23130 –0.41436 –0.33049 0.26567 0.06299 –0.10032 0.16319 0.09653 0.05554
IO 0.07527 –0.57805 –0.32974 0.20819 0.01806 –0.17406 0.00256 –0.21397 –0.17053 0.00329
EL –0.26849 –0.09245 –0.05184 –0.34304 0.19539 0.49703 –0.14735 0.26629 0.13461 –0.01479
IN 0.14076 –0.54969 –0.17904 –0.20525 –0.26817 0.16311 0.15877 0.13640 0.26272 0.06895

Mismatch analysis and Tajima’s D

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; pnorm= 6.75e-05 and pbeta = 0.00 for the interdigitalis clade and D = –6.2531; pnorm = 4.071e-10 and pbeta = 0.00 for C. ngati.

Isolation by distance analyses

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 R2 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 (R2=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 R2adj=0.190 for the interdigitalis clade and p=0.233 and R2adj=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.

Discussion

The lack of correlations between the genetic and geographic (dis)similarity matrices recovered in the Mantel and dbRDA tests for both data sets are consistent with a phylogeographic structure that is not a function of IBD (e.g. Jombart 2021; Chan et al. 2022). This is consistent with insignificant negative Tajima’s D statistics, suggesting the lineages are not undergoing range expansion. However, Teske et al. (2018) noted that mitochondrial data is generally not as reliable as microsatellite data at recovering IBD, especially when individuals from different populations are pooled. Although our analyses did not pool individuals, we are fully aware of this and other issues surrounding the ongoing debates concerning the use of Mantel tests and their application to genetic and spatial distance matrices to recover IBD (see Diniz-Filho et al. 2013, Legendre et al. 2015 and Jombart 2021 and references therein). However, we present our results simply as baseline hypotheses to be tested following the acquisition of genomic data and additional specimens from the ~300–525 km geographic hiatus among some of the individuals of C. ngati.

The 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) is especially perplexing given that the uncorrected pairwise sequence divergence within C. ngati is no higher than 1.3% (Fig. 5). This, coupled with their difference in habitat preference (karst (type series) versus vegetation (the other specimens)), would suggest all members of the interdigitalis clade—based on their phylogenetic relationships (Fig. 5)—retain an ancestral morphology for living in vegetation except for the type series of C. ngati whose divergent morphology may be due to a derived habitat preference for seeking refuge in the crevices of karstic outcroppings (Le et al. 2021). Given that C. ngati from the type locality has a significantly smaller body, a smaller and flatter head with smaller eyes, a shorter snout, and smaller ear openings, as well as shorter hind limbs (Fig. 7C; Table 5), leads us to hypothesize that reduction in these particular morphometric traits (SVL, AG, HL, HW, and HD) may be an adaption for taking refuge in narrow crevices of karst formations (Le et al. 2021; Fig. 8) as seen in other species of Cyrtodactylus from different species groups (Grismer et al. 2016; Nielson and Oliver 2017; Quah et al. 2019; Wood et al. 2020; Kaatz et al. 2021). However, longer limbs which is also a feature of karst associated species, was not observed. Because this is the only known transition of an arboreal species of Cyrtodactylus to rock crevice-dwelling species (Grismer et al. 2021b), no intra-generic comparisons of this specific transition type can be made. However, Grismer (2021) demonstrated that the crevice-dwelling xantusiid lizard Xantusia henshawi Stejneger, differed statistically from the more generalized X. gracilis Grismer and Galvan in also having a relatively smaller head (HW) and eyes (ED) and shorter hind limbs (FemL and TibL). The decoupling of morphological and mitochondrial DNA divergences in karst-adapted species of Cyrtodactylus at the shallow nodes in trees, is not an uncommon theme (e.g. Grismer et al. 2016; Nielson and Oliver 2017; Quah et al. 2019; Wood et al. 2020; Kaatz et al. 2021) and is commensurate with the general ecological plasticity of this genus (Grismer et al. 2020a).

Table 5.

t-test summary statistics for Cyrtodctylus ngati and C. cf. ngati 1 and 2. Shaded cells denote significant p values. Abbreviations are in the Materials and methods.

character Welch’s t Student t p value
SVL 2.9327 0.0189
AG 10.47 3.91E-05
HumL 1.1501 0.2833
ForL 1.7757 0.1137
FemL 13.081 1.235E-05
TibL 2.7255 0.026
HL 5.1505 0.0009
HW 7.6325 0.0001
HD 7.3555 8.20E-05
ED 3.8713 0.0047
EE 1.3404 0.2169
ES 7.1308 0.0107
EN 1.4927 0.1727
IO -64882 0.5346
EL 4.5441 0.0019
IN -1.5426 0.1615

An alternative, but less preferred hypothesis, is that mitochondrial DNA introgression (through hybridization) or incomplete lineage sorting (very recent speciation precluding the accumulation of mutations) is obscuring species boundaries in this group (e.g., McGuire et al. 2007), and the divergent morphology instead reflects non-conspecificity of the types of C. ngati with the Lao and Thai populations despite sharing mitochondrial haplotypes (Fig. 5A). No cases of hybridization (necessary for mitochondrial introgression) are yet known in Cyrtodactylus, although the presence of two divergent mitochondrial lineages within C. ziegleri at its type locality in Vietnam provides for the possibility that mitochondrial DNA may imperfectly track species boundaries in other Cyrtodactylus (see Neang et al. 2020).

These analyses demonstrate the necessity of examining type material (when possible) in the context of a phylogeny in order to correctly identify similarly appearing specimens (see Fig. 10) outside the type locality in order to construct taxonomies that reflect evolutionary history as opposed to those that obscure it. This is especially true for highly specialized species where convergent morphologies do not align with phylogenetic history (e.g. Grismer et al. 2020a; Grismer et al. 2020b; Kaatz et al. 2021). Analyses of morphological data from the paratypes of C. interdigitalis and genetic data from a recently acquired specimen from the type locality at Tham Yai Nam Nao, Nam Nao National Park, Phetchabun Province, Thailand clearly demonstrated that C. interdigitalis is not a widespread species ranging across Thailand and Vietnam as previously thought, but is a range-restricted endemic in the upland regions of the Phethchabun massif. This pattern of upland endemism is common across many Thai species of amphibians and reptiles (e.g. Matsui 2006; Connette et al. 2017; Sumontha et al. 2012, 2017; Wilkinson et al. 2012; Grismer et al. 2017; 2020b; Matsui et al. 2018; Pawangkhanant et al. 2018; Suwannapoom et al. 2018, 2021; Lee et al. 2019; Chomdej et al. 2021 and references therein) and will become even more common as fieldwork continues in these sky-island regions. As climate change continues to have a significant negative impact on range-restricted upland species of amphibians and reptiles world-wide (see Sinervo et al. 2010; Grant et al. 2020 and references therein), it is paramount that field work in these hyper-biodiverse regions of Indochina and Southeast Asia remains an active endeavor and that the publication of new results not be delayed.

Figure 10. 

A Cyrtodactylus interdigitalis from Phetchabun Province, Thailand. Photo by Montri Sumonta. B Cyrtodactylus cf. interdigitalis (ZMMU R-16492) from Phu Hin Rong Kla National Park, Phitsanulok Province, Thailand. Photo by Nikolay A. Poyarkov. C Cyrtodactylus cf. ngati 2 (NCSM 80100) from Houay Wan Stream, tributary of Nam Pha River, Vientiane Province, Laos. Photo by Bryan L. Stuart. D Cyrtodactylus ngati (VNUF R.2020.12) from Pa Thom Cave, Pa Xa Lao Village, Pa Thom Commune, Dien Bien District, Dien Bien Province, Vietnam. Photo by Dzung T. Le. E Cyrtodactylus cf. ngati 1 (NCSM 79472) Houay Liep Stream, Paklay: Ban Pha Liep, Xaignabouli Province, Laos. Photo by Bryan L. Stuart. F Adult male C. cf. rukhadeva from Kaeng Krachan National Park, Phetchaburi Province, Thailand. Photo from Creative Commons Attribution Share Alike.

Acknowledgements

This work (Grant No. RGNS 64-038) was financially supported by Office of the Permanent Secretary, Ministry of Higher Education, Science, Research and Innovation. AR and AA were supported by Kasetsart University Research and Development Institute (KURDI). YC was financial supported by the Research Administration Division (RAD), Khon Kaen University. Wachara Sanguansombat and Sunchai Makchai (Thailand Natural History Museum) made specimens in their care available for study. Sengvilay Seateun, Khampong Thanonkeo, Tanya Chan-ard and David Emmett assisted with collecting NCSM and FMNH specimens in Laos.

References

  • Acinas SG, Klepac-Ceraj V, Hunt DE, Pharino C, Ceraj I, Distel DL, Polz MF (2004) Fine-scale phylogenetic architecture of a complex bacterial community. Nature 430: 551–554. https://doi.org/10.1038/nature02649
  • Agarwal A, El-Ghazawi T, El-Askary H, Le-Moigne J (2007) Efficient hierarchical-PCA dimension reduction for hyperspectral imagery. 2007 IEEE International Symposium on Signal Processing and Information Technology. https://doi.org/10.1109/isspit.2007.4458191
  • Barraclough TG, Birky Jr CW, Burt A (2003) Diversification in sexual and asexual organisms. Evolution 57: 2166–2172. https://doi.org/10.1554/02-339
  • Chan KO, Grismer LL (2022) GroupStruct: an R package for allometric size correction. Zootaxa 5124: 471–482.
  • Chan KO, Hutter CR, Wood Jr PL, Su Y-C, Brown RM (2022) Gene flow increases phylogenetic structure and inflates cryptic species estimations: a case study on widespread Philippine puddle frogs (Occidozyga laevis). Systematic Biology 71: 40–57. https://doi.org/10.1093/sysbio/syab034
  • Chan-ard T, Grossmann W, Gumprecht A, Schulz K-D (1999) Amphibians and Reptiles of Peninsular Malaysia And Thailand; An Illustrated Checklist. Bushmaster Publications, Wuerselen, 240 pp.
  • Chan-ard T, Parr JWK, Nabhitabhata J (2015) A Field Guide to the Reptiles of Thailand. Oxford University Press, New York, 314 pp.
  • Chomdej S, Pradit W, Suwannapoom C, Pawangkhanant P, Nganvongpanit K, Poyarkov NA, Che J, Gao Y, Gong S (2021) Phylogenetic analyses of distantly related clades of benttoed geckos (genus Cyrtodactylus) reveal an unprecedented amount of cryptic diversity in northern and western Thailand. Scientific Reports 11: 2328. https://doi.org/10.1038/s41598-020-70640-8
  • Connette GM, Oswald P, Thura MK, Connette KJL, Grindley ME, Songer M, Zug GR, Mulcahy DG (2017) Rapid forest clearing in a Myanmar proposed national park threatens two newly discovered species of geckos (Gekkonidae: Cyrtodactylus). PLoS One 12: e0174432. https://doi.org/10.1371/journal.pone.0174432
  • Cox MJ, van Dijk PP, Nabhitabhata J, Thirakhupt K (1998) A Photographic Guide to Snakes and Other Reptiles of Peninsular Malaysia, Singapore and Thailand. New Holland Publishers (UK) Ltd., London, 144 pp.
  • David P, Teynié A, Ohler A (2004) A new species of Cyrtodactylus Gray, 1827 (Reptilia: Squamata: Gekkonidae) from southern Laos. Raffles Bulletin of Zoology 52: 621–627.
  • David P, Nguyen TQ, Schneider N, Ziegler T (2011) A new species of the genus Cyrtodactylus Gray, 1827 from central Laos (Squamata: Gekkonidae). Zootaxa 2833: 29–40. https://doi.org/10.11646/zootaxa.2833.1.3
  • Ellis M, Pauwels OSG (2012) The bent-toed geckos (Cyrtodactylus) of the caves and karst of Thailand. Cave and Karst Science 39: 16–22.
  • Frost DR, Hillis DM (1990) Species in concept and practice: herpetological application. Herpetologica 46: 87–104.
  • Fujisawa T, Barraclough TG (2013) Delimiting species using single-locus data and the Generalized Mixed Yule Coalescent approach: a revised method and evaluation on simulated data sets. Systematic Biology 62: 707–724. https://doi.org/10.1093/sysbio/syt033
  • Grismer LL (2008) On the distribution and identification of Cyrtodactylus brevipalmatus Smith, 1923 and Cyrtodactylus elok Dring, 1979. The Rafflles Bulletin of Zoology 56: 177–179.
  • Grismer LL (2021) Comparative ecomorphology of the sandstone night lizard (Xantusia gracilis) and the granite night lizard (Xantusia henshawi). Vertebrate Zoology 71: 425–437. https://doi.org/10.3897/vz.71.e69214
  • Grismer LL, Suwannapoom C, Pawangkhanant P, Nazarov RA, Yushchenko PV, Naiduangchan M, Le MD, Luu VQ, Poyarkov NA (2021c) A new cryptic arboreal species of the Cyrtodactylus brevipalmatus group (Squamata: Gekkonidae) from the uplands of western Thailand. Vertebrate Zoology 71: 723–746. https://doi.org/10.3897/vz.71.e76069
  • Grismer LL, Wood Jr PL, Anuar S, Grismer MS, Quah ESH, Murdoch ML, Muin MA, Davis HR, Aguilar C, Klabacka R, Cobos AJ, Aowphol A, Sites Jr JW (2016) Two new Bent-toed Geckos of the Cyrtodactylus pulchellus complex from Peninsular Malaysia and multiple instances of convergent adaptation to limestone forest ecosystems. Zootaxa 4105: 401–429. https://doi.org/10.11646/zootaxa.4105.5.1
  • Grismer LL, Wood Jr PL, Aowphol A, Cota M, Grismer MS, Murdoch ML, Aguilar C, Grismer JL (2017) Out of Borneo, again and again: biogeography of the Stream Toad genus Ansonia Stoliczka (Anura: Bufonidae) and the discovery of the first limestone cave-dwelling species. Biological Journal of the Linnean Society, 120: 371–395. https://doi.org/10.1111/bij.12886
  • Grismer LL, Wood Jr PL, Le MD, Quah ESH, Grismer JG (2020a) Evolution of habitat preference in 243 species of Bent-toed Geckos (Genus Cyrtodactylus Gray, 1827) with a discussion of karst habitat conservation. Ecology and Evolution 10: 13717–13730. https://doi.org/10.1002/ece3.6961
  • Grismer LL, Wood Jr PL, Poyarkov NA, Le MD, Karunarathna S, Chomdej S, Suwannapoom C, Qi S, Liu S, Che J, Quah ESH, Kraus F, Oliver PM, Riyanto A, Pauwels OSG, Grismer JL (2021b) Karstic landscapes are foci of species diversity in the world’s third-largest vertebrate Genus Cyrtodactylus Gray, (Reptilia: Squamata; Gekkonidae). Diversity 13: 183. https://doi.org/10.3390/d13050183
  • Grismer LL, Wood Jr PL, Poyarkov NA, Minh DL, Kraus F, Agarwal I, Oliver PM, Nguyen S, Nguyen TQ, Karunarathna S, Welton LJ, Stuart BL, Luu VQ, Bauer AM, O’Connell KA, Quah ESH, Chan KO, Ziegler T, Ngo H, Nazarov RA, Aowphol A, Chomdej S, Suwannapoom C, Siler CD, Anuar S, Tri NV, Grismer JL (2021a) Phylogenetic partitioning of the third-largest vertebrate genus in the world, Cyrtodactylus Gray, 1827 (Reptilia; Squamata; Gekkonidae) and its relevance to taxonomy and conservation. Vertebrate Zoology 71: 101–154. https://doi.org/10.3897/vertebrate-zoology.71.e59307
  • Grismer LL, Yushchenko PV, Pawangkhanant P, Naiduangchan M, Nazarov RA, Orlova VF, Suwannapoom C, Poyarkov NA (2020b) A new species of Hemiphyllodactylus Bleeker (Squamata; Gekkonidae) from Peninsular Thailand that converges in morphology and color pattern on Pseudogekko smaragdinus (Taylor) from the Philippines. Zootaxa 4816 2: 171–190. https://doi.org/10.11646/zootaxa.4816.2.2
  • Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS (2018) UFBoot2: Improving the ultrafast bootstrap approximation. Molecular Biology and Evolution 35: 518–522. https://doi.org/10.1093/molbev/msx281
  • Huelsenbeck JP, Ronquist F, Nielsen R, Bollback J.P (2001) Bayesian Inference of phylogeny and its impact on evolutionary biology. Science 294: 2310–2314. https://doi.org/10.1126/science.1065889
  • Husson F, Josse J, Le S, Mazet J (2017) FactoMine R: Exploratory Data Analysis and Data Mining. R package version 1.36.
  • Kaatz A, Grismer JL, Grismer LL (2021) Convergent evolution of karst habitat preference and its ecomorphological correlation in three species of Bent-toed Geckos (Cyrtodactylus) from Peninsular Malaysia. Vertebrate Zoology 71: 367–386. https://doi.org/10.3897/vz.71.e66871
  • Kalyaanamoorthy S, Minh BQ, Wong TK, von Haeseler A, Jermiin LS (2017) ModelFinder: fast model selection for accurate phylogenetic estimates. Nature Methods 14: 587. https://doi.org/10.1038/nmeth.4285
  • Kassambara A, Mundt F (2017) factoextra: Extract and Visualize the Result of Multivariate Data Analyses. R package version 1.0.5.999.
  • Katoh M, Misawa K, Kuma K, Miyata T (2002) MAFTT: A novel method for rapid sequence alignment based on fast Fourier transform. Nucleic Acids Research 30: 3059–3066. https://doi.org/10.1093/nar/gkf436
  • Kumar S, Stecher G, Tamura K (2016) MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets. Molecular Biology and Evolution 33: 1870–1874. https://doi.org/10.1093/molbev/msw054
  • Le DT, Sitthivong S, Tran TT, Grismer LL, Nguyen TQ, Le MD, Ziegler T, Luu VQ (2021) First record of the Cyrtodactylus brevipalmatus group (Squamata: Gekkonidae) from Vietnam with description of a new species. Zootaxa 4969: 492–510. https://doi.org/10.11646/zootaxa.4969.3.3
  • Lee JL, Miller AH, Zug GR, Mulcahy DG (2019) The discovery of Rock Geckos Cnemaspis Strauch, 1887 (Squamata: Gekkonidae) in the Tanintharyi Region, Myanmar with the description of two new species. Zootaxa 466: 40–64. https://doi.org/10.11646/zootaxa.4661.1.2
  • Lin X-L, Stur E, Ekrem T (2018) Exploring species boundaries with multiple genetic loci using empirical data from non-biting midges. Zoologica Scripta 47: 325–341. https://doi.org/10.1111/zsc.12280
  • Lleonart J, Salat J, Torres GJ (2000) Removing allometric effects of body size in morphological analysis. Journal of Theoretical Biology 205: 85–93. https://doi.org/10.1006/jtbi.2000.2043
  • Manthey U, Grossmann W (1997) Amphibien und Reptilien Südostasiens. Natur und Tier Verlag, Münster, Germany, 512 pp.
  • Matsui M, Khonsue W, Panha S (2018) Two New Species of Ansonia from Thailand (Anura: Bufonidae). Zoological Science 35: 39–48. https://doi.org/10.2108/zs170120
  • McCoy MW, Bolker BM, Osenberg CW, Miner BG, Vonesh JR (2006) Size correction: Comparing morphological traits among populations and environments. Oecologia 148: 547–554. https://doi.org/10.1007/s00442-006-0403-6
  • McGuire JA, Linkem CW, Koo MS, Hutchison DW, Lappin AK, Orange DI, Lemos-Espinal J, Riddle BR, Jaeger JR (2007) Mitochondrial introgression and incomplete lineage sorting through space and time: phylogenetics of crotaphytid lizards. Evolution 61: 2879–2897. https://doi.org/10.1111/j.1558-5646.2007.00239.x
  • Miller MA, Pfeiffer W, Schwartz T (2010) Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In: Gateway Computing Environments Workshop (GCE), New Orleans (USA), November 2010, IEEE, 1–8. https://doi.org/10.1109/GCE.2010.5676129
  • Nabhitabhata J, Chan-ard T, Chuaynkern Y (2004) . Checklist of Amphibians and Reptiles in Thailand. Office of Environmental Policy and Planning, Bangkok, 152 pp.
  • Nabhitabhata J, Chan-ard T (2005) Thailand Red Data: Mammals, Reptiles and Amphibians. Office of Natural Resources and Environmental Policy and Planning, Bangkok, Thailand, 234 pp.
  • Neang T, Henson A, Stuart BL (2020) A new species of Cyrtodactylus (Squamata, Gekkonidae) from Cambodia’s Prey Lang Wildlife Sanctuary. ZooKeys 926: 133–158. https://doi.org/10.3897/zookeys.926.48671
  • Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ (2015) IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Molecular Biology and Evolution 32: 268–274. https://doi.org/10.1093/molbev/msu300
  • Nielsen SV, Oliver PM (2017) Morphological and genetic evidence for a new karst specialist lizard from New Guinea (Cyrtodactylus: Gekkonidae). Royal Society Open Science 4: 170781. https://doi.org/10.1098/rsos.170781
  • Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P, Stevens MHH, Szoecs E, Wagner H (2020) Package ‘vegan’. Version 2.5-7. https://cran.r-project.org/web/packages/vegan
  • Pagès J (2015) Multiple Factor Analysis by Example Using R. CRC Press, New York, 272 pp.
  • Pauwels OSG, Chan-ard T (2006) Reptiles of Kaeng Krachan National Park, western Thailand. Natural History Bulletin of the Siam Society 54: 89–108.
  • Pawangkhanant P, Poyarkov NA, Duong TV, Naiduangchan M, Suwannapoom C (2018) A New Species of Leptobrachium (Anura, Megophryidae) from western Thailand. PeerJ 6: e5584. https://dx.doi.org/10.7717/peerj.5584
  • Pons J, Barraclough TG, Gomez-Zurita J, Cardoso A, Duran DP, Hazell S, Kamoun S, Sumlin WD, Vogler AP (2006) Sequence-Based Species Delimitation for the DNA Taxonomy of Undescribed Insects. Systematic Biology 55: 595–609. https://doi.org/10.1080/10635150600852011
  • Quah ESH, Grismer LL, Wood Jr PL, Anuar SAM (2019) The discovery and description of a new species of Bent-toed Gecko of the Cyrtodactylus pulchellus complex (Squamata: Gekkonidae) from the Langkawi Archipelago, Kedah, Peninsular Malaysia. Zootaxa 4668: 51–75. https://doi.org/10.11646/zootaxa.4668.1.3
  • R Core Team. (2018) R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna. Available from: http://www.R-project.org (accessed 29 December 2018).
  • Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior summarization in Baysian phylogenetics using Tracer 1.7. Systematic Biology 67: 901–904. https://doi.org/10.1093/sysbio/syy032
  • Reist JD (1986) An empirical evaluation of coefficients used in residual and allometric adjustment of size covariation. Canadian Journal of Zoology 64: 1363–1368. https://doi.org/10.1139/z86-203
  • Rogers AR, Harpending H (1992) Population growth makes waves in the distribution of pairwise genetic differences. Molecular Biology and Evolution 9: 552–569.
  • Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012) MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61: 539–542. http://doi.org/10.1093/sysbio/sys029
  • Sinervo B, Méndez-de-la-Cruz F, Miles DB, Heulin E, Bastiaans B, Villagrán-Santa Cruz M, Lara-Resendiz R, Martínez-Méndez N, Calderón-Espinosa ML, Meza-Lázaro RN, Gadsden H, Avila LJ, Morando M, De la Riva IJ, Sepulveda PV, Rocha CFD, Ibargüengoyta N, Puntriano CA, Massot M, Lepetz V, Oksanen TA, Chapple DG, Bauer AM, Branch WR, Clobert J, Sites Jr JW (2010) Erosion of lizard diversity by climate change and altered thermal niches. Science 328: 894–899. https://doi.org/10.1126/science.1184695
  • Skalaski JR, Richins SM, Townsend RL (2018) A statistical test and sample size recommendations for comparing community composition following PCA. PLoS One 13: e0206003
  • Smith MA (1935) The fauna of British India, including Ceylon and Burma. Reptilia and Amphibia. Vol. II. Sauria. Taylor and Francis, London, 440 pp.
  • Stuart BL (1999) Amphibians and Reptiles. In: Duckworth JW, Salter RE, Khounboline K (Eds.) Wildlife in Lao PRD: 1999 Status Report. Vientiane, IUCN – The World Conservation Union – Wildlife Conservation Society – Centre for Protected Areas and Watershed Management, 43–67.
  • Sumontha M, Kunya K, Dangsri S, Pauwels OSG (2017) Oligodon saiyok, a new limestone dwelling kukri snake (Serpentes: Colubridae) from Kanchanaburi Province, western Thailand. Zootaxa 4294: 316–328. https://dx.doi.org/10.11646/zootaxa.4294.3.2
  • Sumontha M, Pauwels OSG, Kunya K, Limlikhitaksorn C, Ruksue S, Taokratok A, Ansermet M, Chanhome L (2012) A new species of Parachute Gecko (Squamata: Gekkonidae: genus Ptychozoon) from Kaeng Krachan National Park, western Thailand. Zootaxa 3513: 68–78. https://doi.org/10.11646/zootaxa.3513.1.4
  • Suwannapoom C, Grismer LL, Pawangkhanant P, Naiduangchan M, Yushchenko PV, Arkhipov DV, Wilkinson JA, Poyarkov NA (2021) Hidden tribe: A new species of Stream Tod of the genus Ansonia Stoliczka, 1870 (Anura: Bufonidae) from the poorly explored mountainous borderlands of western Thailand. Vertebrate Zoology 71: 763–779. https://doi.org/10.3897/vz.71.e73529
  • Suwannapoom C, Sumontha M, Tunprasert J, Ruangsuwan T, Pawangkhanant P, Korost DV, Poyarkov NA (2018) A striking new genus and species of cave-dwelling frog (Amphibia: Anura: Microhylidae: Asterophryinae) from Thailand. PeerJ 6: e4422. https://dx.doi.org/10.7717/peerj.4422
  • Talavera G, Dincă V, Vila R (2013) Factors affecting species delimitations with the GMYC model: insights from a butterfly survey. Methods in Ecology and Evolution 4: 1101–1110. https://doi.org/10.1111/2041-210X.12107
  • Teske PR, Golla TR, Sandoval-Castillo J, Emami-Khoyi A, van der Lingen CD, von der Heyden S, Chiazzari B, van Vuuren BJ, Beheregaray LB (2018) Mitochondrial DNA is unsuitable to test for isolation by distance. Scientific Reports 8:8448. https://doi.org/10.1038/s41598-018-25138-9
  • Thorpe RS (1975) Quantitative handling of characters useful in snake systematics with particular reference to intraspecific variation in the Ringed Snake Natrix natrix (L.). Biological Journal of the Linnean Society 7: 27–43. https://doi.org/10.1111/j.1095-8312.1975.tb00732.x
  • Thorpe RS (1983) A review of the numerical methods for recognising and analysing racial differentiation. In: Felsenstein J, ed. Numerical Taxonomy. NATO ASI Series (Series G: Ecological Sciences), vol. 1. Springer-Verlag, Berlin, 404–423. https://doi.org/10.1007/978-3-642-69024-2_43
  • Trifinopoulos J, Nguyen L-T, von Haeseler A, Minh BQ (2016) W-IQ-TREE: a fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Research 44: W232–W235. https://doi.org/10.1093/nar/gkw256
  • Turan C (1999) A note on the examination of morphometric differentiation among fish populations: The Truss System. Turkish Journal of Zoology 23: 259–263.
  • Ulber T (1993) Bemerkungen über cyrtodactyline Geckos aus Thailand nebst Beschreibungen von zwei neuen Arten (Reptilia: Gekkonidae). Mitteilungen aus dem Zoologischen Museum in Berlin 69: 187–200. https://doi.org/10.1002/mmnz.19930690202
  • Welch KRG, Cooke PS, Wright AS (1990) Lizards of the Orient: A Checklist. Robert E. Krieger Publishing Company, Malabar, Florida, 372 pp.
  • Wood Jr PL, Grismer LL, Muin MA, Anuar S, Oaks JR, Sites Jr JW (2020) A new potentially endangered limestone-associated Bent-toed Gecko of the Cyrtodactylus pulchellus (Squamata: Gekkonidae) complex from northern Peninsular Malaysia. Zootaxa 4751: 437–460.