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). Molecular phylogenies recovered multiple monophyletic mitochondrial lineages of individuals (populations) that were 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. In this way, delimiting (phylogeny) and diagnosing (taxonomy) species are not conflated (Frost and Hillis 1990; Frost and Kluge 1994; Hillis 2019).

Field surveys were conducted during the month of August in 2017 and 2022 in Lung Cu Commune, Dong Van District, Ha Giang Province, northeastern Vietnam. After having been photographed alive, newly collected individuals, were anaesthetized and euthanized in a closed vessel with a piece of cotton wool containing ethyl acetate (Simmons 2002), fixed in 85% ethanol, thereafter and subsequently transferred to 70% ethanol for permanent storage. Tissue samples were preserved separately in 70% ethanol. Specimens were deposited in the collection of the Vietnam National University of Forestry (VNUF) and the Institute of Ecology and Biological Resources (IEBR), Hanoi, Vietnam.

Terminology of morphological characters follows Zug (2010), Nguyen et al. (2013), and Grismer et al. (2021). Measurements (morphometrics) were taken with a digital caliper to the nearest 0.1 mm on the left side of the body. Character abbreviations and descriptions are as follows: SVL: snout-vent length (from the tip of snout to the vent); TaL: tail length (from the vent to the tip of the tail, original and regenerated); TrunkL: trunk length (from the posterior margin of the forelimb at its insertion point on the body to the anterior margin of the hindlimb insertion point on the body); HeadL: head length (from the posterior margin of the retroarticular process of the lower jaw to the tip of the snout); HeadW: head width (measured at the angle of the jaws); EyeD: eye diameter (the greatest horizontal diameter of the eyeball); SnEye: snout-eye length (from anteriormost margin of the eyeball to the tip of snout); NarEye: nares-eye length (from the anterior margin of the eyeball to the posterior margin of the external nares); SnW: internarial width (measured between the nares across the rostrum); and EarD: ear diameter (maximum diameter of the ear opening).

The following meristic characters were counted using the same dissecting microscope: Chin: chin scales (the number of scales contacting the infralabials and mental from the juncture of the second and third infralabials on the left side to the juncture of the second and third infralabials on the right side); CN: circumnasal scales (the number of scales abutting the external naris, exclusive of the rostral and first supralabial); SnS: the number of scales between the supranasals; SL: supralabial scales (number of enlarged scales bordering the mouth on the upper jaw from the rostral to a point in line with the posterior margin of the orbit); IL: infralabial scales (the number of enlarged scales bordering the mouth on the lower jaw from the mental to a point in line with the posterior margin of the orbit); VS: ventral scales (the number of longitudinal ventral scales at mid-body contained within one eye diameter); DS: dorsal scales (the number of longitudinal dorsal scales at mid-body contained within one eye diameter); lamellae formulae determined as the number of U-shaped, subdigital lamellae (split and single) on the digital pads on digits II–V of the hands and feet (on the left); SL1F: the number of subdigital lamellae wider than long on the first finger; SL1T: the number of subdigital lamellae wider than long on the first toe; the total number of precloacal and femoral pores (i.e. the contiguous or discontinuous rows of femoral and precloacal pore-bearing scales); and CloacS: the number of cloacal spurs.

The categorical color pattern characters evaluated were the presence or absence of dark markings on the dorsum (BodPartn); presence or absence of a dark pre– and/ or postorbital stripe extending onto at least the neck (PosOStrp); the presence or absence of a linear series of white postorbital and dorsolateral spots on the trunk (PosOTrSpt); and the presence or absence of pale–colored, anteriorly projecting arms of the pale–colored postsacral marking (PosSacrlMakg).

We obtained 1,038 base pairs of the NADH dehydrogenase subunit 2 (ND2) gene from 45 specimens from GenBank and five newly sequenced specimens from Lung Cu Commune, Dong Van District, Ha Giang Province, northeastern Vietnam for the phylogenetic analyses (Table 1). Catalog numbers are abbreviated as follows: CSUFT, Central South University of Forestry and Technology, Changsha, China; ITBCZ, Institute of Tropical Biology, Zoological Collection; KIZ: Kunming Institute of Zoology; LSUHC, La Sierra University Herpetology Collection; NJUN, Nanjing Normal University, Nanjing; MVZ, Museum of Vertebrate Zoology; NUOL, National University of Laos; SYS, The Museum of Biology, Sun Yatsen University (SYS), Guangzhou; ZFMK, Zoologisches Forschunginstitut und Museum Alexander Koenig, Bonn, Germany. Hemiphyllodactylus harterti (Werner, 1900) of the H. harterti group was used to root the tree following Grismer et al. (2013). For molecular phylogenetic analyses, we used the protocols of Le et al. (2006) for DNA extraction, amplification, and sequencing. Tissue samples were extracted using DNeasy blood and tissue kit, Qiagen (California, USA). Extracted DNA from the fresh tissue was amplified by PCR mastermix (Fermentas, Canada). A fragment of the mitochondrial gene ND2 (NADH dehydrogenase subunit 2) was amplified using the primer pair ND2f101A (CAACAGAAGCCACAACAAAAT) and HemiR (GAAGAAGAGGCTTGGKAGGCT) (Nguyen et al. 2013). The PCR volume consisted of 21 μl (10 μl of mastermix, 5 μl of water, 2 μl of each primer at 10 pmol/μl and 2 μl of DNA or higher depending on the quantity of DNA in the final extraction solution). The following temperature profile for PCR was used: 95 °C for 5 min to activate the taq; with 40 cycles at 95 °C for 30 s, 50 °C for 45 s, 72 °C for 60 s; and the final extension at 72 °C for 6 min. PCR products were subjected to electrophoresis through a 1% agarose gel (UltraPureTM, Invitrogen). Gels were stained for 10 minutes in 1×TBE buffer at 2 pg/ml of ethidium-bromide and visualized under UV light. Successful amplifications were purified to eliminate PCR components using GeneJETTM PCR Purification kit (Fermentas, Canada). Purified PCR products were sent to 1^st Base (Selangor, Malaysia) for sequencing.

Maximum likelihood (ML) and Bayesian inference (BI) analyses were used to estimate the phylogenetic relationships among the sampled sequences in the alignment. An ML phylogeny was estimated using the IQ-TREE webserver (Nguyen et al. 2015; Trifinopoulos et al. 2016) preceded by the selection of substitution models using the Bayesian Information Criterion (BIC) in Model Finder (Kalyaanamoorthy et al. 2017), which supported TN+F+G4 as the best fit model of evolution for ND2 codon position 1, TN+F+G4 for position 2 and TIM2+F + G4 for position 3. One thousand bootstrap pseudoreplicates via the ultrafast bootstrap (UFB; Hoang et al. 2018) approximation algorithm were employed and nodes having ML UFB values of 95 and above were considered highly supported (Minh et al. 2013).

The Bayesian inference (BI) analysis was carried out using Mr Bayes 3.2.7. (Ronquist et al. 2012) XSEDE on the CIPRES Science Gateway (Cyberinfrastructure for Phylogenetic Research; Miller et al. 2010), employing the GTR+I+G model of evolution to all partitions. Two independent Markov chain Monte Carlo (MCMC) simulations were performed each with four chains, three hot and one cold. The MCMC simulation was run for 100 million generations, sampled every 10,000 generations, and the first 10% of each run was discarded as burn-in. Convergence and stationarity of all parameters from both runs were checked using Tracer v1.6 (Rambaut et al. 2014) to ensure effective sample sizes (ESS) were above 200. Post-burn-in sampled trees from both runs were combined using the sumt function in MrBayes and a 50% majority-rule consensus tree was constructed. Nodes with Bayesian posterior probabilities (BPP) of 0.95 and above were considered highly supported (Huelsenbeck et al. 2001; Wilcox et al. 2002). Uncorrected pairwise sequence divergences (p-distance) among and within species were calculated in MEGA11 (Kumar et al. 2016) (Table 2).

We compared the Lung Cu population to its closest relatives Hemiphyllodactylus huishuiensis and H. yanshanensis from southern China and to other species in clade 6. Hemiphyllodactylus dushanensis was not evaluated due to the absence of comparative morphological data. Morphological data used for the statistical analyses were obtained from the Lung Cu population (seven specimens) and data from 76 specimens of 12 other species available from previous studies (Ngo et al. 2014; Nguyen et al. 2014; Yan et al. 2016; Sung et al. 2018; Eliades et al. 2019; Nguyen et al. 2020; Do et al. 2020; Zhang et al. 2020; Agung et al. 2022). Raw data for H. typus (four specimens [cat. Nos. CUMZ. R 2013.19; RBWS 19309; JAV bg. R.39; MS92002]) was obtained from the La Sierra University Herpetological Collection, La Sierra University, Riverside, California, USA. Museum Abbreviations: CUMZ - Museum of Zoology, Chulalongkorn University, Bangkok, Thailan; JAV-J.-M. VINSON; MS-Montri Sumontha; RBWS: Royal Belgian Institute of Natural Sciences. The raw morphological data are provided in Tables 3, 4.

All morphological analyses were conducted in R v.4.2.1 (R Core Team, 2021). An ANOVA was conducted on morphometric and meristric characters with statistically similar variances (i.e. p values > 0.05 in a Levene’s test) to search for the presence of statistically significant mean differences (p < 0.05) across the data set. Characters bearing statistically significant differences were subjected to a TukeyHSD test to ascertain which population pairs differed significantly from each other for those particular characters. Boxplots were generated in order to visualize the range, mean, and degree of differences among species bearing statistically different mean values for sets of characters.

Morphometric characters used in the statistical analyses were SVL, TrunkL, HeadL, HeadW, SnEye, NarEye, EyeD, and SnW. Tail metrics were not used due to the high degree incomplete sampling (i.e., regenerated, broken, or missing). To remove potential effects of allometry on morphometric traits (sec. Chan and Grismer 2022), we used 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 into a single data set so as not to conflate potential intra- with interspecific variation (Reist 1986; McCoy et al. 2006) (Table 3). Meristic characters (scale counts) used in statistical analyses were Chin, CN, SL, IL, VS, DS, SL1F, and SL1T (Table 4). Precloacal and femoral pores were excluded from the multivariate analyses due to their absence in females. Categorical characters analyzed were the presence or absence of BodPartn; PosOStrp; PosOTrSpt; and PosSacrlMakg (Table 3).

Morphospatial relationships based on the normalized morphometrics and meristics were visualized using principal component analysis (PCA) from the ADEGENET package in R (Jombart et al. 2010) to determine if their positioning was consistent with the putative species boundaries delimited by the molecular phylogenetic analyses and defined by the univariate analyses (see above). PCA, implemented using the “prcomp()” command in R, is an indiscriminate analysis plotting the overall variation among individuals (i.e. data points) while treating each individual independently (i.e. not coercing data points into pre-defined groups). Subsequent to the PCA, a discriminant analysis of principle components (DAPC) was used to test for corroboration and further discrimination of morphospatial differences among the putative species. DAPC a priori groups the individuals of each predefined population inferred from the phylogeny into separate clusters (i.e. plots of points) bearing the smallest within-group variance that produce linear combinations of centroids having the greatest between-group variance (i.e. linear distance; Jombart et al. 2010). DAPC relies on standardized data from its own PCA as a prior step to ensure that variables analyzed are not correlated and number fewer than the sample size. Principal components with eigenvalues accounting for 90–95% of the variation in the data set were retained for the DAPC analysis according to the criterion of Jombart et al. (2010).

To test and further corroborate the PCA and DAPC analyses, we conducted a multiple factor analysis (MFA) on the above mentioned morphological characters plus the categorical color pattern charactrers for a near total evidence data set (see Tables 5, 6). The MFA was implemented using the mfa () command 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. In an MFA, each individual is described by a different set of variables (i.e. characters) which are structured into different data groups in a global data frame–in this case, quantitative data (i.e. meristics and normalized morphometrics) and categorical data (i.e. color pattern). In the first phase of the analysis, separate multivariate analyses are carried out for each set of variables–principal component analyses (PCA) for the quantitative data sets and a multiple correspondence analysis (MCA) for categorical data. The data sets are then normalized separately by dividing all their elements by the square root of their first eigenvalues. For the second phase of the analysis, the normalized data sets are concatenated into a single matrix for a global PCA. Standardizing the data in this manner prevents one data type from overleveraging another. In other words, the normalization of the data in the first phase prevents data types with the highest number of characters or the greatest amount of variation from outweighing other data types 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 the contribution of each data type and character to the overall variation in the data set (Pagès 2015; Kassambara and Mundt 2017).

The ML and BI analyses returned well-supported, nearly identical topologies for all the known species of clade 6 (Agung et al. 2022) and were concordant with genus-wide phylogenies from previous analyses (Grismer et al. 2013; Grismer et al. 2020; Agung et al. 2022). The specimens from Lung Cu Commune, Dong Van District, Ha Giang Province, northeastern Vietnam, formed a monophyletic lineage in the typus group (Grismer et al. 2013) embedded within a larger lineage of Indochinese species where it is nested within clade 6. Herein, it was recovered as the well-supported (BI 0.99/ML 91) sister species to a lineage comprised of H. huishuiensis and H. yanshanensis (Fig. 2). The Lung Cu lineage differs from H. huishuiensis and H. yanshanensis by having uncorrected average pairwise sequence divergences of 5.3% and 6.8%, respectively (Table 2). In terms of genetic distance, the Lung Cu population is similar to H. nahangensis and H. dushanensis, but it has a genetic distance of 4.6% and 4.8% from H. nahangensis and H. dushanensis, respectively. Additionally, significant differences of morphometric and meristic traits among them are indicated in the comparisons. As such we hypothesize this population represents a new evolutionary lineage and is described as a new species below.

Figure 2. 

Maximum likelihood consensus tree of clade 6 (excepting Hemiphyllodactylus harterti used to root the tree following Grismer et al. 2013) of Hemiphyllodactylus with Bayesian posterior probabilities (BPP) and Ultrafast Bootstrap support (UFB) value coding at the nodes, respectively.

The ANOVA and TukeyHSD post hoc tests of the adjusted morphometric and meristic characters were consistent with the phylogenetic results and their high degree of pairwise genetic distance among the species in recovering a number of statistically significant mean differences between the Lung Cu population and all other species (Tables 5, 6). Variation in all morphometric and metric characters are visualized in Figs 3, 4.

Figure 3. 

Boxplot comparisons of significantly different meristic characters between the Lung Cu population and other Hemiphyllodactylus species of clade 6. Pale blue circles are means and the black horizontal bars are medians.

Figure 4. 

Violin plots of the significantly different morphometric characters between the Lung Cu population and other Hemiphyllodactylus species of clade 6 overlain with box plots showing the range, frequency, mean (white dot), and 50% quartile (black rectangle) of the size-adjusted morphometric and meristic characters. New species in bold.

The clustering of species in the PCA of clade 6 showed the first two principal components (PC1 and PC2) recovered 52.5% of the variation in the normalized morphometric and meristic data set (Fig. 5A). PC1 accounted for 34.6% of the variation in the data set and loaded most heavily for trunk length (TrunkL), head width (HeadW), snout-eye length (SnEye), nares-eye length (NarEye), and eye diameter (EyeD). PC2 accounted for an additional 17.9% of the data set and loaded most heavily for head length (HeadL), circumnasal scales (CN), dorsal scales (DS), the number of subdigital lamellae wider than long on the first finger (SL1F), and the number of subdigital lamellae wider than long on the first toe (SL1T) along PC2 (Fig. 5A, Table 7). The PCA recovered Hemiphyllodactylus lungcuensis sp. nov. to be widely separated from most other species and only overlapping with the distantly related H. nahangensis and minimally overlapping with H. dupanglingensis. The new species is well-separated from most other species in the DAPC and only slightly overlapping with the distantly related H. dupanglingensis (Fig. 5B). More importantly, its 66% confidence ellipse does not overlap with that of any other species.

Figure 5. 

A principal component analysis (PCA) of Hemiphyllodactylus species of clade 6 showing their morphospatial relationships along the first two principal components B discriminant analysis of principal components (DAPC) based on retention of the first five PCs accounting for 98% of the variation.

The MFA analysis recovered all species as being separated from one another including Hemiphyllodactylus lungcuensis sp. nov. and H. nahangensis (Fig. 6A). The normalized morphometric data contributed to approximately 50% of the variation along Dim-1 followed by the meristic and categorical data. For Dim-2, the categorical data contributed 60% of the variation followed by meristic and mensural data. Dim-3 showed that meristic data contributed 45% of the variation followed by categorical and mensural data. Finally, Dim-4 showed that meristic data contributed 70% of the variation followed by categorical and mensural data (Fig. 6B).

Figure 6. 

A MFA scatter plot showing the morphospatial relationships among the Hemiphyllodactylus species of clade 6 B bar graphs showing the percent contribution of each data type to the overall variation in the data set for the first four dimensions. The dashed red line in the bar graphs indicates the expected average value if the contributions of each data type were equal.

The authors have declared that no competing interests exist.

No ethical statement was reported.

This research is funded by the Vietnam Academy of Science and Technology (Grant Number CSCL09.02/22-23). Research of Thuong Huyen Nguyen was funded by the Master, PhD Scholarship Programme of Vingroup Innovation Foundation (VINIF), code VINIF.2022.TS125. Research of Vinh Quang Luu in Herpetology Laboratory, Department of Biology, La Sierra University, U.S. was supported by the Fulbright Program.

All authors contributed to this work. Methodology: L. Lee Grismer, Vinh Quang Luu, Jess L. Grismer; collecting data: Thuong Huyen Nguyen, Quyen Hanh Do, Vinh Quang Luu, Tuoi Thi Hoang; data analysis: Vinh Quang Luu, Tuoi Thi Hoang, L. Lee Grismer, Jesse L. Grismer, Minh Duc Le; writing—original draft preparation, Thuong Huyen Nguyen, Quyen Hanh Do, Vinh Quang Luu, Tuoi Thi Hoang; writing—review and editing: L. Lee Grismer, Minh Duc Le, Cuong The Pham, Thomas Ziegler, Truong Quang Nguyen, Vinh Quang Luu; supervision: L. Lee Grismer, Thomas Ziegler, Truong Quang Nguyen.

Vinh Quang Luu https://orcid.org/0000-0002-0634-1338

Thuong Huyen Nguyen https://orcid.org/0000-0002-4303-525X

Quyen Hanh Do https://orcid.org/0000-0002-9437-4673

Cuong The Pham https://orcid.org/0000-0001-5158-4526

Truong Quang Nguyen https://orcid.org/0000-0002-6601-0880

Minh Duc Le https://orcid.org/0000-0002-2953-2815

Thomas Ziegler https://orcid.org/0000-0002-4797-609X

Jesse L. Grismer https://orcid.org/0000-0002-2542-5430

L. Lee Grismer https://orcid.org/0000-0001-8422-3698

All of the data that support the findings of this study are available in the main text or Supplementary Information.