Geographical rang evolution of whipping frogs based on a phylogeny of genus Polypedates (Anura: Rhacophoridae)

The genus Polypedates has a wide range, covering eastern India, Southeast Asia, South China, and eastward to the Philippines and Japan. Because of their poor marine dispersal capabilities, they are ideal organisms to infer geological and environmental history. Based on a large amount of data set for the genus Polypedates that we compiled from a dozen previous studies and partial mitochondrial and nuclear genes collected in this study, we calculated thorough statistical biogeographical analyses. We can confirm the genus’ Sundaland origin and showed its northward dispersal into Indochina and mainland Asia from the Late Oligocene to Middle Miocene. And the Red River did not mediate species exchange between Southeast Asia and mainland Asia until the end of the Miocene, with a sudden onset of northward dispersal in several clades independently at that time. The widespread P. leucomystax dispersed from Indochinese Peninsula southward to Malay Peninsula, Sumatra and Java, afterwards to Borneo and Sulawesi. Our biogeographical analysis supports the view of a recent introduction of P. leucomystax to the Philippines and the Ryukyus as previously suggested. Our results indicate that climate fluctuations have a profound impact on species diversification within genus Polypedates . The Red River did not act as a diffuse filter barrier for the species exchange until the end of the Miocene. And the widespread P. leucomystax dispersed in Southwest possibly facilitated by the freshwater plume of the emergent South Sunda River System.

likely facilitated the colonization of its wide range.
During the last decades, research on this genus mostly concentrated on the molecular phylogenetic relationships among Polypedates species, and particularly on the taxonomic study of Polypedates leucomystax species complex [5][6][7][8][9][10]. Kuraishi et al. [4] conducted a comprehensive phylogenetic study of the P. leucomystax complex based on both mtDNA and nuDNA, and including divergence time estimation, but with a limitation of biogeographic patterns analysis. Polypedates species were once regarded as extraordinary widely distributed from India to eastern Indonesia [11][12][13]. But after a long time of taxonomic confusion, reclassification of the genus led to the exclusion of P. mutus, P.
impresus, P. braueri and P. megacephalus in south China [8] and Taiwan [14]. Brown et al. [5] suggested recent population expansion of P. leucomystax complex in Philippines and Sulawesi. This leads to the hypothesis that human-mediated dispersal between oceanic islands. Kuraishi et al. [15] investigated the origin of P. leucomystax on the Ryukyus, also concluding that this was facilitated by accidental human transport from Southeast Asia. Blair et al. [16] published a study on the evolutionary history of the P. leucomystax in south China and on the Indochinese Peninsula including estimation of divergence times and population genetic analyses, which suggested a northern origin and southward dispersal. Their result refuted that the demographic expansion of this species and the wide spread sympatry of lineages originated from human-mediated dispersal between insular and mainland populations, and found the Red River as a partial barrier for gene flow. Buddhachat et al. [6] revealed allopatric and sympatric speciation under climatic pressure in P. leucomystax complex.
However, samples of these studies were only from Sunda Islands or Indochina mainland. A complete biogeographic study for the genus Polypedates is still unavailable.
As a widely distributed treefrog group, the genus Polypedates occupied in south China, India, most islands of Southeast Asia and Japan. And their current distribution area includes the biodiversity hotspots of Himalayas, Hengduan Mountains, Indo-Burma, Sundaland, part of Wallacea and Philippines [17]. From previous research, it manifested that this genus diversed from Oligocene [18], and remains to be genetic divergence in the very close recent. Therefore, it seems to be that the dispersal process of Polypedates is in the wake of paleotectonism and paleoclimate events such as the Indian-Eurasian collision, the uplift of Himalayas, and the Pleistocene sea-level fluctuations from glacial and interglacial. It is still necessary to present convincing evidence based on previous achievement to make clear how geological events affected the evolution of these frogs, meanwhile mediately verify when the geological events occurred with the time estimate of the living Polypedates.
Here we compiled a large data set of currently recognized species within the genus Polypedates from the literature and partial mitochondrial and nuclear genes collected in this study to (1) investigate its phylogenetic relationships, reinterpreted the distribution of each species particularly P. leucomystax,

Phylogenetic analyses and divergence time estimates
Both Maximum Likelihood and Bayesian phylogenetic inference resulted are consistent ( Fig. 1 and Additional file 1: Fig. S1). Each major clade was strongly supported in both analyses. The species P. cruciger, P. maculatus and P. pseudocruciger which are all from India belong to the same clade with high support. Most species contained in P. leucomystax complex such as P. macrotis, P. mutus, P. braueri, P. impresus, and P. megacephalus clades were strongly supported as valid species. Species dispersal between India and other place happened twice, once with the species P. cruciger, P. maculatus and P. pseudocruciger, and once with the species P. teraiensis (Fig. 2).

Biogeography analysis
The two biogeographical analyses methods equally support a Sunda islands origin for the genus Polypedates (Fig. 2). However, more precisely, we could infer as the most likely scenario that P.   [19] which was abled from middle Oligocene to early Pliocene [20]. P. leucomystax dispersed to the Philippines very recently (Holocene); subsequently reaching Japan.
The MRCA P. leucomystax dated to the late Miocene and early Pliocene (95% HPD 6-3 Myr) with an ancestral range comprising Malay Peninsula, Sumatra and Java (Fig. 3). In a very recent epoch P.
leucomystax (samples 86-101) showed high gene flow between Borneo and the Malay Peninsula. The spread of P. leucomystax from Java to Sulawesi was detected only once, but to Sumatra twice, firstly from Indochina and second from Java. The Philippines' samples (104 and 106) along with the Japanese samples (105 and 107) nested within the Java clade, indicating that P. leucomystax in of the Philippines originated from Java and subsequently dispersed to Japan.

Discussion
Late Oligocene to Early Miocene origin of the Polypedates.
Amphibians are often considered excellent model system for biogeographic studies on overseas dispersal events [18,19,[21][22][23][24][25][26][27]. Many biogeographic studies suggested Southeast Asia is a cradle for amphibian speciation [21,22,25,26]. In this study, the phylogenetic analyses supported the sister relationship of Taruga and Polypedates. Taruga eques and Taruga fastigo were accommodated by Meegaskumbura et al. [28]. As Taruga in this study occurred in India and Rhacophoridae definitely origined in India [18], it can be inferred that the MRCA of genus Taruga and Polypedates originated from India. However, the biogeographical analyses of the present, compiled data set confirm Sundaland as the origin of the genus Polypedate, with subsequent dispersed to Indochina (Fig. 2).
These patterns elucidated that the MRCA of genus Polypedates diverged from the MRCA of genus Taruga probably in Sundaland, and differentiation dispersed to more northern regions, which in line with Kuraishi et al. [4]. On this occasion, the lineages from India in genus Polypedates should be considered further. Although locality specimens and sequences from India are lack in this study, it can be inferred that dispersal between India and Southeast Asia does not seen to be only once. P. cruciger, P. maculatus, and P. pseudocruciger dispersed from Sundaland to India  first time, coincided with the global Mid Miocene Climatic Optimum [29]. And P. teraiensis dispersed in sympatric regions in Indochina and India second time. Additionally, cryptical species distributed in India might indicate ancestor dispersal from India to Southeast Asia, such as other species [18,[30][31][32], which would change the original biogeographical history of genus Polypedates. Additional sampling from India is required to test these predictions.

Colonization from South to North
After originated in Sundaland, genus Polypedates colonized into the Indochina (14-8 Myr) and then dispersed northward to mainland Asia during the Pliocene and Pleistocene. During these periods, genus Polypedates underwent range expansions at different times. Similarity, a recent study of genus Leptobrachella, which is also distributed in Southeast Asia, showed that range expansions from Sundaland to mainland Aisa occurred in this genus from Mid-Miocene to Pliocene [21]. Hughes et al. [33] suggested the Isthmus of Kra unlikely to be significant biogeographical boundaries, and also indicated the climate zone in Indochina during Last Glacial Period is the explanation for species distinct distribution. We assumed that genus Polypedates expanded Indochina from Sundaland, leading to diversification of species, under the influence of climate changes.
During the Middle Miocene, no divergence could be detected between the Indochina (i.e., southwest of the Red River) and South China (i.e., northeast of the Red River). The Red River zone has a complex geological history [34,35]. The formation of the Red River between southern China and northern Vietnam, initiated at 27 Myr and terminated locally at about 21 Myr, not exceeding 17 Myr [36][37][38], which could act potential geographic barriers for species dispersal. During this period, tectonic movements did not cause species divergence. As many taxa occur only one side of the Red River (such as Sophora davidii in Fan et al. [39]), we could hypothesise that change in the runoff regime might have caused a stronger isolation effect of the Red River by then, with the reservation that a denser population genetic sampling of the respective species will be necessary to draw sound conclusions.
megacephalus. Our result of the unencumbered of the Red River from late Miocene to early Pliocene for the genus Polypedate was similarity to the viewpoint from Blair et al. [16], whereas it consummated his consequence with a more larger dataset. Blair et al. [16] found that the Red River forms a partial barrier for some lineages and cannot explain the overall pattern, most likely due to the indistinct taxonomy of P. leucomystax and a circumscribed lineage divide. In the view of several studies showed genetic diversification between species on both side of the Red River [41,42], Red River zone might conceivably function as diffuse filter barriers for the genus Polypedates.
Range evolution throughout oceanic islands within P. leucomystax.
From early Pliocene till Holocene, P. leucomystax dispersed back to Sundaland from the Indochina (Fig. 2) and spread to the Malay Peninsula, the Greater Sunda Islands, the Philippines and to Japan (Fig. 3). The north-to-south dispersal pattern was also consistent with a genealogy and demographic analysis [16]. Not surprising, P. leucomystax shows considerable genetic differentiation across its vast range, in line with previous studies [5]. In view of the peculiar paleoclimate and tectonics through the continent of Sundaland and nearby area in Pleistocene [43], the wide spread of P. leucomystax might be driven by eustatic sea-level fluctuations. The high gene flow between Borneo and the Malay Peninsula samples conformed to the contiguous terrestrial areas resulted from glaciation. The Pleistocene colonization of Sulawesi from Java, despite the absence of a continuous land connection between Sulawesi and Sundaland [44], is not surprising and fits to the frequent dispersal across Wallace's line in other organisms including freshwater crabs, frogs, reptiles and mammals [45].
Possibly, the freshwater plume of the South Sunda River System enabled dispersal on driftwood from Java to Sulawesi during an eustatic sea-level lowstand, as hypothesised for freshwater crabs [46]. P. leucomystax spread to the Philippines from Java and then reached Japan close to the Holocene.
Considering the massive distance and biogeographical barriers for amphibians that make a natural dispersal scenario highly unlikely, we support the hypotheses that transportation of agricultural products between islands facilitated the range expansion of P. leucomystax to the Philippines [5] and that accidental transport of few individuals with military cargo from somewhere around the Philippines established the populations on the Ryukyus as suggested by Kuraishi et al. [15]. We can not reject the hypothesis that Pleistocene sea-level dynamics have effects on the diversification of this species as an explanation for the range expansion with P. leucomystax. However, further evidence is required to address this possibility.

Conclusions
Over all, our skyline analyses revealed the evolutionary history of the genus Polypedates. The genus Polypedates mostly originated in insular Southeast Asia during the Late Oligocene. In the process of northward dispersal, climate changes influenced expansion and species diversification from Sundaland to Indochina. The Red River did not act as a diffuse filter barrier for the species exchange until the end of the Miocene. Since the early Pliocene P. leucomystax dispersed all over the Sundaland with Java as a source area. Our results make an effort to the combination of tectonics and climate changes and biological evolution with the genus Polypedates, corroborate the hypothesis about climate changes during Last Glacial Period and Red River fault zone with other fauna and flora. We also expect more samples from India to enrich the database, which may rewrite the biogeographic explanation within this group.

DNA extraction, PCR amplifications and sequencing
In order to researching the genus Polypedates, candidate DNA fragments, including three mtDNA and three nuclear DNA for 15 individuals from 6 species were newly acquired in this study. Genomic DNA was extracted from either muscle, or liver tissues initially preserved in either 95% ethanol. The Ezup Column Animal Genomic DNA Purification Kit (Sangon Biotech, China) was used for genomic DNA extraction following minor modifications to the manufacturer's protocol. We amplified and sequenced mitochondria gene (12S rRNA, tRNA Val and 16S rRNA) and three nuclear gene fragments with the indicated primer pairs, including proopiomelanocortin (POMC), Exon 1 of tyrosinase (RHOD) and exon 1 of tyrosinase (TYR) (Additional file 1: Table S2). Double stranded polymerase chain reaction (PCR) amplifications were conducted in total volume of 25 µl, including 12.5 µl Taq PCR Master Mix (2 X blue dye), and 9.5 µl ddH 2 O and 1 µl FS01 and 1 µl Rend (C = 10 µl/L), and 1 µl DNA templat. Cycle conditions for mitochondrial DNA were an initial denaturation of 3 min at 95℃, followed by 94℃ denaturation for 1 min, 55℃ annealing for 1 min, and for 1 min extension at 72℃. Final extension at 72℃ was conducted for 10 min. For TYR, RHOD, and POMC, the same procedure as for mitochondrial DNA was used, but with annealing at 56, 52, and 52℃, respectively.

Phylogenetic Analyses
The best partitioning scheme and model selection were estimated with PartitionFinder 2 [58] based on the Bayesian Information criterion. The data set was accordingly partitioned into two, a mitochondrial (12S rRNA gene, tRNA val , 16S rRNA gene) and a nuclear (TYR, RHOD, POMC) partition.
And the best fitting substitution model was GTR + I + G for the mitochondrial and TRNEF + I + G for the nuclear partition. We used two methods to construct phylogenetic trees, including Bayesian inference (BI) and maximum-likelihood (ML). BI was implemented in MrBayes 3.2.0 [59], the following settings were applied: number of Markov chain Monte Carlo (MCMC) generations = 5000000 and sampling frequency = 100. The first 12500 sampled trees were discarded as a conservative burn-in. ML phylogeny of the concatenated mitochondrial and the nuclear partition was calculated using RAxML 8.2.10 [60] under a GTR + Gamma model with 1000 bootstrap replicates; switching to the per-site rate category model.

Divergence time estimation
We used BEAST 2.5 [61] to estimate both a phylogenetic tree and divergence times. Site models were estimated with jModelTest 2.1.6 [62] based on the Akaike information criterion. And the best fitting substitution model was GTR + I + G for the mitochondrial and SYM + I + G for the nuclear partition. Few stratigraphic events can be associated with the divergence of species in Polypedates. Consequently, in the light of previously published data [18], we implemented second calibration points for the most recent common ancestor of Polypedates, Feihyla, Ghatixalus and Taruga, implemented as a Normal calibration density (Mean = 34.8, Sigma = 2.3) and the divergence point between genus Polypedates and Taruga, implemented as a Normal calibration density (Mean = 24, Sigma = 2.1). We applied a relaxed uncorrelated lognormal clock and a Yule tree prior. The analysis was repeated two times for 50 million generations separately, sampling every 5000th iteration on CIPRES Science Gateway (https: //www.phylo.org/portal2/login). We checked for stationarity of the Markov chain and potential autocorrelation (effective sample sizes > 200 for all sampled parameters) in Tracer v1.6 [63]. The first 25% of samples were discarded as burn-in, and the samples of both runs combined in LogCombiner v2.5. This file was used in TreeAnnotator v2.5 to identify and annotate the maximum clade credibility tree. Given the inherent uncertainty of divergence time estimates, we rounded the estimates (95% HPD interval boundaries) to the full million.

Estimation of ancestral areas
We conducted two kinds of biogeographical analyses (1) a likelihood method under DIVALIKE + J model (2) a Bayesian approach under the Bayesian Binary MCMC (BBM) method (Ronquist and Huelsenbeck, 2003) by the program RASP4 [64]. We defined the geographic areas based on established biogeographic and continental boundaries, such as the Red River fault between China and IndoChina Peninsula [41], the Isthmus of Kra between Indochina and the Malay Peninsula [65] (Fig. 2).
The DIVALIKE + J approach was conducted with the R-package "BioGeoBEARS" [66]. Firstly, we compared the six models adopted in the R-package "BioGeoBEARS" based on their AIC weight. Our ancestral range evolution analyses supported DIVALIKE + J as the model with the best fit to our data (LnL= -80.48, AIC = 167, see Additional file 1: Table S3). We set the maximum number of areas of a species' range to three. For the BBM method, we set the maximum number of areas was three, and used F81 as the states frequencies, and set Gamma distribution as among-site rate variation. Other parameters are default.
To investigate the natural history of insular P. leucomystax, we pruned the tree from BEAST maximum clade credibility tree and only left P. leucomystax samples with the R-package "Ape" [67], and applied the same DIVALIKE + J approach on it. The new defined geographic areas were: S-Sumatra, M-Malay Peninsula, B-Borneo, A-Java, P-the Philippines, J-Japan, W-Sulawesi (Fig. 3). Before the DIVALIKE + J analysis, we are grouping the samples with every clade and separated sample 85 from other samples in the same clade as one group, and sample 102, 103 from other samples in the same clade as one group additionally, then estimated the net evolutionary divergence between each group with MEGA7 [56] (Fig. 3). The result showed that the distance between group Ⅳ and Ⅴwere 0.01, Ⅵ and Ⅶ were 0.004 (see Additional file 1: Table S4), therefore we considered the species in each clade as a cryptic lineages and continue the following DIVALIKE + J analysis.

Declarations
Here we would like to thank Yunyun Lv and Claudio Sabatelli for technical assistance and insightful comments and suggestions.

Authors' contributions
J-T.L. and D-C.J. conceived the study; X-L.D. and L-M.Y. performed data analyses; X-L.D., L-M.Y. and S.K. prepared the initial manuscript draft. All authors worked on the final manuscript.

Availability of data and materials
The datasets analysed during the current study are available in the GenBank repository, https://www.ncbi.nlm.nih.gov/

Ethics approval and consent to participate
This study was performed in accordance with the approval of Experimental Animal Ethics Committee of Chengdu Institute of Biology, Chinese Academy of Sciences.

Consent for publication
Not applicable.   S1. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 3
Estimation of ancestral areas for the species Polypedates leucomystax based on the Bayesian maximum clade credibility tree given in Fig. 2. 95% credibility interval of the time 24 estimates was given in square brackets. The square labels on the tree are range estimates according to the DIVALIKE+J model. The Roman numerals are used for the group distance estimation. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. Collection sites of species in this study. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.