Do Cryptic Species Exist in Hoplobatrachus rugulosus? An Examination Using Four Nuclear Genes, the Cyt b Gene and the Complete MT Genome

he Chinese tiger frog Hoplobatrachus rugulosus is widely distributed in southern China, Malaysia, Myanmar, Thailand, and Vietnam. It is listed in Appendix II of CITES as the only Class II nationally-protected frog in China. The bred tiger frog known as the Thailand tiger frog, is also identified as H. rugulosus. Our analysis of the Cyt b gene showed high genetic divergence (13.8%) between wild and bred samples of tiger frog. Unexpected genetic divergence of the complete mt genome (14.0%) was also observed between wild and bred samples of tiger frog. Yet, the nuclear genes (NCX1, Rag1, Rhod, Tyr) showed little divergence between them. Despite this and their very similar morphology, the features of the mitochondrial genome including genetic divergence of other genes, different three-dimensional structures of ND5 proteins, and gene rearrangements indicate that H. rugulosus may be a cryptic species complex. Using Bayesian inference, maximum likelihood, and maximum parsimony analyses, Hoplobatrachus was resolved as a sister clade to Euphlyctis, and H. rugulosus (BT) as a sister clade to H. rugulosus (WT). We suggest that we should prevent Thailand tiger frogs (bred type) from escaping into wild environments lest they produce hybrids with Chinese tiger frogs (wild type).


Introduction
Vertebrate mitochondrial DNA (mtDNA) is a closed circular genome which is approximately 16-20 kb long [1]. This genome typically contains 37 genes (2 rRNAs, 22 tRNAs, and 13 protein-coding genes) and a long non-coding region called the control region or D-loop region [1,2]. The mitochondrial (mt) genome has several valuable characteristics, including small size, fast evolutionary rate, relatively conserved gene content and organization, maternal inheritance, as well as limited recombination [3][4][5][6]. Mitochondrial genomes are useful molecular markers in cryptic species identification because of the differences in compositional features, divergence of protein-coding genes, number and size of non-coding regions, and gene arrangement [7][8][9].
According to the Amphibian Species of the World 5.6, an online reference (19 Mar. 2014) [10], about 6344 species of Order Anura exist worldwide. Among them, 319 species and subspecies can be found in China [11]. The genus Hoplobatrachus with five species is important among the Dicroglossini of Dicroglossinae [12][13][14]. Only one of these species, known as the Chinese tiger frog occurs naturally in wild environments in China [10,11,15]. This frog has been identified as Hoplobatrachus rugulosus [13], although some Chinese researchers insist that it should be named Hoplobatrachus chinensis, believing this to be a senior synonym of H. rugulosus [15]. The Chinese tiger frog is listed in Appendix II of CITES as the only Class II nationally-protected frog in China. A bred tiger frog introduced to China from Thailand is called the Thailand tiger frog by Chinese, and has also been identified as H. rugulosus according to Frost's taxonomic methods [10]. Thailand tiger frogs are bred in many farmers for local meat consumption. But more and more Thailand tiger frogs have been captured in the field after escaping from farms, which may affect the diversity of local Chinese tiger frogs. Alam et al. [16] found high divergences in Hoplobatrachus using Cyt b, 12S rRNA, and 16S rRNA genes. They suggested that H. chinensis (= H. rugulosus) may be subdivided into more than one species. Pansook et al. [17] found two distinct clades in H. rugulosus from Thailand using Cyt b gene, and also suggested that two distinct species may be present in H. rugulosus. In this study, we also found a high divergence between Chinese tiger frog and Thailand tiger frog using Cyt b gene. A cryptic species complex is a group of organisms that are typically very closely related yet their precise classification and relationships cannot be easily determined, although some can be separated by DNA sequence analyses [18][19][20]. Using the complete mt genome, Alam et al. [21] and Yu et al. [22] found duplications of the ND5 gene and the control region in Hoplobatrachus tigerinus and H. rugulosus. The current work aimed to investigate further the differences between Chinese tiger frog and Thailand tiger frog, as well as to determine whether cryptic species are present in H. rugulosus. Accordingly, we determined the complete mt genome sequences of Thailand tiger frog (bred type (BT)), and then compared the differences in gene arrangement, base compositional features, and genetic divergences of mt genes between Chinese tiger frog (wild type (WT)) and Thailand tiger frog (BT). The protein structures of the ND5 genes in Hoplobatrachus were also compared. Additionally, we sequenced nuclear genes (NCX1, Rag1, Rhod, and Tyr) of wild and bred tiger frogs to examine their genetic divergence. We also performed molecular phylogenetic analyses to discuss the relationship between Chinese tiger frog and Thailand tiger frog, and all available dicroglossids including Occidozyga martensii, based on the 11 mt protein-coding genes using five Ranidae species as out-groups. We follow the names proposed by Frost et al. [12] to avoid taxonomic confusion.

Ethical statement
Although the Chinese tiger frog is a protected species, the Chinese tiger frog samples (wild tiger frogs) we used were donated by the Committee of Forest Administrative Bureau of Jinhua (CFABJ), People's Republic of China in 2000 and 2010. The officers of Forest Administrative Bureau of Jinhua have seized these wild frogs, which died during illegal captivity. The Thailand tiger frog samples (bred tiger frogs) were purchased from various farms from Jinhua, Zhejiang province and sacrificed using ether in our laboratory. The husbandry and breeding procedures of the Thailand tiger frogs in the farms were carried out under the Animal Husbandry Law of the People's Republic of China. The study protocol was reviewed and approved by the Committee of Animal Research Ethics of Zhejiang Normal University.
Sample and DNA extraction. All samples were stored at-70°C in the Institute of Ecology, Zhejiang Normal University. Information for all samples is shown in S1 Table. Whole genomic DNA was extracted from frozen tissue sample of thigh muscle of tiger frogs using a standard proteinase K/SDS digest extraction method followed by phenol-chloroform isolation and ethanol precipitation [23]. A sample of the Thailand tiger frog (No. THW1) was used to amplify the complete mt genome; other samples were used to amplify the partial Cyt b gene.
All PCRs were performed using a MyCycler Thermal Cycler (Bio-Rad, Hercules, CA, USA). TaKaRa Ex-Taq and LA-Taq Kits (Takara Biomedical, Dalian, China) were used for normal PCR and LA-PCR reactions respectively. The normal PCR was carried out in a 50 μl reaction mixture containing 5 μl of buffer (10× concentration), 4 μl of MgCl 2 (25 mM), 4 μl of dNTP (2.5 mM), 2 μl of each primer (20 μM), 0.25 μl of Ex-Taq polymerase, 30.75 μl of sterile distilled water, and 2 μl of template DNA. The PCR reactions consisted of an initial denaturation at 95°C for 4 min; 35 cycles of denaturation at 94°C for 40 s plus annealing at 48°C to 60°C for

Sequence assembly and analysis
Sequences were checked and assembled using SeqMan (Lasergene version 5.0) [32]. The locations of the 13 protein coding genes and 2 rRNA genes were determined by comparing the homologous sequences of other anurans using Clustal W in Mega 5.0 [25,33]. All tRNA genes, except tRNA Ser (AGY) and tRNA Cys genes, were identified by their cloverleaf secondary structure in tRNA-scan SE 1.21 [34] using the default parameters, and tRNA Ser (AGY) and tRNA Cys genes were determined by comparing the homologous sequences of other anurans. The organization of the H. rugulosus (BT) mt genome was formed using GenomeVx (Fig 1)

Molecular phylogenetic analyses
A total of 55 sequences of the Cyt b gene, including 19 from this study, in-group species from Pansook et al. [17] and out-group species (S1 Table), were used to evaluate the divergence between the Chinese tiger frog and the Thailand tiger frog. All 55 sequences yielded 565 bp of Cyt b gene fragment, including 231 variable and 197 parsimony-informative sites. Phylogenetic relationships were constructed by neighbor joining (NJ) analysis in Mega 5.0 [33] using Kimura 2-parameter model, gamma distributed, gamma parameter 6, with both transitions and transversions included. Statistical support was evaluated using 1000 bootstrap replicates. Using the data with 20 sequences yielded 1624 bp of Rag1, Rhod, and Tyr genes (excluding NCX1 gene as no squeence of Hoplobatrachus is reported in GenBank), including 36 variable sites. The outgroup was combined H. occipitalis (Rag 1 gene HM163613, Tyr gene AJ564729) and H. tigerinus (Rhod gene AB489039). Phylogenetic relationships were also constructed by neighbor joining (NJ) analysis in Mega 5.0 [33] with parameter set as above.
To confirm further the phylogenetic relationships of the Chinese tiger frog and the Thailand tiger frog among dicroglossids, 16 available complete mt genomes of anurans based on the addition of two Hoplobatrachus mt genomes, including five species of Ranidae ((Babina adenopleura, Pelophylax nigromaculata, Pelophylax plancyi, Odorrana ishikawae, and Odorrana tormotus) [22,36,37]) as out-groups, were retrieved from GenBank (S2 Table). We constructed the phylogenetic trees based on a concatenated set of 11 mt protein-coding genes excluding ND5 because of the high sequence divergence of its two copies in H. rugulosus BT and ATP8 because of its low number of informative sites in anurans.
The nucleotide sequences for each of the 11 mt protein-coding genes were aligned using Clustal W in Mega 5.0 [25,33]. To select the conserved regions of the sequences, each alignment was analyzed with Gblocks 0.91b [38] using default settings. We concatenated the alignments of the 11 mt protein-coding genes, and recovered an alignment consisting of 2791 amino acid residues. An alignment of 8373 nucleotides was obtained using the amino acid alignment as the backbone reference. A saturation analysis was performed for subsets with first, second, and third codon positions using DAMBE 4.2.13 [39]. The results showed that the third codon positions were saturated. Consequently, they were excluded from the final nucleotide alignment and an alignment of 5582 nucleotides was obtained. Maximum parsimony (MP) analysis with the nucleotide dataset was performed using PAUP Ã 4.0b10 [40]. A total of 1000 bootstrap replications were generated, each with 10 replications with random taxon order.
Model selection for the nucleotide dataset was performed with Modeltest version 3.7 [41]. The GTR+I+G model was chosen for the likelihood and Bayesian analyses. Maximum likelihood (ML) analysis of the nucleotide dataset was performed using PAUP Ã 4.0b10 [40] with 1000 bootstrap replications. Bayesian inference (BI) analysis of the nucleotide dataset was performed with MrBayes 3.1.2 [42]. Eight chains were run in parallel for 10 000 000 generations, with trees sampled every 1000 generations. The burn-in sizes for both nucleotide datasets were determined by checking convergences of −log likelihood (−lnL) values, and the first 50 000 generations were discarded. Bayesian posterior probabilities were calculated according to the remaining set of trees. All Markov chain Monte Carlo runs were repeated twice to confirm consistent estimation of the posterior parameter distributions.

Genome organization and arrangement
The complete mt genome of H. rugulosus (BT) was 20 926 bp long and contained 14 proteincoding genes (including the extra copy of the ND5 gene), 2 ribosomal RNAs, 23 transfer RNA genes (including the extra copy of the tRNA Met gene), and 2 non-coding regions (including the extra copy of the control region (D-loop)) ( Table 2). In H. rugulosus (BT), the distinctive features included a modified cluster of rearranged tRNA genes (TPF tRNA gene cluster), the tandem duplication of tRNA Met genes (Met1 and Met2), the translocation of tRNA Leu (CUN) and ND5 genes, and the two copies of D-loop-ND5 regions. The tRNA Leu (CUN) gene was located between the two D-loop-ND5 regions, rather than in the typical LTPF tRNA cluster (Fig 1 and Table 2). The first D-loop-ND5 region was located between the Cyt b and tRNA Leu (CUN) genes, and the second D-loop-ND5 region was located between the tRNA Leu (CUN) and tRNA Thr genes (Fig 1 and Table 2). Most of these genes were coded on the H-strand, except for ND6 and eight tRNA genes, as reported in other anurans (Fig 1).

Protein-coding genes
There were six reading frame overlaps in the mt genome of H. rugulosus (BT): COI and tRNA Ser (UCN) shared nine nucleotides; ATP8 and ATP6 shared seven nucleotides; ND4L and ND4 shared seven nucleotides; ND2 and tRNA Trp shared two nucleotides; COIII and tRNA Gly shared one nucleotide; and ND3 and tRNA Arg shared two nucleotides. Other overlaps are shown in Table 2. Two ND5 genes (with 84.1% similarity) were found in H. rugulosus (BT) as well as H. tigerinus (two identical ND5 genes) and H. rugulosus (WT) (two identical ND5 genes). The first ND5 gene (ND5-1) was located between the D-loop and tRNA Leu (CUN) , and the second ND5 gene (ND5-2) was located between the D-loop and tRNA Thr (Fig 1 and Table 2). The dataset comparing the two ND5 genes in H. rugulosus (BT) included 288 variable sites in a total of 1848 aligned nucleotide sites and 131 variable sites in amino acid sequence over a total of 616 alignment amino acid sites. The codon frequency of the two ND5 genes in H. rugulosus (BT) is shown in Table 3.
Protein-coding genes in H. rugulosus (BT) begin with ATG as the start codon, except COI with ATA and the ND6 gene, for which H. rugulosus (WT) used ACA-H. rugulosus (BT) used ATG. ND5-1, ATP8, COIII, ND3, ND4L, and Cyt b genes are terminated with TAA as the stop codon, COI and ND6 end with AGG, ND2 and ND5-2 end with TAG, and the other four protein-coding genes end with an incomplete stop codon (a single stop nucleotide T) ( Table 2). In  (Table 2). Three protein-coding genes surprisingly showed different lengths between H. rugulosus WT and BT: ND3 differed by 3 bp (one amino acid), ND5-1 differed by 6 bp (two amino acids), and ND5-2 differed by 12 bp (four amino acids).

Ribosomal and transfer RNA genes
In H. rugulosus (WT), the 12S rRNA gene (940 bp long) was located between tRNA Phe and tRNA Val genes, and the 16S rRNA gene (1587 bp long) was located between tRNA Val and tRNA Leu (UUR) genes. The mt genome of H. rugulosus (BT) contained 23 tRNA genes (including the extra copy of tRNA Met gene) (Fig 2), that were interspersed in the genome and ranged in size from 65 bp to 73 bp. As reported in anurans, all tRNA genes can be folded into typical cloverleaf secondary structures with the known exception of the tRNA Ser(AGY) gene (Fig 2).
The divergence of 12S rRNA and 16S rRNA genes between H. rugulosus (WT) and H. rugulosus (BT) was 7.6% and 6.6%, respectively, using the uncorrected p-distance model. Table 4 shows the divergences (uncorrected p-distances) of the 23 tRNA genes, 12S rRNA genes, and 16S rRNA genes among H. tigerinus and H. rugulosus (WT and BT).
The putative origin of the L-strand replication (O L ) [51] was located in the WANCY tRNA gene cluster between the tRNA Asn and tRNA Cys genes (Fig 1). This region was 30 bp long and had the potential to fold into a stem-loop secondary structure with a stem formed by 9 paired nucleotides and a loop of 12 nucleotides. In this study, we analyzed 12 sequences of O L in Dicroglossidae, and found that the sequence of H. rugulosus (BT) was similar to those of E. hexadactylus, H. tigerinus, H. rugulosus (WT), Nanorana pleskei, and Quasipaa spinosa (Fig 3).   The divergence between two clades of H. rugulosus was 11.5% using the uncorrected p-distance model. BI, ML, and MP phylogenetic analyses based on the nucleotide dataset of 11 protein-coding genes had a similar topology (Fig 5), which is consistent with Zhang et al. [52]. In this study, we recovered topological relationships among dicroglossid clades with high bootstrap and posterior probability (Fig 5). O. martensii (Occidozyginae: Occidozygini) occupied the basal phylogenetic position among the dicroglossid frogs (posterior probabilities 1.00 in BI, bootstrap value 100% in ML and MP). The nucleotide dataset also favored a topology that placed Hoplobatrachus as a sister clade to Euphlyctis (posterior probabilities 1.00 in BI, bootstrap value 100% in ML and MP), and H. rugulosus (BT) was a sister clade to H. rugulosus (WT) (posterior probabilities 0.93 in BI, bootstrap value 54% in ML and 83% in MP).

Phylogenetic analyses in Hoplobatrachus
Analyses in Hoplobatrachus using the nuclear genes.   gene including H. tigerinus (AB277358) and H. occipitalis (AJ564729) from GenBank with 14 variable sites of 532 nucleotides, were used to evaluate the divergence between the Chinese tiger frog and the Thailand tiger frog, respectively. S3-S6 Tables show the nucleotide genetic divergences among the samples with NCX1, Rag1, Rhod, and Tyr genes, respectively. The nucleotide mean genetic divergences between H. rugulosus (BT) and H. rugulosus (WT) using NCX1, Rag1, Rhod, and Tyr genes were 0.7%, 0.3%, 0.1%, 0.1%, respectively. We found fixed different nucleotide sites in the Rag1 gene: C and G in H. rugulosus (BT) but T and A in H. rugulosus (WT) on 386th and 695th site, respectively. Using NJ analysis based on the data of Rag1, Rhod, and Tyr genes, we also found two clades (clade 1 and 2) in Hoplobatrachus (S1 Fig).

Evidence of cryptic species
The divergences of partial Cyt b gene and mt genome in H. rugulosus. The taxonomic problem of H. rugulosus has been debated by a number of researchers [10,13,16,53]. Alam et al. [16] and Pansook et al. [17] found high sequence divergences in H. rugulosus. In our study, we found high divergence (13.8%) between H. rugulosus WT and BT using the Cyt b gene, and high divergence (11.5%) between two different clades (clade 1 and 2) of the phylogenetic relationship in H. rugulosus (Fig 4). Based on the Cyt b gene divergence we support the suggestion [17] that H. rugulosus is a cryptic species complex.
Nucleotide divergence within the complete mt genome was 14.0% between H. rugulosus WT and BT, which is very high but not beyond the genetic divergence found in cryptic species of diverse other taxa such as Ciona intestinalis, Friesea grisea, and Taenia taeniaeformis [7][8][9]. Deep nucleotide divergence in protein-coding genes (11.2-24%) was also found using the uncorrected p-distance model between H. rugulosus (WT) and H. rugulosus (BT), which is similar to the values for cryptic species in some frogs [18,20].
The two ND5 genes in H. rugulosus (WT) and H. tigerinus were identical within each species whereas the two ND5 genes in H. rugulosus (BT) had only 84.1% similarity sequence. Although the transmembrane domains and three-dimensional structure analyses of ND5 proteins had differences among H. tigerinus, H. rugulosus (BT), and H. rugulosus (WT), the two ND5 genes in H. rugulosus (WT) had an identical tertiary structure (Fig 7). The predicted results (Fig 7) showed that ND5-1 and ND5-2 of H. rugulosus (BT) had a similar structural plan, which suggested a common origin in the two ND5 genes. A comparison of all four types of ND5 genes in Hoplobatrachus, the protein structure of the two ND5 genes in H. rugulosus (BT) differed from those in H. rugulosus (WT) and H. tigerinus, which suggested that H. rugulosus (BT) diverged from H. rugulosus (WT) in evolution. According to the different nucleotide divergences, the transmembrane domains, three-dimensional structure, and gene arrangement of the ND5 gene between H. rugulosus (WT) and H. rugulosus (BT), H. rugulosus (WT) and H. rugulosus (BT), two different duplication of ND5 genes may have been independently formed in phylogenetic evolution.
Based on the mt genome divergence of the total nucleotide and protein-coding genes (Table 4), as well as the character of ND5 gene in H. rugulosus (BT), we can also conclude that H. rugulosus is a cryptic species complex. Yet, we failed to find the high genetic divergence between H. rugulosus (BT) and H. rugulosus (WT) using NCX1, Rag1, Rhod, and Tyr genes. mt genome rearrangement in H. rugulosus (BT). In all published dicroglossid sequences, tRNA Leu (CUN) , tRNA Thr , tRNA Pro , and tRNA Phe genes were translocated from their original position of Archaeobatrachia and formed a LTPF tRNA gene cluster the upstream of the 12S rRNA gene. However, the members and arrangement of this tRNA gene cluster are slightly modified in some taxa [21,54]. In our study, the tRNA Leu(CUN) gene was found between the duplicated D-loop-ND5 regions as previously observed in H. tigerinus [21] and H. rugulosus (WT) [22]. A TPF tRNA gene cluster at the upstream of 12S rRNA gene was formed, which was consistent with H. tigerinus [21] and H. rugulosus (WT) [22]. We analyzed the four samples of Thailand tiger frogs by amplifying sequences from ND5 to 12S RNA genes, and the TPF tRNA gene cluster was also found. Thus, the TPF tRNA gene cluster can be regarded as a synapomorphy of Hoplobatrachus, which was consistent with the finding of Alam et al. [21]. Zhou et al. [5], Alam et al. [21], and Chen et al. [55] suggested that the tandem duplication of tRNA Met gene can be regarded as a synapomorphy of Dicroglossinae. In our results, tandem tRNA Met genes were found to exist in H. rugulosus (BT).
Identical tandem duplication of ND5-D-loop was observed in H. tigerinus and H. rugulosus (WT), but only a similar tandem duplication of ND5-D-loop was observed in H. rugulosus (BT). The possible rearrangement pathway of the two different ND5 genes in H. rugulosus (BT) can be explained by the random duplication model as well as H. tigerinus [21]. Yet, the random loss has not happened in two ND5-D-loop regions. The mutation in two ND5 genes of H. rugulosus (BT) maybe happened dependently by some unknown selection pressures during the evolution history. The different arrangement pathways of the ND5 gene and the three-dimensional structure of ND5 proteins in H. rugulosus (BT) also suggested that H. rugulosus represents a cryptic species complex comprising at least two well supported lineages (BT and WT).

mt genome rearrangement and phylogeny of Dicroglossidae
The phylogenetic relationships among dicroglossids were mostly similar to the previous molecular phylogeny [5,52,56]. Although the phylogenetic topology is the similar to Zhang et al [52], the bootstrap values of MP, ML and NJ and the posterior probabilities of BI are well supported (Fig 5)In the complete mtDNA of dicroglossids retrieved from GenBank, ND5 gene translocation was observed in O. martensii, Fejervarya limnocharis, Fejervarya cancrivora, E. hexadactylus, H. rugulosus (WT), and H. tigerinus. In this study, the location of two ND5 genes of H. rugulosus (BT) was consistent with H. tigerinus and H. rugulosus (WT) [21,22]. We identified all different mt genomes of anurans in GenBank. ND5 gene translocation occurred in three distinct lineages of focusing on families closely related to Dicroglossidae: the lineage to O. martensii (Occidozyginae), the lineages to Mantellidae and Rhacophoridae, as well as the lineages to Fejervarya, Euphlyctis, and Hoplobatrachus (Dicroglossinae).
Generally, mt genome arrangements are believed to reflect phylogenetic relationships [57][58][59]. Using such information, three major lineages can be separated in Dicroglossidae: translocation of ND5, LTPF tRNA gene cluster, and tandem duplication of tRNA Met , which were observed in O. martensii (Occidozyginae). The LTPF tRNA gene cluster and the tandem duplication of tRNA Met were observed in Quasipaa, Nanorana, and Limnonectes. The translocation of ND5, the modified LTPF tRNA gene cluster, and the tandem duplication of tRNA Met were observed in Fejervarya, Hoplobatrachus, and Euphlyctis. The evolutionary relationships of dicroglossid taxa indicated in the phylogenetic trees, which was based on the concatenated sequences of the 11 protein coding genes (Fig 4), were similar to the traditional classification [14,54,[60][61][62][63]. Based on the results of genome rearrangement and phylogenetic relationship in dicroglossids, we found that the genome rearrangement of Dicroglossidae was consistent with the results of phylogenetic analyses (Fig 4).
Although NJ analysis based on the data of Rag1, Rhod, and Tyr genes, two clades (clade 1 and 2) in Hoplobatrachus (S1 Fig) were supported, the genetic divergences of those genes is low. Our findings highlight the need for further nuclear gene and crossbreeding studies of H. rugulosus (WT) and H. rugulosus (BT). Whether cryptic species exist or not in H. rugulosus, Thailand tiger frogs (H. rugulosus (BT)) should be prevented from escaping into wild environments to protect Chinese tiger frogs. We strongly advocate that Thailand tiger frogs in the farms should be strictly managed and that release of these frogs into the wild should be prohibited lest they produce hybrids with Chinese tiger frogs.

Conclusion
Based on the Cyt b gene divergence (13.8%) between H. rugulosus WT and BT, the mt genome divergence of the total nucleotide (14.0%) between H. rugulosus WT and BT, different three-dimensional structure of ND5 proteins and different rearrangement pathways of the ND5 gene between H. rugulosus WT and BT, all suggest that H. rugulosus is a cryptic species complex. Although using the nuclear gene (NCX1, Rag1, Rhod, and Tyr genes), we failed to find the high genetic divergence between H. rugulosus WT and BT, we found two clades (clade 1 and 2) in Hoplobatrachus (S1 Fig) using NJ analysis based on the data of Rag1, Rhod, and Tyr genes. Because of the genetic difference between H. rugulosus WT and BT, we suggest that we should prevent Thailand tiger frogs from escaping into wild environments lest they produce hybrids with Chinese tiger frogs.  Table. List of species used in this study, along with GenBank accession numbers and A+T content of total mitochondrial genome and control region (D-loop). (DOC) S3 Table. List of the nucleotide genetic divergences among 19 samples of NCX gene in species used in this study. The meaning of symbols is shown in S2 Table. (XLS) S4 Table. List of the nucleotide genetic divergences among 21 samples of Rag 1 gene in species used in this study including H. rugulosus (HM163612) and H. occipitalis (HM163613). The meaning of symbols is shown in S2 Table. (XLS) S5 Table. List of the nucleotide genetic divergences among 21 samples of Rhod gene in species used in this study including H. rugulosus (AJ564731) and H. tigerinus (AB489039). The meaning of symbols is shown in S2 Table. (XLS) S6 Table. List of the nucleotide genetic divergences among 21 samples of Tyr gene in species used in this study including H. tigerinus (AB277358) and H. occipitalis (AJ564729). The meaning of symbols is shown in S2 Table. (XLS) Natural Science Foundation of China (No. 31172116 and 31472015). We also thank the Academic Editor, Dr. Don Colgan and the anonymous reviewers for their important comments on the manuscript.