Results
I. Phylogenetics and population structure of B. gargarizans
A total of 77 mtDNA (ND2) haplotypes were detected for B. gargarizans , with a mean intraspecific pairwise distance of 0.0535 (0.0018–0.0722) (Table 3 ). There were 153 polymorphic sites, of which 63 were singleton-variable and 90 were parsimony-informative sites. The BI analysis recovered eight clades (Fig. 2 ). For intra-clade genetic distances, the highest values were in clade G, while the lowest values were in clade A (Table 1 ). Compared to results of Zhan and Fu (2011), a new clade (Clade H) was identified, containing samples exclusive to Northeast Asia (Korean Peninsula, Northeast China [Heilongjiang Province], and Russia). Clade H contained 18 haplotypes (Table 2 ). Pairwise distances between the Clade H and the other seven clades ranged between 0.0451–0.0722 (Table 1 ).
Haplotypes of individuals from the four geographic regions identified by Zhan and Fu (2011) were found in multiple genetic clades: W region had haplotypes in four genetic clades (A, C, F, G), with clade F being exclusive to this region; C region had haplotypes in five genetic clades (A, B, C, E, G), with clade E being exclusive to this region; the SE region had haplotypes in four genetic clades (B, C, D, G); and the NE region had haplotypes in four genetic clades (A, C, D, H), with clade H being exclusive to this region (Table 2 ).
Genetic diversity indices are summarized in Table 3 . The number of singleton variable site (S) was 32 and nucleotide diversity (π) was 0.04412. Areas with high genetic diversity indices were the W region (68 polymorphic sites, 49 parsimony-generic sites, 19 single-variable sites) and the C region (nucleotide diversity = 0.03793). Pairwise distances and nucleotide diversity were lowest in the NE region (Table 3 ).
The multilocus NJ phylogram from POFAD recovered seven major groups, with one being unique to our study (Fig. 3 ). For Zhan and Fu (2011)’s six groups, there was no correlation with either region or altitude, with individuals from the W region included in all groups. Members of the newly identified group in this study (Fig. 3 ) came from the SE region (Fujian and Zhejiang Provinces in China) and the NE region (Heilongjiang Province in China, Korea, and Russia).
II. B. gargarizansin Northeast Asia
The ND2 -based genetic network was divided into two major groups (Fig. 4 ). Cluster A included all haplotypes of Clade H (18 haplotypes), while Cluster B included three clades: Clade A (hap30, hap43, hap44), Clade C (hap15), and Clade D (hap31). There was clear genetic differentiation between Cluster A and B, despite their close geographic distance. In particular, B. gargarizans from China’s Heilongjiang Province and Russia were included in Cluster A, although they were closer geographically to some individuals in Cluster B (Jilin and Liaoning Provinces). The mean pairwise difference between the two groups was 0.0503 (0.0451–0.0533) (Table 1 ).
III. Divergence dating analysis
In general, we recovered broad confidence intervals of the divergence time estimates (Fig. 5 , Table S6 ). Bufo gargarizansis estimated to have diverged 7.29 Ma (2.83–12.69 Ma), while the major clades within the species are estimated with mean ages of 2.25–4.30 Ma. Our analysis inferred that the divergence pattern of B. gargarizans is from West (mean ages 4.76 Ma, 4.30 Ma) to Southeast (3.46 Ma) in the Pliocene, followed by Northeast (2.25 Ma) in the early Pleistocene. Divergence within Clade H occurred during the Holocene period (1.07 Ma, 0.61 Ma).