Species Delimitation and Evolutionary History of Tree Frogs in the Hyla chinensis Group (Hylidae, Amphibian)

Species are the cornerstone in many domains of biology research, which make accurate species delimitation critically important. In this study, the systematics and biogeography of the Hyla chinensis group were analyzed based on phylogeny, species delimitation, and ancestral area reconstruction methods. The phylogenetic results showed that six specific clusters existed in the H. chinensis group. Bayesian Phylogenetics and Phylogeography (BPP) analysis indicated that six distinct species exist due to the high probability values (>0.95), which were also supported by the Bayes factor (BF) analysis. The divergence time of the H. chinensis group was estimated to date back to 18.84 million years ago (Mya) in the early Miocene. Combining the results of ancestral area reconstruction, the H. chinensis group might have originated from Guangxi-Hainan, then spread eastwardly and reached Nanling Mountains, Wuyi Mountains, Huangshan Mountain, and Taiwan. In right-about colonization, it was gradually extended to the Yunnan-Guizhou Plateau, Sichuan Basin, Qinling Mountains, and Dabie Mountains. Considering the geological movement from early Miocene to Pliocene, the colonization pattern of the H. chinensis group may be closely related to the progressive uplift of the Qinghai-Tibetan Plateau (QTP) and historical climate change. Our study provided evidence for species delimitation and speciation process within the H. chinensis group. Our study supported the hypothesis that the evolutionary divergence in this species group was a consequence of the progressive uplift of the QTP and environmental change.


INTRODUCTION
For biogeography, abiotic factors (e.g., climate changes and tectonic events) and biological factors (e.g., interspecific or intraspecific interactions, competition, and predation) act as the major drivers temporally and geographically for biological evolution and diversification (Benton, 2009). Generally, for mountainous landscapes, the interactions of those factors provided beneficial conditions for the various microhabitats. Herein, those species endemic to mountain habitats often exhibit special phylogeographic patterns, such as the relatively small populations with well-defined geographical boundaries (Shepard and Burbrink, 2011;Huang et al., 2017;Pan et al., 2019). In southern China, many mountains (e.g., Hengduan Mountains, Qinling Mountains, Daba Mountains, Wuyi Mountains, and Dabie Mountains) are scattered, which form potential spatially isolated sky islands, providing various microhabitats with beneficial conditions for the speciation process of endemic species (Gao et al., 2015;Zhen et al., 2016). For example, due to the various microhabitats under climate and tectonic events, the Qinghai-Tibetan Plateau (QTP) had significant influence on the evolution of many animal groups (Päckert et al., 2012;Favre et al., 2015).
On the other hand, Li et al. (2015) had demonstrated that the Hyla originated from North America, then diffused to China via Beringia during the Middle Eocene to Early Oligocene (Smith et al., 2005;Wiens et al., 2006), which may be inferred that the speciation of H. chinensis group may be from northern China to southern China. However, the phylogenetic tree in the study by Li et al. (2015) disclosed that the base clades of H. chinensis group were all located in southern China, which may be a hint of another expansion route of the H. chinensis group.
Using genetic data and multiple analysis methods to solve taxonomic uncertainties enables us to disclose phylogenetic topology and speciation process. Here, we reveal a phylogeny of the H. chinensis group based on multiple mitochondrial and nuclear genes covering currently described species or subspecies within the H. chinensis group (Li et al., 2015). On the basis of species delimitation methods, we aim to clarify systematic and taxonomic matters bound up with species within the H. chinensis group. Meanwhile, we evaluate whether orogeny and climate oscillations affected the speciation and evolutionary history of H. chinensis group.

Ethics Statement
In this study, the sample collection of H. tsinlingensis and H. chinensis was conducted by a long-term investigation project on amphibian diversity in Dabie Mountains and Huangshan Mountain. This investigation project and sample collection were approved by the Anhui Normal University Academic Ethics Committee, Anhui Province, China.

Taxon Sampling
Based on previous study, we embraced almost all currently recognized species (76 individuals) within the H. chinensis group (Li et al., 2015) and chose two species (Hyla arborea, Hyla orientalis) as outgroups. Additionally, our own specimens (17 H. tsinlingensis individuals and two H. chinensis individuals) were collected from Dabie Mountains and Huangshan Mountain during 2011 to 2014, all samples were non-invasive sampling and the specimens were stored in School of Life Sciences, Anhui University, China (Figure 1). Details on specimen vouchers and GenBank accession numbers and specimen sites are listed in Table 1.

Laboratory Methods
The proteinase K digestion and phenol/chloroform extraction method were used to extract total genomic DNA (Sambrook et al., 1989). For combined previous sequence data (Li et al., 2015), the same genes were selected based on published primers and new primers (Supplementary Table S1), including four mitochondrial genes (12S ribosomal small subunit gene/12S rRNA, NADH dehydrogenase subunit 1 gene/ND1, tRNA-Leu, and the partial 16S ribosomal large subunit gene/16S) and one nuclear protein-coding gene (proopiomelanocortin A/POMC) .
All PCRs were performed within the same conditions in 30 µl volume: 10-40 ng of genomic DNA, 15 µl 2 × EasyTaq PCR SuperMix polymerase (containing 1 U Ex Taq, 0.4 mM dNTP, 3 mM Mg 2+ ; TransGen Biotech) and 0.2 µM of primers. PCRs were performed by the following protocol: an initial denaturing step of 5 min at 94 • C, followed by 32 cycles with denaturing 30 s at 94 • C, annealing 30 s at 50 • C and 55 • C (for mitochondrial gene and nuclear gene, respectively), extending 40 and 100 s (for mitochondrial gene and nuclear gene, respectively) at 72 • C, and a final extension step of 10 min conducted at 72 • C. PCR samples were checked on a 1% agarose gel. Subsequently, PCR products were purified by EasyPure PCR Purification Kit (TransGene), and each fragment was sequenced in both directions on the ABI 3730 semiautomated sequencer (PE Applied Biosystems). sequences of all genes (Burland, 1999) with default parameters, and the nucleotide sequences were checked by eyes. All the genes were concatenated for analysis and aligned in MEGA v.7.0 (Tamura et al., 2011). Aligned sequences had a total length of 2,474 bp (12S rRNA, 815 bp; 16S+tRNA+ND1, 1,172 bp; POMC, 487 bp). Two datasets were applied in phylogenetic analyses: (1) a data set consisting of the combined mtDNA genes (12S rRNA+16S+tRNA+ND1) was used to conduct species tree, Bayes factor (BF) delimitation (BFD) analyses (Sullivan and Joyce, 2005), infer divergence times, phylogenetic network, and genetic distance analysis; (2) the entire set of mitochondrial and nuclear genes (12S rRNA+16S+tRNA+ND1+POMC) was used to conduct the phylogenetic reconstruction [maximum likelihood (ML); Bayesian] and Bayesian Phylogenetics and Phylogeography (BPP) analysis (Rannala and Yang, 2003;Yang and Rannala, 2010).
Before phylogenetic analysis, the software jModeltest v.2 (Darriba et al., 2012) was used to find the best-fit nucleotide substitution model of each gene using Bayesian information criterion (BIC), and this optimal model (GTR+G, 12S; GTR+I+G, 16S+tRNA+ND1+POMC) was selected and implemented in all downstream analyses. Bayesian phylogenetic analysis was performed on different partitions of mitochondrial and nuclear datasets with a mixed-model approach separated into using MrBayes v.3.2.2 (Ronquist and Huelsenbeck, 2003). The homologous sequence of H. arborea and H. orientalis was used as outgroups. Two independent runs of Markov Chain Monte Carlo (MCMC) analyses for 10 million generations were conducted. The run was sampled every 1,000 generations, and 10% of the initial samples were discarded as "burn-in." The ML tree was generated with RAxML v.7.0.3 (Stamatakis, 2008) using the GTR model for mitochondrial and nuclear datasets. Support of nodes was calculated with 1,000 bootstrap replicates with the fast bootstrapping algorithm. Aside from the above analysis, we also operated "net between putative species mean distance" between the H. chinensis group species with 1,000 bootstrap replicates by Kimura two-parameter model on mitochondrial genes in MEGA v.7.0 (Tamura et al., 2011).

Divergence Time Estimation
Mitochondrial genes were used to estimate divergence times among H. chinensis group in BEAST v.1.8.0 . An MCMC approach with uncorrelated lognormal relaxed molecular clock for rate variation was set. Two independent runs were performed, consisting of 10 million generations, each run sampling every 1,000 generations with a burn-in set to 10% of  "−" represents no molecular data. The distribution areas were consistent with that within Figure 3.
Frontiers in Ecology and Evolution | www.frontiersin.org the samples. Tracer v.1.6 was used to check the stationarity of results (Rambaut and Drummond, 2007). TreeAnnotator v.1.8.0 (Rambaut and Drummond, 2007) and FigTree v. 1.4.2 (Rambaut, 2014) was used to annotate and visualize tree information. In the absence of appropriate fossils, we selected several calibration points information from previous work (Li et al., 2015), assuming a normal distribution for the divergence time between H. arborea group and the H. chinensis group, with a mean of 23 million years ago (Mya) and standard deviation of 3.04 (18-28 Mya).

Bayes Factor Delimitation
The BF is a common species delimitation selection tool in phylogenetics (Sullivan and Joyce, 2005) based on the marginallikelihood estimates (MLE) via path-sampling (PS) and steppingstone sampling (SS) analyses Xie et al., 2011;Li and Drummond, 2012). The scopes of BF are as follows: 0 < BF < 2 is not worth more than a bare mention, 2 < BF < 6 is positive evidence, 6 < BF < 10 is strong support, and BF > 10 is decisive .
Coupled with the former studies (Hua et al., 2009;Li et al., 2015) and the above phylogenetic analysis inference, six ingroup species in the H. chinensis group were assumed and four species delimitation scenarios (True, Lump, Split, and Reassign) were generated to disclose the inner species number in * BEAST (Heled and Drummond, 2010). For "True" scenario, individuals were assigned to six ingroup species in the H. chinensis group as prior set. For the "Lump" scenario, we inferred that two ingroup species were regarded as a single species, corresponding to the ingroup number of species from six to five. In contrast, the "Split" scenario suggested two ingroup species each split into two species, which indicated the total number of ingroup species from six to eight. As to the "Reassign" scenario, a total of three individuals were "incorrectly" reassigned to different ingroup species than in the "True" tree. PS and SS analyses were each run, totaling 10 8 generations with a chain length of 10 6 generations for 100 path steps.

Bayesian Phylogenetics and Phylogeography
In contrast to the results of our BFD method and to a commonly used method of species delimitation, we performed species delimitation analysis with the phased dataset for the two mitochondrial loci and one nuclear locus implemented in BPP v.3.0 (Rannala and Yang, 2003;Yang and Rannala, 2010). This method utilizes a reversible jump MCMC (rjMCMC) algorithm to calculate the posterior probabilities to speciation events that contain more or less lineages (Yang and Rannala, 2010). Between all BPP analyses, probability values ≥0.95 were considered a strong support in favor of a speciation event (Leaché and Fujita, 2010). The guide tree was generated from the species tree. The priors of ancestral population size (θ) and root age (τ) are directly related to the posterior probabilities of each result for the models. For example, the combination of large values for θ and small values for τ is assumed to be the most conservative, leading to a lower number of speciation events (Leaché and Fujita, 2010;Yang and Rannala, 2010). We evaluated three schemes of the prior of the θ and τ: (1) θ = G (1:2,000) and τ = G (1:10), (2) θ = G (1:2,000) and τ = G (1:100), (3) θ = G (2:2,000) and τ = G (1:10). The parameters of the rjMCMC analyses were set as 500,000 generations with sampling every 50 steps and 100,000 burn-in steps.

Species Tree Inference
The coalescent-based method implemented in * BEAST was used to construct the species tree (Heled and Drummond, 2010) based on mitochondrial genes. Two independent runs of 20 million generations were conducted. The sample frequency was set as 10,000 generations, and 10% of the total samples were discarded as burn-in. The other models and prior specification applied were set as follows: the nucleotide substitution model: HKY; Relaxed Uncorrelated Lognormal Clock (estimate); Yule process of speciation; random starting tree; alpha Uniform (0, 10). The convergence was checked by examining trace plots and histograms in Tracer. Runs were combined using LogCombiner. In addition, we tended to construct an uncorrected p-distances phylogenetic network with heterozygous ambiguities averaged and normalized by Splitstree v. 4.13.1 (Huson and Bryant, 2006). The neighbor-net ordinary least squares variance and equal angle algorithm were used, and 1,000 bootstrap replicates were calculated to assess branch support.

Ancestral Area Reconstruction
Geographical regions were delimited in terms of the current distribution area of the sequenced species of the H. chinensis group, at the same time, the information coming from the relevant literatures (Tang and Zhang, 1984;Li and Yang, 1985;Fei et al., 2009;Frost, 2014). The five areas were as follows: N, the southern China (Guangxi-Hainan provinces), which is a main distribution area of H. zhaopingensis; W, Eastern China; S, the southern Guangxi province in China, which is an important distribution range about H. sanchiangensis; Y, the eastern of the Tibetan Plateau (Yunnan-Guizhou Plateau and Sichuan Basin); Z, the Tsinling-Dabie Mountains (Figures 1, 3). LAGRANGE is a deep-time biogeographical model-based method that allows the incorporation of paleogeographical data (Ree and Smith, 2008;Chacón and Renner, 2014). Taking the effect of the LAGRANGE model components into account, we designed experiments that transform the adjacency matrix, hence, resulting in a total of three experiments (M0, M1, and M2). This is according to the assumption that the H. chinensis group, like all organisms, has a lower possibility to disperse over non-adjacent areas than adjacent areas. For this reason, no ranges were forbidden for M0; WZ, SZ, and NZ were forbidden for M1; WZ, NZ, SZ, and NY were forbidden for M2. To select the optimal model, we compared their log-likelihood (the data presented by LAGRANGE); meanwhile, we used the standard cutoff value of two log-likelihood units as indicating a conspicuous imbalance between models, with the less negative likelihood being preferred (Chacón and Renner, 2014). Ancestral areas were reconstructed by dispersal-extinction-cladogenesis model (Garzione et al., 2000) as carried out in the software Lagrange v.20130526 (Ree and Smith, 2008), and the chronogram obtained in BEAST was the starting component of the analyses.

RESULTS
Phylogenetic analysis of concatenated sequences (mtDNA data and nuclear gene) recovered a well-resolved tree with six major clades (labeled A-F) within the H. chinensis group (Figures 2, 3). Clade A corresponds to H. tsinlingensis and H. annectans chuanxiensis, which are mainly located in the Qinling-Dabie mountains; Clade B included H. annectans, H. a. wulingensis,  and H. a. jingdongensis, while H. a. gongshanensis and H. a. tengchongensis are within clade C, and they are all distributed in the Yunnan-Guizhou Plateau and Sichuan Basin. The remaining clade D (i.e., H. sanchiangensis), E (i.e., H. chinensis), and F (i.e., H. zhaopingensis) are located in the Guangxi province, Hainan province, and the Eastern China, respectively (Figures 2, 3). The phylogenetic network of the H. chinensis group showed the consistent groupings compared with the phylogenetic methods (Figure 4). Dating analyses indicated that the most recent common ancestor (MRCA) of the H. chinensis group dates back to 18.84 Mya [95% of the highest posterior density (HPD), 19.50-17.18 Mya] in the mid-Miocene. The divergence time between clades within the H. chinensis group was taken place from late Miocene (11.88 Mya) to the early Pleistocene (around 4.82 Mya) (Figure 3).
The BFD based on PS (BF, 12.62) and SS (BF, 20.84) models decisively supported six species in the H. chinensis group ( Table 2), corresponding to the six clades disclosed by the phylogenetic tree (Figures 2, 3). BPP methods also supported six separated species due to higher than 0.95 probability values ( Table 3). Species tree, consistent with BPP tree topology, recovered a concordant, robust phylogenetic topology (Figure 5). Pairwise sequence divergence (p uncorrected distance) between hidden species in the H. chinensis group ranged from 2.1% (A vs. B) to 11.4% (E vs. F) ( Table 4).
In the ancestral area reconstruction, the best model for the H. chinensis group was M2, which supported that it was dispersed from southern China to the Qinling-Dabie mountains and to the Eastern of the Tibetan Plateau was restricted ( Table 5). The analyses supported the Southern China (Guangxi-Hainan provinces, Area N) and Eastern China (Area W) as the ancestral area of the H. chinensis group, and most speciation events were attributed to the progressive uplift of the Himalayas (Figure 3 and Supplementary Table S2). Additionally, H. tsinlingensis originated from the Eastern of the Tibetan Plateau (Yunnan-Guizhou Plateau and Sichuan Basin, Area Y).

DISCUSSION
The phylogenetic analysis identified all the individuals of the H. chinensis group formed into six genetically distinct population clusters (i.e., Clades A-F) (Figures 2-4). Based on BF and BPP methods, the species delimitation suggested that these six genetically distinct clades could be regarded as hidden separated species in the H. chinensis group, which also got the support from genetic distance ( Table 4). In brief, clade A corresponds to H. tsinlingensis and H. annectans chuanxiensis; clade B included H. annectans, H. a. wulingensis, and H. a. jingdongensis; while H. a. gongshanensis and H. a. tengchongensis are within clade C. The remaining clades (D-F) corresponded to three putative species (H. sanchiangensis, H. chinensis, and H. zhaopingensis), respectively (Figures 2, 3), which were supported by the various species delimitation analyses. Overall, compared with the study of Li et al. (2015), some minor differences exist: H. annectans chuanxiensis belongs to H. tsinlingensis, not H. annectans; two subspecies of H. annectans (H. a. gongshanensis and H. a. tengchongensis) were regarded as separated species. As for the species not involved (H. simplex), we cannot give corresponding conclusions due to missing key samples in our research. As mentioned by previous researchers, further clarification of the relationship of H. simplex and H. zhaopingensis can be resolved only based on the H. simplex from its type locality in Vietnam and in-depth analysis (Li et al., 2015).  (Zachos et al., 2001(Zachos et al., , 2008. (A-F) represents the six clades within the H. chinensis group.    In addition, for the distribution pattern of the clade D, this clade includes samples from Anhui and Guangxi; it seems that a large geographic gap exists between these places, which Hyla is a small, arboreal, and semiaquatic frog; prefers to live in warm and damp environments, which widely inhabits boscages, paddy fields, or edges of rivers; and breeds in still water in ponds or paddy fields (Fei, 2005). During this period, the southern China humid climate is conducive to the survival of the species. For example, paleobotanical data indicated that the southeastern of the QTP has a warm and humid climate, was dominated by subtropical vegetation during the Miocene (Jacques et al., 2011), which had provided an opportunity for the first stage of speciation in the H. chinensis group. The second stage of speciation in the H. chinensis group occurs in the Southwest of China (Yunnan Province and Sichuan Province) and the Qinling Mountains-Dabie Mountains in China from the late Miocene to Pliocene (5.57∼4.82 Mya). During the Late Miocene to Pliocene, the progressive uplift of the QTP particularly at its eastern and northern margins (mostly province of Yunnan, Sichuan, and Qinghai) led to the formation of some rivers; the Hengduanshan hot spot of biodiversity was composed of those areas (Favre et al., 2015). In addition, the upheaval of the QTP had a significant impact on the atmospheric circulation in Asia and promoted the development of the Asian monsoon system (Kutzbach et al., 1993;Tang et al., 2013). The East Asian Monsoon system has controlled China's climate at that time, and this condition brought moisture air from the ocean to East China (Liu et al., 2009). Combining the geological events, they may be contributed to the second stage of speciation in the H. chinensis group.
In conclusion, although a previous study had demonstrated that the Hyla originated from North America, then diffused to China via Beringia during the Middle Eocene to Early Oligocene (Smith et al., 2005;Wiens et al., 2006;Li et al., 2015), it did not mean that the speciation of the H. chinensis group may be from northern China to southern China. Under the framework of speciation time within the H. chinensis group, the rapid uplifting mountain ranges (the Tibetan Plateau and its adjacent mountain) formed a blocky orographic barrier for many endemic species (Favre et al., 2015), which also played an important role in the formation of the Asian monsoon system (Guo et al., 2008;Song et al., 2010;Tang et al., 2013). Additionally, three East Asian monsoon intensification periods (∼15 Ma, ∼8 Ma, and 4-3 Ma) (Wan et al., 2007;Molnar et al., 2010;Jacques et al., 2011) also had urged the formation of humid and warm climates in south China (Sun and Wang, 2005), which was favorable for geographical dispersal, especially for amphibians (Che et al., 2010;Wu et al., 2013;Yan et al., 2013;Ye et al., 2013). More dispersal events often means that these species had more opportunities for allopatric divergence, which greatly affected the high levels of inter-population genetic divergence and unique patterns of genetic structure (Che et al., 2010;Wu et al., 2013;Yan et al., 2013;Ye et al., 2013;Favre et al., 2015). Therefore, based on those results, unlike the history of the Hyla evolution (North America→Beringia→China) (Smith et al., 2005;Wiens et al., 2006;Li et al., 2015), we can infer that the speciation and diffusion in the H. chinensis group had been from Guangxi-Hainan provinces to Guangxi province and Eastern China, and then to the Yunnan-Guizhou Plateau and Sichuan Basin, finally spread to Qinling-Dabie mountains. The diversification and speciation in the H. chinensis group also may be related to the special geological deformations and the climatic history.

CONCLUSION
As one of the species complexes in Hyla, the determined species number in the H. chinensis group was full of competing. Until now, no research focuses on the species delimitation based on the genetic data. In this study, different species delimitation approaches revealed that multiple species exist in the H. chinensis group. These methods indicated that there are six distinct species (from Clades A to F, respectively) in this species group. The progressive uplift of QTP and climate change led to the dispersal progress and formation of hidden species diversity in the Hyla chinensis group. Nevertheless, diagnostic morphological characters and other ecological evidences are still needed to supply for providing the integrative revision of this species group.

DATA AVAILABILITY STATEMENT
The accession number(s) can be found in the table within this article.

ETHICS STATEMENT
The animal study was reviewed and approved by the Anhui Normal University.

AUTHOR CONTRIBUTIONS
BZ, JL, and XW conceived the study. PY and TP contributed to sample collection. PY, GW, TP, XK, IA, and WZ carried out the laboratory work. PY, GW, TP, and XK analyzed the data and wrote the manuscript with contributions from BZ, XW, JL, PY, and TP. IA corrected the language. All authors approved the final version of the manuscript and agreed to be held accountable for its content.