Field surveys were conducted in northern Vietnam by: 1) A. V. Pham and M. A. Vang in Sa De Phin Commune, Sin Ho District, Lai Chau Province in May 2015, and in Xuan Nha Nature Reserve, Van Ho District, Son La Province on 15 June 2016; 2) H. N. Ngo et al. in Phu Canh Nature Reserve, Da Bac District, Hoa Binh Province on 11 June 2016; 3) T. D. Le et al. in Xuan Son National Park, Du Village, Xuan Son Commune, Tan Son District, Phu Tho Province on 7 July 2016; and 4) T. S. Nguyen in Xuan Lien Nature Reserve, Bat Mot Commune, Thuong Xuan District, Thanh Hoa Province in July 2015 (Fig. 2).

Figure 2. 

Distribution map of the new populations of Tylototriton from North Vietnam, based on the following symbols: square (taxon 3, this study) the population from Sin Ho District, Lai Chau Province; diamond (taxon 1, this study) the upper one identifies the population from Van Ho District, Son La Province, and the lower one identifies the population from Thuong Xuan District, Thanh Hoa Province; triangle (taxon 2, this study) the upper one identifies the population from Tan Son District, Phu Tho Province, the middle one identifies the population from Da Bac District, Hoa Binh Province, and the lower one identifies the population from Lac Son District, Hoa Binh Province. The two populations identified by the circles represent T. asperrimus sensu stricto from China. High resolution remote sensing land cover information was extracted from “GLAD-UMD and SERVIR-Mekong, Natural annual tree canopy structure and surface water dynamics products, 2017” (lower left panel). Bioclimatic variables (right side) were extracted from remote sensing data provided by Deblauwe et al. (2016).

Specimens were anaesthetized and euthanized in a closed vessel with a piece of cotton wool containing ethyl acetate (Simmons 2002), fixed in 80% ethanol for five hours, and subsequently transferred to 70% ethanol for permanent storage. Tissue samples were preserved separately in 70% ethanol prior to fixation. Specimens were subsequently deposited in the collections of the Institute of Ecology and Biological Research (IEBR), Hanoi, the Tay Bac University (TBU), Son La Province, Vietnam, and the Zoologisches Forschungsmuseum Alexander Koenig (ZFMK), Bonn Germany.

Tissue samples from muscle of preserved specimens were extracted using the DNeasy blood and tissue kit, Qiagen (California, USA). A fragment of a mitochondrial gene, the NADH dehydrogenase subunit 2 (ND2), was amplified by PCR mastermix (Fermentas, Burlington, ON, Canada) using the primer pair, Sal_Nd2_F1 (5’- AAGCTTTTGGGCCCATACC-3’) (Nishikawa et al. 2013b) and a newly design primer TyloR1 (5’- GGTCTTTGGTCTYATTATCCTAA -3’). 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 minutes to activate the taq; with 40 cycles at 95 °C for 30 s, 58 °C for 45 s, 72 °C for 60 s; and the final extension at 72 °C for 6 minutes. PCR products were subjected to electrophoresis through a 1 % agarose gel (UltraPure™, Invitrogen, La Jolla, CA). Gels were stained for 10 min in 1 X TBE buffer with 2 pg/ml ethidium-bromide and visualized under UV light. Successful amplifications were purified to eliminate PCR components using a GeneJET™ PCR Purification kit (Fermentas). Purified PCR products were sent to FirstBase Malaysia for sequencing. We included 12 new samples from five populations distributed in north and north central Vietnam to another 21 known species’ samples of Tylototriton (Table 1). Additionally, five species were selected as outgroups: Echinotriton andersoni, E. chinhaiensis, Lyciasalamandra atifi, Notophthalmus viridescens, and Pleurodeles waltl, to root the tree (Qian et al. 2017; Wang et al. 2018).

The sequences were aligned in Clustal X v2 (Thompson et al. 1997) with default settings. Data were analyzed using maximum parsimony (MP) and maximum likelihood (ML) as implemented in PAUP 4.0b10 (Swofford 2001), and Bayesian analysis in MrBayes 3.2 (Ronquist et al. 2012). For 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) (Felsenstein 1985) was calculated using 1,000 pseudo-replicates and 100 random taxon addition replicates. All characters were equally weighted and unordered. For ML analysis, we used the optimal evolution model as selected by ModelTest v3.7 (Posada and Crandall 1998). To estimate BP in the ML analysis, a simple taxon addition option and 100 pseudo-replicates were employed. We considered BP values of ≥ 70 % to represent strong support (Hillis and Bull 1993).

For Bayesian analyses, we used the optimal model, GTR+I+G as selected by Modeltest v3.7, for ML and combined Bayesian analyses. Two simultaneous analyses with four Markov chains (one cold and three heated) were run for 10 million generations with a random starting tree and sampled every 1,000 generations. Log-likelihood scores of sample points were plotted against generation time to determine stationarity of Markov chains. The cutoff point for the burn-in function was set to 21, equivalent to 21,000 generations, in the Bayesian analysis, as -lnL scores reached stationarity after 21,000 generations in both runs. Nodal support was evaluated using Bootstrap in PAUP and posterior probability in MrBayes v3.2. Uncorrected pairwise divergences were calculated in PAUP*4.0b10.

We selected the relaxed-clock method (Drummond et al. 2006) to estimate divergence times. The obtained dataset was used as input for the computer program BEAST v1.8.0 (Drummond and Rambaut 2007). A priori criteria for the analysis were set in the program BEAUti v1.8.0. One calibration point, the split between the clade containing Tylototriton vietnamensis + T. panhai and the clade consisting of T. asperrimus and other related species, estimated for 12.4 ± 2.3 million years ago (MYA) (Wang et al. 2018), was used to calibrate the phylogeny. A general time-reversible (GTR) model using gamma + invariant sites with four gamma categories was employed along with the assumption of a relaxed molecular clock. As for the priors, we used all default settings, except for the Tree Prior category that was set to Yule Process, as recommended for species-level analyses. The codon-partitioned dataset was used for a single run. In addition, a random tree was employed as a starting tree. The length chain was set to 10⁷, and the Markov chain was sampled every 1,000 generations. After the dataset with the above settings was analyzed in BEAST, the resulting likelihood profile was then examined by the program Tracer v1.6 to determine the burn-in cutoff point. The final tree with calibration estimates was computed using the program TreeAnnotator v1.8.0 as recommended in the BEAST program manual.

All specimens were sexed by evaluating the size of the opening of the cloacal fissure: females show a puncture-like opening and males a wider slit-like opening. The holotype of T. asperrimus (ZMB 34089), collected from Guangxi Province, China, was loaned from the Zoologisches Museum Berlin (Museum für Naturkunde Berlin) and evaluated as a female (Fig. 3). In addition we investigated two other Vietnamese female specimens, one from IEBR: JJLR01195 from Pu Hoat Nature Reserve, Nghe An Province (T. notialis) and another from the Vietnam Forestry University (VFU) in Hanoi: VFUA.2009.8 (also known as voucher Tao1214 in Nishikawa et al. [2013b]) from Thuong Tien Nature Reserve, Hoa Binh Province (T. cf. asperrimus). Morphological comparisons were only performed among animals of the same sex, and only males had a sufficiently large number of specimens (N) to perform statistical analysis.

Figure 3. 

Holotype of Tylototriton asperrimus (ZMB 34089). In sequence: dorsal view; ventral view; lateral view with detail of ovaries; and detail of dorsal view of the head. Photographs T. Ziegler.

A total of 23 morphological characters were measured following Bernardes et al. (2017) to the nearest 0.01 mm with a digital caliper as follows: snout-vent length (SVL); head length (HL); head width (HW) measured behind the eyes and before the beginning of the parotoids; maximum head width (MHW); parotoid width (PW); maximum parotoid height (PH); eye length (EL); inter-eye distance (IE); inter-narial distance (IN); eye-narial distance (EN); lower jaw length (LJL) from tip of lower jaw to jaw angle; maximum upper eyelid length (UEL); humerus length (HUM); radius length (RAD); femur length (FEM); tibia length (TIB); axilla to groin (AG); trunk length (TkL) from wrinkle of throat to anterior tip of vent; length of the 5^th anterior dorsal nodule (L5N); width of vertebral cord (WVr) measured at the height of the 5^th nodule; cloaca length (ClL) length of cloaca muscle; tail length (TL); tail height (TH). The following ratios were calculated based on the measures above: total forelimb length (FORE); total hindlimb length (HIND); hind-limb to forelimb lengths (HIND/FORE); the relative length of radius to humerus (RAD/HUM); tibia to femur (TIB/FEM); and tail length to tail height (TL/TH).

The morphological comparison between the new taxa and their congeners were based on the specimen examination and the following literature: Fei et al. (1984), Böhme et al. (2005), Stuart et al. (2010), Shen et al. (2012), Nishikawa et al. (2013a), Nishikawa et al. (2013b), and Yang et al. (2014). When measurements were involved, only the ones taken in similar ways were found suitable for comparison and used as reference.

We first compared the morphological characters of males between the two clades originating on both sides of the Da River: the western clade from Son La and Thanh Hoa provinces (referred to as taxon 1) and the eastern clade from Hoa Binh and Phu Tho provinces (referred to as taxon 2; for reference see Fig. 2). Subsequently we compared the above-mentioned males (jointly referred to as T. cf. asperrimus) and the males originating from Lai Chau Province (referred to as taxon 3).

The statistical analyses had to be conducted on different subsets of morphological characters according to data availability. Morphological characters that could not be obtained for all the species had to be excluded from the overall analysis. These included: PW, PH, EL, IE, UEL, AG, and ClL. Whether the measured morphological characters showed a linear increase with body size was analyzed through correlation analyses (see Suppl. material 1). Accordingly, measurements of morphological characters and character ratios were standardized by SVL (R[character]: % SVL) to exclude the effect of body size, and log-transformed. A Principal Component Analysis (PCA) was tested by a one-way Analysis of Variance (ANOVA) between populations. Because morphological traits within individuals are not independent of each other, comparisons between different morphological traits of species were based on Multivariate Analysis of Variance (MANOVA) and proceeding to ANOVA and Tukey HSD test only if the MANOVA yielded a significant result (i.e., ‘protected ANOVA’ (van Ende 2001). Roy’s Greatest Root was chosen as test of significant differences among groups in the MANOVA procedure.

Significance levels were set to 95 %. All statistical analyses were performed in R v 3.1.2, the vegan package was used to calculate PCA (Oksanen et al. 2015).

Climatic information at the sample sites were extracted from remote sensing data (Deblauwe et al. 2016). Representing averages across last decades with a spatial resolution of 0.1°, the following bioclimatic variables were available: Annual Mean Temperature BIO1, Mean Diurnal Range BIO2, Isothermality BIO3, Temperature Seasonality BIO4, Max Temperature of Warmest Month BIO5, Min Temperature of Coldest Month BIO6, Temperature Annual Range BIO7, Mean Temperature of Wettest Quarter BIO8, Mean Temperature of Driest Quarter BIO9, Mean Temperature of Warmest Quarter BIO10, Mean Temperature of Coldest Quarter BIO11, Annual Precipitation BIO12, Precipitation of Wettest Month BIO13, Precipitation of Driest Month BIO14, Precipitation Seasonality BIO15, Precipitation of Wettest Quarter BIO16, Precipitation of Driest Quarter BIO17, Precipitation of Warmest Quarter BIO18, and Precipitation of Coldest Quarter BIO19.

The combined matrix contained 1036 aligned characters. Of those, 370 were parsimony informative. MP analysis of the dataset recovered 2 most parsimonious trees with 1400 steps (CI = 0.54; RI = 0.65). Our phylogenetic analyses recovered the Vietnamese T. cf. asperrimus as a sister taxon to T. asperrimus from China with strong support values from all analyses (MP_BP = 90, ML_BP = 88, PP = 100) (Fig. 4). The genetic differences between Vietnamese populations and the Chinese lineage were 3.3 to 3.6 % for the population from Son La Province; 3.2 to 3.4 % for the population from Thanh Hoa Province; 3.3 to 3.6 % for the population from Phu Tho Province; 3.2 to 3.6 % for the population from Da Bac District, Hoa Binh Province; and 3.4 to 3.5 % for the population from Lac Son District, Hoa Binh Province, respectively (Table 2).

Figure 4. 

Phylogram based on the Bayesian analysis. Number above and below branches are MP/ML bootstrap values and Bayesian posterior probabilities (> 50 %), respectively. Dashes represent values < 50 %. Sample AB769531 is from Nishikawa et al. 2013b.

Furthermore, our genetic analyses identified different lineages within the Vietnamese clade of T. cf. asperrimus. The genetic variation between taxon 1 and taxon 2 varied between 2.5 % (between Thanh Hoa and Hoa Binh populations) and 3.1 % (between Son La and Phu Tho populations). In contrast, within-population differences were only 0.0 to 0.6 % in taxon 1 and 0.1 to 0.9 % in taxon 2.

The population from Lai Chau Province turned out to be a distinct and basal lineage within a weakly supported clade, including T. notialis, T. asperrimus from China, and taxon 1 and taxon 2 from Vietnam (Fig. 4). In this case the genetic differences of taxon 3 to the topotypical population of T. asperrimus ranged between 4.1 to 4.2 % to taxon 1 between 3.6 to 4.0 %, and to taxon 2 between 4.1 to 4.5 % (see Table 2 for genetic distances). Our time estimates are very similar to those generated by Wang et al. (2018), and the results show that T. asperrimus from China split from taxon 1 about 2.5 MYA (95% highest posterior densities – 95% HPD = 1.4–3.7), while taxon 3 diverged from the two taxa approximately 3.4 MYA (95% HPD = 2.3–4.8) (see Suppl. material 2).

Vietnamese species compared to the Chinese holotype

This comparison was only based on three female specimens: the holotype of T. asperrimus, one from Hoa Binh Province (taxon 2), and one from Nghe An Province (T. notialis) (Table 3). Due to the lack of replicates it was not possible to perform statistical analyses between the Chinese and the Vietnamese clades. After correcting the absolute measures to ratios of snout-vent length, the most prominent differences between the female of taxon 2 and the female holotype of T. asperrimus from China were: a wider and longer head (MHW = 28.99, HL = 29.15 in taxon 2 vs. MHW = 25.35, HL = 26.60 in T. asperrimus), a longer lower jaw (LJL = 17.40 in taxon 2 vs. 14.85 in T. asperrimus), and higher values for most of the measured head features (including the distance between the eyes) for taxon 2. The exceptions were found in the distance between eye and nostril (EN = 6.37 in T. asperrimus vs. 4.86 in taxon 2) and head width (HW = 18.54 in T. asperrimus vs. 13.67 in taxon 2) which in these cases the values were higher in T. asperrimus. The female from taxon 2 also had higher values for tail length (TL = 85.62 in taxon 2 vs. 77.28 in T. asperrimus), cloacal muscles (ClL = 11.39, ClW = 7.63 in taxon 2 vs. ClL = 7.01, ClW = 3.65 in T. asperrimus), and vertebra width (WVr = 4.46 in taxon 2 vs. 2.80 in T. asperrimus). The female from China had a longer trunk length (TkL = 74.58 in T. asperrimus vs. 69.67 in taxon 2). The female from Nghe An Province differed by having the smallest eye length, the shortest distance between both eyes, the smallest glandular warts and by having the longest limbs, while other measurements did not separate it from other lineages.

Comparisons within T. cf. asperrimus from Vietnam

The comparison between taxon 1 and taxon 2 included only males. Absolute measures and ratios of species’ morphological traits corrected by snout vent length are shown in Table 4. Taxon 1 and taxon 2 did not differ in their respective SVL (t-test = -1.55, DF = 18, p = 0.14). Taxon 1 presented wider head than taxon 2 (MHW = 27.37 ± 1.67, HW = 19.86 ± 0.95 for taxon 1 vs. MHW = 25.11 ± 0.81, HW = 18.75 ± 0.99 taxon 2). The ratios of EN and IN also differed between lineages, with taxon 1 having a relatively longer snout (EN = 6.16 ± 0.68) than taxon 2 (5.50 ± 0.37) and taxon 2 having a wider snout (IN = 8.58 ± 0.57) than taxon 1 (7.94 ± 0.85). Taxon 2 showed the highest variation range in limb data. The ratio of FEM was longer in taxon 2 (12.59 ± 0.78) than in taxon 1 (10.99 ± 0.67), as well as the ratio of TIB (23.69 ± 2.35 in taxon 2 vs. 22.41 ± 0.61 in taxon 1), which together also resulted in longer hind limbs (HIND = 36.28 ± 2.89 in taxon 2 vs. 33.40 ± 0.99 in taxon 1). The ratio of the fore-limbs on the other hand was alike between lineages. Although taxon 2 presented lower minimum values for both RAD (18.80) and HUM (9.30) than in taxon 1 (RAD = 20.08 and HUM = 11.04). The tail in taxon 1 was longer (87.06 ± 5.25) and less high (11.79 ± 1.14) than in taxon 2 (TL = 84.61 ± 5.20, TH = 12.43 ± 2.46). Taxon 1 showed longer trunk (70.28 ± 3.41) than taxon 2 (68.78 ± 1.24).

The statistical analysis was based on nine males of taxon 1 and ten males of taxon 2. A PCA analysis resulted in six principal components (PC) explaining 87 % of the total variation. The first two PCs accounted for 52 % of the variation. The scatterplot between PC1 and PC2 showed a clear separation of the two clades, with only a small overlap area (Fig. 5A).

The head related data (MANOVA: F_{1, 17} = 11.75, DF = 6, p < 0.001), and the limb related data (MANOVA: F_{1, 17} = 5.10, DF = 9, p = 0.01) were significantly different between the two lineages. Tail and dorsal morphological traits were not significantly different (MANOVA: F_{1, 17} = 1.42, DF = 6, p = 0.3). Our results identified MHW, HW, EN, IN, RAD/HUM, FEM, TIB/FEM, HIND, and HIND/FORE as important traits separating both lineages (Table 5). Taxon 1 has a wider head (both as MHW +8.3 %; and as HW +5.6 %) and a longer snout (EN +10.7 %). Taxon 2 has a wider snout (IN +7.5 %) (Fig. 6).

Regarding the limb data, FEM was 12.7 % longer on taxon 2, as well as the overall hind-limb length (HIND +7.9 %) and the ratio of HIND to FORE (+5.6 %). On the contrary, the ratios of tibia to femur (TIB/FEM +10.2 %) and radius to humerus (RAD/HUM +26.2 %) were larger in taxon 1 (Fig. 7).

Figure 5. 

Scatterplot between PC1 and PC2 of the morphological characters corrected to SVL and log–transformed, for A taxon 1 and taxon 2 of the Vietnamese Tylototriton cf. asperrimus B the head- and dorso- related data of taxon 3 from Lai Chau Province and T. cf. asperrimus from Vietnam sensu lato; and C the limb related data of taxon 3 from Lai Chau Province and T. cf. asperrimus from Vietnam sensu lato. In the graphics T. cf. refers to T. cf. asperrimus.

Figure 6. 

Boxplot of the most differing characters related to head and dorsal values between taxon 1, taxon 2, and taxon 3. Characters were corrected to SVL and log–transformed. For abbreviations see Materials and methods.

Figure 7. 

Boxplot of the most differing characters related to limb values between taxon 1, taxon 2, and taxon 3. Characters were corrected to SVL and log–transformed. For abbreviations see Materials and methods.

Comparison of taxon 3 from Lai Chau Province with taxon 1 and taxon 2

This analysis is only based on males. Absolute measures and ratios of species’ morphological traits corrected by snout vent length are shown in Table 4. All three taxa had similar measures for SVL, TL, and TH. The narrowest head was recorded in taxon 3 (HW 15.29 ± 1.04; 19.86 ± 0.95 in taxon 1; 18.75 ± 0.99 in taxon 2) and its maximum values were still below the minima recorded for taxon 1 and taxon 2 (max HW = 16.40 in taxon 3; min HW = 18.44 in taxon 1 and = 16.98 in taxon 2). The snout length was longer in taxon 3 (EN 6.77 ± 0.78), than in taxon 1 (6.16 ± 0.68) or taxon 2 (5.50 ± 0.37). HUM was longer in taxon 3 (13.89 ± 1.69) and showed a maximum range (11.99 to 16.57) not repeated in taxon 1 (11.04 to 13.89) nor in taxon 2 (9.30 to 13.54). Consequently, the sizes of the fore-limbs were also longer in taxon 3 (35.67 ± 2.13, range = 33.88 - 38.9) than in taxon 1 (33.67 ± 0.89, range = 32.43–35.35) and in taxon 2 (33.68 ± 2.11, range = 30.36 - 36.97). In taxon 3 the width of the vertebral cord (WVr 3.53 ± 0.25) and the length of the rib nodules (L5N 3.31 ± 0.87) were wider than in taxon 1 (3.06 ± 0.37, and 2.67 ±0.39, respectively) and taxon 2 (2.90 ± 0.37, and 2.91 ±0.41, respectively). Trunk length, on the other hand was shorter in taxon 3 (TkL 66.22 ± 3.86), than in taxon 1 (70.28 ± 3.41) or taxon 2 (68.78 ± 1.24) (Table 4).

The data set of head and dorsal morphological traits was based on 19 observations of taxon 1 and taxon 2 together and three observations of taxon 3 from Lai Chau. A PCA identified five principal components (PCs) which together explained 84 % of the morphological variation (cumulative explanation of the first 3 PCs = 66 %; of the first 4 PCAs = 75 %). The first two PCs accounted for 48 % of the variation graphically showing a clear separation of the two clades (Fig. 5B). HW, EN, WVr, and L5N were identified as the characters differentiating between the species (MANOVA: F_{1, 20} = 20.52, p < 0.001) (Table 5). Head width (HW) was 21 % smaller in taxon 3 than in taxon 1 and taxon 2 (F_{2, 19} = 36.79, p < 0.001), and the size of the rib nodules (L5N) was 15 % longer in taxon 3 than in taxon 1 and taxon 2 (F_{2, 19} = 6.59, p < 0.01). The two remaining characters were only different between taxon 3 and taxon 2. Both the snout length (EN) and the width of the vertebral cord (WVr) were longer in taxon 3 than in taxon 2, by 17 % (F_{2, 19} = 7.21, p < 0.01) and 16 % (F_{2, 19} = 3.45, p < 0.05), respectively (Fig. 6).

The limb data included 21 observations of taxon 1 and taxon 2 together and four of taxon 3 and resulted in a PCA with three PCs explaining 88 % of the variation. The overall MANOVA (F_{1, 23} = 1.92, p = 0.13) was not significantly different between both lineages (Fig. 5C).

Macroclimatic comparison

Our data show that T. asperrimus in Guangxi, China experiences the lowest temperatures during the coldest months (3–6 °C) than any of the remaining three taxa in North Vietnam (12 °C). This species also shows the highest amount of precipitation during the coldest (169–233 mm vs. 38–80 mm for the remaining three taxa) and driest (170–180 mm vs. 38–80 mm for the three remaining taxa) quarter of the year, as well as in the driest month (26–44 mm vs. 4–10 mm for the three remaining taxa) (Table 6).

Genetic and morphological differences found in this study support the taxonomic separation between T. cf. asperrimus from Vietnam and T. asperrimus sensu stricto (from China), thus confirming the distinctness of the Vietnamese clade. Furthermore, we uncovered genetic and morphological variations within the Vietnamese T. cf. asperrimus clade. However, based on our current knowledge these should be evaluated with caution regarding taxon 1 and taxon 2. Therefore these taxa are treated herein at the subspecies level until further evidence is presented. In addition, due to distinct morphological and molecular divergence, the population from Lai Chau Province was revealed to be distinct at the species level.