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).

Four samples of the newly collected specimens were analyzed, two from Gia Lai Province, Vietnam (VNUF R.2020.22 – field number GL.02, VNUF R.2020.33 – field number GL.03) and two from Vientiane Province, Laos (VNUF R.2022.42 – field number MF.02, VNUF R.2022.52 – field number MF.03). We used the protocols of Nguyen et al. (2021) for DNA extraction, amplification, and sequencing. The complete NADH dehydrogenase subunit 2 (ND2) gene with six partial or complete adjacent tRNAs, approximately 1200 bp long, respectively, were amplified and sequenced using the primer pair, MetF1(5’-AAGCTTTCGGGCCCATACC-3’) and COIR1(5’-AGRGTGCCAATGTCTTTGTGRTT-3’) (Macey et al. 1997). Genomic DNA was extracted from all liver tissues stored in ethanol following the standard protocols of DNeasy blood and tissue kit, Qiagen (California, USA). The PCR volume consisted of 20 µl (1 µl each primer, 7 µl water, 10 µl of Taq mastermix and 1 µl DNA template). PCR conditions were: 95 °C for 5 min, followed by 42 cycles: 95 °C for 30 s, 50 °C for 45 s and 72 °C for 60 s with a final elongation step for 6 min at 72 °C. PCR products were visualized using electrophoresis through a 1.2% agarose gel, marker 100 bp, 1X TAE and stained with RedSafe Nucleic Acid Staining Solution and photographed under UV light of Geldoc system (Quantum CX5, Villber, France). Successful amplifications were purified using innuPREP Gel Extraction Kit (Analytik Jena, Germany). Cleaned PCR products were sent to 1^st Base (Malaysia) for sequencing in both directions.

We obtained 1,444 base pairs of NADH dehydrogenase subunit 2 gene (ND2) sequence data and the flanking tRNAs from 29 ingroup individuals of Dixonius representing 13 nominal species including the new samples from Vietnam and Laos. Heteronotia spelea was used as an outgroup to root the tree based on the phylogenetic results generated by Gamble et al. (2015). Sequence data for other species were acquired from GenBank. Newly generated sequences were deposited in GenBank (Table 1).

Maximum likelihood (ML) and Bayesian inference (BI) were used to estimate phylogenetic trees. Best-fit models of evolution determined in IQ-TREE (Nguyen et al. 2015) using the Bayesian information criterion (BIC) implemented in ModelFinder (Kalyaanamoorthy et al. 2017) indicated that F81+F was the best-fit model of evolution for the tRNAMET and K2P+I, and HKY+F+G4 were the best models of evolution for tRNAs2 and ND2, respectively. 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 v. 3.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 80,000,000 generations with the cold chain sampled every 8000 generations and the first 10% of each run being discarded as burn-in. The posterior distribution of trees from each run was summarized using the sumt function in MrBayes v. 3.2.4 (Ronquist et al. 2012). Stationarity was checked with Tracer v. 1.6 (Rambaut et al. 2014) to ensure effective sample sizes (ESS) for all parameters were well above 200. We considered Bayesian posterior probabilities (BPP) of 0.95 and above and ultrafast bootstrap support values (UFB) of 95 and above as an indication of strong nodal support (Huelsenbeck et al. 2001; Minh et al. 2013). Uncorrected pairwise sequence divergences (p-distance) were calculated in MEGA 11 (Kumar et al. 2016) using the complete deletion option to remove gaps and missing data from the alignment prior to analysis.

A time-calibrated Bayesian phylogenetic tree was estimated using BEAST 2 (Bayesian Evolutionary Analysis by Sampling Trees) v. 2.7.3 (Drummond et al. 2012) implemented in CIPRES (Cyberinfrastructure for Phylogenetic Research; Miller et al. 2010) where the ingroup node subtending the split between Dixonius aaronbaueri and the remaining species was given a 24.04 mya prior with an offset range of 20.23–27.68 mya following Gamble et al. (2015). The split between Heteronotia and Dixonius was set at 45.0 mya with an offset range of 33.3–56.8 mya (Gamble et al. 2015). An input file was constructed in BEAUti (Bayesian Evolutionary Analysis Utility) v. 2.7.3. An optimized relaxed clock with unlinked site models, linked clock models and linked trees, and a calibrated Yule prior were employed for the species level. BEAST Model Test (Bouckaert and Drummond 2017), implemented in BEAST, was used to numerically integrate over the uncertainty of substitution models while simultaneously estimating phylogeny using Markov chain Monte Carlo (MCMC). MCMC chains were run for 80 million generations and logged every 8,000 generations. The BEAST log file was visualized in Tracer v. 1.7.2 (Rambaut et al. 2014) to ensure effective sample sizes (ESS) were above 200 for all parameters. A maximum clade credibility tree using mean heights at the nodes was generated using TreeAnnotator v. 2.7.3 (Rambaut and Drummond 2013) with a burn-in of the first 10% of each run. Nodes with Bayesian posterior probabilities (BPP) of 0.95 and above were considered strongly supported (Huelsenbeck et al. 2001; Wilcox et al. 2002).

The morphological data set comprised six closely related species including six type specimens of Dixonius minhlei from Dong Nai Province, Vietnam (IEBR A.0801-02, VNMN R.2016.1-2, ZFMK 97745-46), three type specimens of D. lao from Khammouane Province, Laos (VNUF R.2016.2, IEBR A.2016.5-6), eight specimens of D. siamensis from Pursat Province, Cambodia (LSUHC 07328, 07378, 08420, 08487, 08491, 08522, 09284, 09289), five type specimens of D. somchanhae from Vientiane Capital, Laos (VNUF R.2020.1-5), four specimens of D. sp. from Gia Lai Province, Vietnam, and 12 specimens of D. vietnamensis from Nha Trang Province, Vietnam (ZRC 2.6024-27, IEBR R.2016.1, 2016.3, 2016.4, VNMN R.2016.3-4, ZFMK 97747-49).

Morphological data included both meristic and morphometric characters. Morphological characters were taken from the 44 specimens following Bauer et al. (2004) and Ngo and Ziegler (2009). Morphometric characters were taken after preservation with a digital caliper to the nearest 0.1 mm under a zoom stereomicroscope on the right/left of the body. Recorded data included: SVL: snout-vent length (taken from the tip of the snout to the vent), TaL: tail length (taken from the vent to the tip of the tail, original or partially regenerated), TW: tail width (taken at the base of the tail immediately posterior to the postcloacal swelling), BW: body width (greatest width of torso, taken at the level of midbody), HL: head length (the distance from the posterior margin of the retroarticular process of the lower jaw to the tip of the snout), HW: head width (the distance from the posterior margin of the retroarticular process of the lower jaw to the tip of the snout), HD: head depth (the maximum height of head measured from the occiput to base of the lower jaw), EL: ear length (greatest oblique length across the auditory meatus), TBL: Tibia length (taken on the ventral surface from the posterior surface of the knee while flexed 90° to the base of the heel), AG: axilla to groin length (taken from the posterior margin of the forelimb at its insertion point on the body to the anterior margin of the hind limb at its insertion point on the body), FA: forearm length (taken on the ventral surface from the posterior margin of the elbow while flexed 90° to the inflection of the flexed wrist), ED: eye diameter (the greatest horizontal diameter of the eye-ball), EN: eye nostril distance (measured from the anterior margin of the bony orbit to the posterior margin of the external nares), ES: eye snout distance (measured from anteriormost margin of the bony orbit to the tip of snout), EE: eye ear distance (measured from the anterior edge of the ear opening to the posterior edge of the bony orbit), IN: internarial distance (measured between the external nares across the rostrum), IO: interorbital distance (measured between the dorsal-most edges of the bony orbits).

Meristic data taken were: V: ventral scales (counted transversely across the abdomen midway between limb insertions from one ventrolateral fold to the other), DTR: longitudinal rows of dorsal tubercles (counted transversely across the body midway between the limb insertions from one ventrolateral body fold to the other), PV: paravertebral scales (counted in a paravertebral row from first scale posterior to parietal scale to last scale at the level of vent opening), PV`: paravertebral scales (counted in a row between limb insertions), T4: lamellae under fourth toe (counted from the distal scale containing claw to basal scale that broadly contacts adjacent fragmented scales), IOS: Interorbital scales (counted at narrowest point between orbits), ICS: interciliary scales (counted between supraciliaries at midpoint of orbit), SPL: supralabials (counted from the largest scale at the corner of the mouth to the rostral scale), IFL: infralabials (counted from termination of enlarged scales at the corner of the mouth to the mental scale), MO: number of supralabial at midorbital position, PP: precloacal pores in males.

Color pattern on dorsum including the presence or absence of canthal stripes (CanthStrp), the presence or absence of strong darkly barred lips (LipBar), the presence or absence of dark-colored round blotches on the top of the head (RdHdBlch) and dorsum (RdBodBlch), and the presence or absence of two regularly arranged whitish tubercles on flanks (Tub). The raw morphological data for all characters and specimens are presented in Tables 2, 3.

All statistical analyses were performed using R v. 4.2.1 (R Core Team, 2021). Morphometric characters used in the statistical analyses were SVL, BW, HL, HW, HD, EL, ED, EN, ES, EE, IN, IO, FAr, TBLr, and AGr. 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). All data were scaled to their standard deviation to ensure they were analyzed on the basis of correlation and not covariance. Meristic characters (scale counts) used in statistical analyses were SPLr/l, IFLr/l, MO, IOS, ICS, V, DTR, and T4r/l. Precloacal and femoral pores were omitted from the analyses due to their absence in females. Categorical characters analyzed were CanthStrp, LipBar, RdHdBlch, RdBodBlch, and Tub.

A Levene’s test for normalized morphometric and meristic characters was conducted to test for equal variances across all groups. Analyses of variance (ANOVA) were conducted on meristic and normalized morphometric characters (see below) with statistically similar variances to search for the presence of statistically significant mean differences (p < 0.05) among species across the data set. Characters bearing statistical differences were subjected to a TukeyHSD test to ascertain which species pairs differed significantly from each other for those particular characters. Boxplots were generated for discrete meristic characters in order to visualize the range, mean, median, and degree of differences between pairs of species bearing statistically different mean values and violin plots were generated for continuous morphometric characters to visualize the same.

Morphospatial positions 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 differences 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 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 of the data. 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 each data type’s and each character’s contributions to the overall variation in the data set (Pagès 2015; Kassambara and Mundt 2017).

The results of ML, BI, and BEAST analyses produced trees with identical topologies and strong support at nearly every node (Figs 2, 3). The molecular analyses suggest that Dixonius aaronbaueri is the sister species to a clade containing all other species of Dixonius. Additionally, all analyses recovered the newly discovered population from Chu Se District, Gia Lai Province, Vietnam as the strongly supported (1.00/100) sister species of D. minhlei and the newly discovered population from Muangfuang District, Vientiane Province, Laos as the strongly supported (1.00/100) sister species of D. lao (Figs 2, 3). Uncorrected pairwise sequence divergences among Dixonius species ranged from 2.57–18.84% (Table 4). Ranges for the new species described (see below) are as follows: new species from Vietnam 3.60–15.73%, being most similar to D. minhlei and most distant to D. aaronbaueri and new species from Laos 3.10–18.17%, being most similar to D. lao and most distant to D. aaronbaueri.

Figure 2. 

Maximum likelihood topology of the Dixonius species from Vietnam and Laos with ultrafast bootstrap values (UFB) and Bayesian posterior probabilities (BPP) at the nodes, respectively.

Figure 3. 

BEAST chronogram of the Dixonius species from Vietnam and Laos. Numbers at the nodes are mean ages in millions of years. Bars represent 95% highest posterior densities.

The time-calibrated BEAST analysis places the divergence between Dixonius aaronbaueri and the remaining species of Dixonius at approximately 24.04 mya (20.23–27.68 highest posterior density [HPD]). Within the Vietnam’s lineages, D. gialaiensis sp. nov. and D. minhlei diverged from each other at approximately 3.19 mya (0.79–5.78 HPD) and within the Lao lineages, D. muangfuangensis sp. nov. and D. lao diverged approximately 3.47 mya (1.37–6.16 HPD) (Fig. 3).

The first two principal components (PC1 and PC2) of the PCA analysis recovered 56.6% of the variation in the morphometric and meristic data set (Fig. 4A) and loaded most heavily for body width (BW), head width (HW), eye nostril distance (EN), eye snout distance (ES), and eye ear distance (EE) along PC1 and interorbital distance (IO), supralabials (SPLr/l), number of supralabial at midorbital position (MO), and ventral scales (V) along PC2 (Table 5). The PCA recovered D. gialaiensis sp. nov. and D. muangfuangensis sp. nov. to be widely separated from most other species with D. muangfuangensis sp. nov. only overlapping with the distantly related D. siamensis. The two distantly related new species are well-separated from most other species in the DAPC but each overlaps with one other species in their 67% inertia ellipses (Fig. 4B).

Figure 4. 

A principal component analysis (PCA) of Dixonius species showing their morphospatial relationships along the first two components based on normalized morphometric and meristic characters B discriminant analysis of principal components (DAPC) based on retention of the first five PCs with 67% inertia ellipsoids.

The MFA analysis recovered all species to be separated from one another including Dixonius muangfuangensis sp. nov. and D. siamensis (Fig. 5A). The morphometric data contributed to approximately 40% of the variation along Dim-1 followed by the categorical and meristic data. For Dim-2, the categorical data contributed 80% of the variation followed by morphometric and meristic data. Dim-3 showed that meristic data contributed 70% of the variation followed by morphometric and categorical data (Fig. 5B).

Figure 5. 

A MFA scatter plot showing the morphospatial relationships among the Dixonius species based on normalized morphometric, meristic, and color pattern characters B bar graphs showing the percent contribution of each data type to the overall variation in the data dimensions 1–4. The dashed red line in the bar graphs indicates the expected average value if the contributions of each data type were equal.

The ANOVAs and subsequent TukeyHDS tests demonstrated that Dixonius gialaiensis sp. nov. bears statistically different mean values between it and all other species in various combinations of characters (Tables 6, 7) and differs significantly from its sister species D. minhlei in head length (HL mean = 1.07 vs. 0.86, p = 0.000, respectively), in head width (HW mean = 0.89 vs. 0.82, p = 0.005, respectively), and in axilla to groin length, (AGr mean = 1.23 vs 1.32 p = 0.022) (Fig. 5; Tables 6, 7). Dixonius muangfuangensis sp. nov. also differed significantly from all other species in various combinations of characters and from its sister species D. lao it differs in head length (HL mean = 1.15 vs. 1.11, p = 0.004, respectively) and numbers of infralabials (IFL mean = 6.50 vs. 7.83, p = 0.026, respectively), and in numbers of supralabial at midorbital position (MO mean = 6.00 vs. 7.50, p = 0.00001, respectively) (Fig. 5, Tables 6, 7). Variation in all metric characters are visualized in Figs 6, 7.

Figure 6. 

Boxplot comparisons of meristic characters among the Dixonius species where interspecific statistical differences were recovered (see Table 7). Pale blue circles are means and the black horizontal bars are medians.

Figure 7. 

Violin plots of the normalized morphometric characters overlain with box plots showing the range, frequency, mean (white dot), and 50% quartile (black rectangle) of characters where interspecific statistical differences were recovered (see Table 7). New species in bold italics.