﻿A new species of the Cyrtodactyluschauquangensis species group (Squamata, Gekkonidae) from Lao Cai Province, Vietnam

﻿Abstract We describe a new species of the genus Cyrtodactylus based on five adult specimens from Bac Ha District, Lao Cai Province, northern Vietnam. Cyrtodactyluslucisp. nov. is distinguished from the remaining Indochinese bent-toed geckos by a combination of the following morphological characteristics: medium size (SVL up to 89.5 mm); dorsal tubercles in 17–19 irregular transverse rows; ventral scales in 32–34 longitudinal rows at midbody; precloacal pores present in both sexes, 9 or 10 in males, 8 or 9 in females; 12–15 enlarged femoral scales on each thigh; femoral pores 9–12 in males, 5–10 in females; postcloacal tubercles 2–4; lamellae under toe IV 21–23; dorsal pattern consisting of 5 or 6 irregular dark bands, a thin neckband without V-shape or triangle shape in the middle, top of head with dark brown blotches; subcaudal scales transversely enlarged. Molecular phylogenetic analyses recovered the new species as the sister taxon to C.gulinqingensis from Yunnan Province, China, with strong support from all analyses and the two taxa are separated by approximately 8.87–9.22% genetic divergence based on a fragment of the mitochondrial ND2 gene. This is the first representative of Cyrtodactylus known from Lao Cai Province.


Introduction
The Cyrtodactylus chauquangensis species group is broadly distributed in the northern Indochina-Burma region, from northern Thailand and Laos to north central and northwestern Vietnam and to southwestern China (Uetz et al. 2023).Taxa within the group are almost exclusively adapted to karst ecosystems.Le et al. (2016) suggested that the group included at least ten species.Grismer et al. (2021aGrismer et al. ( , 2021b) ) provided a taxonomic review and analyzed phylogenetic relationships of 17 species and one undescribed form from northern Thailand.The group currently contains 23 recognized species with several taxa recently discovered from Yunnan Province, southern China (Grismer et al. 2021a(Grismer et al. , 2021b(Grismer et al. , 2021c;;Liu andRao 2021, 2022).
During our recent field trip in northern Vietnam, we collected five specimens of an unnamed gekkonid species from Bac Ha District, Lao Cai Province, which can be assigned to the Cyrtodactylus chauquangensis group based on molecular data.However, the population from Lao Cai Province can be distinguished from congeners by morphological differences and genetic divergence.Therefore, we describe it as a new species in the following.

Sampling
Field surveys were conducted in Bac Ha District, Lao Cai Province, Vietnam in June 2022 and October 2023 (Fig. 1).After being photographed in life, specimens were anesthetized and euthanized in a closed vessel with a piece of cotton wool containing ethyl acetate (Simmons 2002), fixed in 85% ethanol and subsequently stored in 70% ethanol.Specimens were subsequently deposited in the collections of the Institute of Ecology and Biological Resources (IEBR), Hanoi, Vietnam.

Molecular data and phylogenetic analyses
DNA was extracted using DNeasy Blood and Tissue kit (Qiagen, Germany) following manufacturer's instructions.Extracted DNA was amplified by HotStar Taq Mastermix (Qiagen, Germany) with 21 µl volume (10 µl of mastermix, 5 µl of water, 2 µl of each primer at 10 pmol and 2 µl of DNA).PCR conditions were: 95 °C for 15 min to active the taq; with 40 cycles at 95 °C for 30 s, 52 °C for 45 s, 72 °C for 60 s; and the final extension at 72 °C for 6 min.A fragment of the mitochondrial gene, NADH dehydrogenase subunit 2 (ND2), was amplified using the primer pair MetF1 (5'-AAGCTTTCGGGCCCATACC-3') and COIR1 (5'-AGRGTGCCAATGTCTTTGTGRTT-3') (Arevalo et al. 1994;Macey et al. 1997).PCR products were visualized using electrophoresis through a 2% agarose gel stained with ethidium bromide.Successful amplifications were purified to eliminate PCR components using GeneJET TM PCR Purification kit (ThermoFischer Scientific, Lithuania).Purified PCR products were sent to FirstBase (Malaysia) for sequencing in both directions.We included two samples of the newly discovered population from Lao Cai Province, one of Cyrtodactylus bichnganae, one of C. bobrovi, one of C. cucphuongensis, one of C. huongsonensis, one of C. ngoiensis, one of C. sonlaensis, one of C. taybacensis, and one of C. vilaphongi along with all available GenBank sequences of these species and other members of the Cyrtodactylus chauquangensis group.Two species, C. hontreensis and C. septimontium, of the C. intermedius group, were selected as outgroups (Grismer et al. 2021b).In the end, we were able to incorporate all ingroup taxa (Table 1).
After sequences were aligned by Clustal X v.2. 1 (Thompson et al. 1997), data were analyzed using maximum likelihood (ML) as implemented in IQ-TREE (Nguyen et al. 2015), maximum parsimony (MP) implemented in PAUP*4.0b10(Swofford 2001) and Bayesian inference (BI) as implemented in MrBayes v.3.2.7   (Ronquist et al. 2012).For the MP analysis, heuristic analysis was conducted with 100 random taxon addition replicates using tree-bisection and reconnection (TBR) branch-swapping algorithm, with no upper limit set for the maximum number of trees saved.Bootstrap support (BP) was calculated using 1000 pseudo-replicates and 100 random taxon addition replicates.All characters were equally weighted and unordered.For the ML analysis, we used IQ-TREE v. 1.6.8 (Nguyen et al. 2015) with a single model and 10000 ultrafast bootstrap replications (UFB).The optimal model for nucleotide evolution was determined using jModelTest v. 1.2.4 (Darriba et al. 2012).
For the BI analysis, we used the optimal model determined by jModelTest with parameters estimated by MrBayes v.3.2.7.Two independent analyses with four Markov chains (one cold and three heated) were run simultaneously for 10 7 generations with a random starting tree and sampled every 1000 generations.Loglikelihood scores of sample points were plotted against generation time to detect stationarity of the Markov chains.Trees generated prior to stationarity were removed from the final analyses using the burn-in function.The posterior probability values (PP) for all nodes in the final majority rule consensus tree were provided.We regard BP ≥ 70% and UFB and PP of ≥ 95% as strong support and values of < 70% and < 95%, respectively, as weak support (Hillis and Bull 1993;Ronquist et al. 2012;Minh et al. 2013).
The optimal model for nucleotide evolution was set to GTR+I+G for ML and BI analysis.The cut-off point for the burn-in function was set to 60, or 0.6% of the total number of trees generated, in the Bayesian analysis, as -lnL scores reached stationarity after 60,000 generations in both runs.Uncorrected pairwise divergences were calculated in PAUP*4.0b10.

Morphological characters
Measurements were taken with a digital calliper to the nearest 0.1 mm.Abbreviations are as follows: SVL: snout-vent length, measured from tip of snout to vent; TaL: tail length, measured from vent to tip of tail (* = regenerated); HL: head length, measured from tip of snout to retroarticular process of jaw; HW: head width, maximum width of head; HH: head height, from occiput to underside of jaws; OrbD: orbital diameter, greatest diameter of orbit; SE: snout to eye distance, from tip of snout to anterior-most point of eye; EE: eye to ear distance, from anterior edge of ear opening to posterior corner of eye; NE: nares to eye distance, from anterior-most point of eye to posterior-most point of nostril; ED: ear length, longest dimension of ear; ForeaL: forearm length, from base of palm to tip of elbow; CrusL: crus length, from base of heel to knee; TrunkL: trunk length, distance from axilla to groin measured from posterior edge of forelimb insertion to anterior edge of hindlimb insertion; BW: body width, the widest distance of body; Internar: internarial distance, distance between nares; Interorb: interorbital distance, shortest distance between left and right supraciliary scale rows.
Scale counts were taken as follows: SL: supralabials, counted from the first labial scale to corner of mouth; IL: infralabials, counted from the first labial scale to corner of mouth; N: nasal scales surrounding nare; IN: postrostrals or internasals; PM: postmentals; GST: granular scales surrounding dorsal tubercles; V: ventral scales in longitudinal rows at midbody; SLB: number of scales along the midbody from mental to anterior edge of cloaca; FP: femoral pores; PP: precloacal pores; PAT: postcloacal tubercles; TubR: tubercle, number of dorsal longitudinal rows of tubercles at midbody between the lateral folds; EFS: enlarged femoral scales, number of enlarged femoral scale beneath each thigh; NSF IV: number of subdigital lamellae on the fourth finger; NST IV: number of subdigital lamellae on the fourth toe.Bilateral scale counts were given as left/right; above sea level (asl).

Multiple Factor Analysis (MFA)
The MFA was also applied in this study using morphometric and meristic characteristics, including SVL, HL, HW, HH, OrbD, SE, EE, ED, ForeaL, CrusL, TrunkL, Internar, Interob and SL, IL, GST, V, TubR, EFS, FP, PP, PAT, NSF IV, NST IV.Other morphological characteristics were not used due to the limitation of available morphometric and meristic data or incomplete sampling (regenerated tail).All statistical analyses were performed using R Core Team (2023).The MFA used six quantitative groups -"SVL", "Head" (including HL, HW, HH), "Eye" (consist of OrbD, SE, EE, ED), "FT" (including ForeaL and CrusL), "TrunkL", "Inter" (consist of Internar and Interorb) and eight qualitative groups -"SpeciesInfor" (including Name of species and ID), "SL-IL" (consist of SL and IL in both sides), "GST_PAT_ TubR" (including GST, PAT in both sides and TubR), "V", "EFS" in both sides, "FP" in both sides, "PP", "LIV" (consist of NSF IV and NST IV in left side).To remove the effects of allometry, morphometric data were also normalized to adjust raw data of morphometrics through the allom() function in R package GroupStruct (available at heep://github.com/chankinonn/GroupStruct).Accordingly, the allometric formula is X adj = log 10 (X) -ß[log 10 (SVL)-log 10 (SVL mean )], where X adj = adjusted value; X = measured value; ß = unstandardized regression coefficient for each population and SVL mean = overall average SVL of two populations (Thorpe 1975(Thorpe , 1983;;Turan 1999;Lleonart et al. 2000;Grismer et al. 2021a;Chan and Grismer 2022).The ordination test was performed using packages Factoextra (Kassambara and Mundt 2017) and FactoMineR (Le et al. 2008) in the software R. The approach was applied to identify active groups and to explain phenotypic variance by estimating the first two Dim values-eigenvalue proportions.Similar coded colors in the MFA scatter plot, surrounded with convex hulls, were presented to visualize the phenotypic spaces of the new species and the most closely related species from China, namely Cyrtodactylus gulinqingensis Liu, Li, Hou, Orlov & Ananjeva, 2021; spaces were shown within a spatial coordinate of dimension axes (Dim1 and Dim2).To evaluate the overlap, the loadings of Dim1 and Dim2 of each Cyrtodactylus individual were extracted to identify the difference between the two species using the T-test.For all the tests, we applied a significance level of p < 0.05.

Phylogenetic analysis
The matrix of molecular data contained 1300 aligned characters, of which 580 were parsimony informative.The MP analysis produced a single most parsimonious tree (tree length = 2359, consistency index = 0.49, retention index = 0.66).Tree topologies from three analyses, ML, MP, and BI were similar and the Cyrtodactylus from Bac Ha District, Lao Cai Province was recovered with strong statistical support in all analyses as the sister taxon to C. gulinqingensis (BP = 94%; UBP = 100%; PP = 1.00) (Fig. 2).In terms of genetic divergences, the new species is separated from C. gulinqingensis by 8.87-9.22%based on a fragment of the mitochondrial ND2 gene.Genetically, it is also significantly divergent from other species within the C. chauquangensis group with a pairwise divergence of 12.32-23.85%(Suppl.material 1).

Morphological analysis
Morphologically, the new species from Bac Ha District, Lao Cai Province is closely similar to C. gulinqingensis from Yunnan Province, China, however, they plotted separately from each other in MFA (Fig. 3A) and there was a significant difference between two species (p < 0.05).The MFA also identified the data set of SVL, Head, Eye, FT, TrunkL, Inter, SL-IL, GST_PAT_TubR, V, EFS, FP, PP as active groups (Fig. 3B).The Eye, FT, Head, Inter, SVL and Trunk groups were the most important in both the first and second multi-factorial dimensions (Fig. 3C, D).Diagnosis.The new species can be distinguished from other members of the genus Cyrtodactylus by a combination of the following characteristics: Size medium (SVL up to 89.5 mm); dorsal tubercles in 17-19 irregular transverse rows; ventral scales in 32-34 longitudinal rows at midbody; precloacal pores present in both sexual, 9 or 10 in males, 8 or 9 in females; 12-15 enlarged femoral scales on each thigh; femoral pores 9-12 in males, 5-10 in females; postcloacal tubercles 2-4; lamellae under toe IV 21-23; dorsal pattern consisting of 5 or 6 irregular dark bands, a discontinuous thin neckband without V-shape or triangle shape in the middle, dorsal head surface with dark brown blotches; subcaudal scales transversely enlarged.

Taxonomy
Description of holotype.Adult male, snout-vent length (SVL) 86.3 mm; body relatively short (TrunkL/SVL 0.4); head distinct from neck, moderately long (HL/ SVL 0.28), relatively wide (HW/HL 0.69), slightly depressed (HH/HL 0.41); eye slightly large (OrbD/HL 0.24), pupils vertical; upper eyelid fringe with spinous scales; ear opening below the postocular stripes, obliquely directed and oval, small in size (ED/HL 0.06); two enlarged supranasals, separated from each other anteriorly by one internasal; nares oval, surrounded by supranasal, ros- tral, first supralabial and three postnasals; loreal region and frontal concave; snout long (SE/HL 0.41), round anteriorly, longer than diameter of orbit (OrbD/ SE 0.58); snout scales small, round, granular, larger than those in frontal and parietal regions; rostral wider than high with a medial suture, bordered by first supralabial on each side, nostrils, two supranasals and one internasal; mental triangular, wider than high; postmentals two, enlarged, in contact posteriorly, bordered by mental anteriorly, first infralabial laterally, and an enlarged chin scale posteriorly; supralabials 11/10; infralabials 11/10.Dorsal scales granular; dorsal tubercles round, keeled, conical, four or five times larger than the size of adjoining scales, each surrounded by 10 granular scales, tubercles forming 17 irregular longitudinal rows at midbody; ventral scales smooth, medial scales 2-3 times larger than dorsal granules, round, subimbricate, largest posteriorly, in 32 longitudinal rows at midbody; lateral folds present, without interspersed tubercles; gular region with homogeneous smooth scales; ventral scales between mental and cloacal slit 170; precloacal groove absent; three rows of enlarged scales present in posterior region of pore-bearing scales; ten precloacal pores arranged in a chevron; 12 or 13 enlarged femoral scales beneath thighs continuous with pore-bearing precloacal scales; femoral pores present on each enlarged femoral scales (except one on right thigh), 24 in total; precloacal pores large, horizontal elongated, positioned in posterior margin of scales; femoral pores small, round, positioned in the center of scales.
Coloration in life.Ground color of dorsal surface of head, neck, body, limbs and tail light brown.Dorsal surface of head with some dark brown blotches; labial region brown with yellowish cream stripes; skin above the eye gray; eyelid with light yellow color; iris yellow copper with black marking; pupil vertical, elliptical, black; nuchal loop dark brown, discontinous, extending from posterior corner of eye to the neck; tubercles on head, limbs, dorsum light brown to yellow; dorsum with five irregularly-shaped transversal bands and additional irregular smaller blotches; upper surface of limbs with irregular brown marks; six dark brown irregular bands on original part of tail while regenerated part of tail dark gray; chin, throat, chest, belly, lower limbs and ventral surface of tail cream.
Coloration in preservative.The overall color scheme slightly fades in 70% alcohol; yellow color disappeared in preservation while main characteristics are still clearly discernible; dorsal ground color of head, neck, body, limbs and tail grayish brown; color of chin, throat, chest, belly and lower limbs did not change noticeably in preservation.
Sexual dimorphism and variation.The males differ from females in the shape of precloacal pores (larger in males), and the presence of hemipenial swellings at the tail base.For other morphological characteristics see Table 2, Figs 4, 5.
Distribution.Cyrtodactylus luci sp.nov. is currently known only from the type locality in Bac Ha District, Lao Cai Province, Vietnam (Fig. 1).
Etymology.The species was named after the zoologist from the Vietnam National Museum of Nature, Vietnam Academy of Science and Technology, late Associate Professor Doctor Luc Van Pham, who contributed greatly to the biodiversity study in Vietnam.For the common names, we suggest Luc's Bent-toed Gecko (English) and Thạch sùng ngón lực (Vietnamese).Natural history.The bent-toed geckos were collected between 19:00 and 22:00, both on limestone cliffs and on trees, about 1.0-1.8m above the ground.The surrounding habitat was secondary karst forest of medium and small hardwoods mixed with shrubs and vines (Fig. 6).Air temperature was 25.9 °C and relative humidity was 92%.
Comparisons.Cyrtodactylus luci sp.nov. is distinguishable from all other members of the C. chauquangensis species group by a unique combination of morphological characteristics.

Discussion
The new species from Bac Ha District, Lao Cai Province, is most similar to Cyrtodactylus gulinqingensis, a recently described species from Muguan County, Wenshan Prefecture, Yunnan Province of China (Liu et al. 2021).In terms of geographic distribution, the type locality of C. luci is approximately 40 km distant from that of its sister species in China.However, they are distinguished from each other by morphological differences as well as a genetic divergence of 8.87-9.22%(ND2 gene).
Our tree topology (Fig. 2) is similar to that reported in Grismer et al. (2021b).However, while C. auribalteatus is recovered as a member of the clade including C. dumnuii, C. wayakonei and other taxa in this study, it is grouped with the lineage consisting of C. sonlaensis, C. huongsonensis and C. soni in Grismer et al. (2021b).According to our phylogenetic analyses, the new species and C. gulinqingensis from Yunnan cluster with the latter clade with strong nodal support provided only by BI (Fig. 2).In addition to C. luci and C. gulinqingensis, the other species in the group occur in Son La (C.sonlaensis) and Ninh Binh (C.soni) provinces and the suburb of Ha Noi City (C.huongsonensis), northwestern Vietnam.
In the Cyrtodactylus chauquangensis group, except for C. doisuthep, a species known from dry evergreen and deciduous dipterocarp forests in Thailand (Kunya et al. 2014), all 23 remaining species are karst dwellers, comprising three species from Yunnan Province of China, five species from northern Laos, four species from northern Thailand, and 12 species from northern Vietnam (Uetz et al. 2023, this study).In terms of altitudinal distribution range, the members of this species group are found at elevations from 17 m (C.soni) to 1660 m (C.doisuthep) but most of them occur at elevations between 300 and 800 m a.s.l (Kunya et al. 2015;Le et al. 2016).The new species is the 24 th species of the C. chauquangensis group, the first species from Lao Cai Province and the eastern side of the Red River in Vietnam, and the 53 rd species of Cyrtodactylus known from Vietnam (Ngo et al. 2022;Uetz et al. 2023).

Figure 2 .
Figure 2. Phylogram based on the Bayesian analysis.Number above and below branches are ML/MP bootstrap and ultrafast bootstrap values and Bayesian posterior probabilities (≥ 50%), respectively.Asterisk and hyphen denote 100% and > 50% values, respectively.

Figure 3 .
Figure 3.A MFA of Cyrtodactylus luci sp.nov.from Vietnam and C. gulinqingensis from China B scatterplot the groups of all variables for Dim1 and Dim2 axes in the MFA, green triangles as inactive groups of variables, red triangles as active groups of variables C bar plot of groups' contribution to the first axes (Dim1) in the MFA D bar plot of groups' contribution to the second axes (Dim2) in the MFA.

Table 1 .
Species of Cyrtodactylus used in the phylogenetic analysis including localities and GenBank accession numbers of the mitochondrial NADH dehydrogenase subunit 2 (ND2) fragment gene (-: data unavailable).