The reanalysis of biogeography of the Asian tree frog, Rhacophorus (Anura: Rhacophoridae): geographic shifts and climatic change influenced the dispersal process and diversification

Ethical statement

The collection of samples was performed within a long-term investigation project on amphibians of Dabie mountains. This investigation project and the sample collection were approved by the Animal Research Ethics Committee of Anhui University (Animal Ethics number: AHU3110) and Anhui Tianma National Nature Reserve, Anhui Province, China. Field experiments were approved by Anhui Tianma National Nature Reserve, Anhui Province, China.

Data collection

For the phylogenetic analyses, sequences of about half the species of Rhacophorus were used in combination with sequences of two outgroup species, Polypedates megacephalus (Rhacophoridae, Polypedates) and Spinomantis peraccae (Mantellidae, Spinomantis) (Li et al., 2013; Li et al., 2012a; Li et al., 2012b). Sequence data were obtained from GenBank (the GenBank Accession numbers are given in Table S1). In total, there were 149 individuals of 57 species of Rhacophorus involved (Fig. 1 and Table S1). All the taxonomic revisions within Rhacophorus were follow previous studies (Biju et al., 2013; Li et al., 2013; Orlov et al., 2012).

Sampling, DNA extraction, PCR amplification, and sequencing

Between 2012 to 2015, nine specimens of R. zhoukaiyae were collected from the Dabie Mountains, China (Pan et al., 2017). Muscle tissue from each individual was sampled and preserved in 100% ethanol for DNA extraction. Total DNA was extracted from the samples using a standard proteinase K/phenol-chloroform protocol (Sambrook, Fritsch & Maniatis, 1989). An EasyPure PCR Purification Kit (TransGene, Strasbourg, FR) was used to purify the DNA extractions. The sequences of 12S and 16S ribosome RNA (rRNA) of R. zhoukaiyae were collected from Pan et al. (2017). In addition, we also amplified and sequenced five nuclear gene fragments with the indicated primer pairs (Table S2), including brain-derived neurotrophic factor (BDNF), proopiomelanocortin (POMC), recombination activating gene 1 (RAG-1), rhodopsin exon 1 (RHOD) and tyrosinase exon 1 (TYR) (Bossuyt & Milinkovitch, 2000; Li et al., 2009; Van Der Meijden et al., 2007; Vieites, Min & Wake, 2007; Wiens et al., 2005). Polymerase chain reactions (PCR) were performed using a reaction mixture (25 µL) containing 1 µL genomic DNA (concentration 10–50 ng/µL), 2.5 µL 10× buffer, 1 µL of 2.5 mM MgSO4, 2 µL of 2 mM dNTPs, 1 U Taq polymerase (Meridian Bioscience, Singapore) and 0.3 mM of each of the primers. Pure molecular biology grade water was added to reach the appropriate volume. The amplification protocol included an initial denaturation step of 95 °C for 5 min; this was followed by 32 cycles of denaturation at 95 °C for 30 s, primer annealing at 51 °C –57 °C for 30 s, and an extension at 72 °C for 40 s–80 s, with a final extension at 72 °C for 10 min. PCR products were purified using an EasyPure PCR Purification Kit (TransGene) and sequenced using previous primers and the BigDye Terminator v3.0 Ready Reaction Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) following the manufacturer’s instructions on an ABI Prism 3730 automated sequencer. All the sequences obtained in this study were deposited into GenBank (Table S1). For the analyses, the sequences were trimmed to match data downloaded from GenBank, then all the sequences were aligned automatically using Clustal X version 1.83 (Thompson et al., 1997), followed by visual confirmation and manual adjustments. Nucleotide sites with ambiguous alignments were removed from the analyses, and gaps were analyzed as missing data.

Phylogenetic analyses

Two different datasets were generated for the different analyses. Dataset 1 was used for a phylogenetic analysis of Rhacophorus by Maximum Likelihood (ML) and Bayesian methods, and was comprised of the 12S and 16S rRNA gene together with the complete t-RNA for the valine sequence of the Rhacophorus species and the outgroups (Table S1). Dataset S2 contained more genes (12S, 16S, Val, BDNF, POMC, RAG-1, RHOD, TYR) of more individual and species than Dataset S1 (Table S1). However, it was only used to calculate a Bayesian consensus tree. The best-fit model of evolution was calculated by MrModeltest 1.0 b under the AIC criterion (Nylander, 2003). ML analyses were performed in RAxML version 8 (Stamatakis, 2014) and a general time reversible model of nucleotide substitution under the Gamma model of rate heterogeneity (i.e., GTRCAT). Support for the internal branches for the best-scoring tree was evaluated via the bootstrap test with 1,000 iterations. A Bayesian inference of phylogeny was performed using the MrBayes 3.1.2 software program (Huelsenbeck & Ronquist, 2005), using the best-fit substitution model. Two Markov Chain Monte Carlo (MCMC) models were run to provide additional confirmation of the convergence of posterior probability distributions. Analyses were run for 3,000,000 generations. Chains were sampled every 1,000 generations. The first 25% of the total trees were discarded as “burn-in” and the remaining trees were used to generate a majority-rule consensus tree and to calculate Bayesian posterior probabilities.

Divergence time analyses

To estimate divergence times of Rhacophorus, we applied a Bayesian MCMC method with mitochondrial genes (Dataset S1), which employs a relaxed molecular clock approach, as implemented in BEAST 1.7.4 (Drummond et al., 2012). We assumed a relaxed uncorrelated log normal model of lineage variation and a Yule Process prior to the branching rates based on the GTR + I + G model as recommended by MrModeltest 1.0 b (Nylander, 2003). Four replicates were run for 10,000,000 generations with tree and parameter sampling every 1,000 generations. The first 25% of samples were discarded as burn-in. All parameters were assessed by visual inspection using Tracer v. 1.5 (Rambaut & Drummond, 2007). The tree was generated and visualized with TreeAnnotator v. 1.6.1 (Rambaut & Drummond, 2010) and FigTree v. 1.3.1 (Rambaut, 2009), respectively. Calibration points were taken from Li et al. (2013) (Table 1). In addition, to visualizing the temporal accumulation of species, a log-transformed lineage-through-time (LTT) (Nee, May & Harvey, 1994) plot was constructed and compared with the null distribution for the LTT line simulated under the empirical pure-birth model. For visualizing diversification rate changes, we plotted the number of newly appearing species against the fixed time intervals of 2 million years (Ma) (Venditti, Meade & Pagel, 2010).

Ancestral area reconstructions

Ancestral area reconstructions were inferred by the program RASP 3.2 (Yu et al., 2015) for speciational evolution in phylogenetic trees, using the Bayesian Binary MCMC (BBM) method (Ronquist & Huelsenbeck, 2003) and the statistical dispersal-vicariance method (S-DIVA) (Yu, Harris & He, 2010). To reconstruct ancestral areas on the basis of the topography, the distributional range of Asian Rhacophorus was divided into four regions, W, X, Y and Z (Fig. 2). W represents Southeast Asia, including the Indochinese Peninsula, Sundaland and the south margin of the Tibetan Plateau, X contains the Hengduan mountains and the mountains around the Sichuan Basin, Y refers to South China and Japan and Z represents India (Fig. 2). The tree data sets and the condensed tree were generated by BEAST 1.7.4 (Drummond et al., 2012). The distribution of each species was collected from http://maps.iucnredlist.org. For all analyses, the maximum number of ancestral areas at each node was constrained to three. The frequencies of an ancestral range at a node were averaged over all trees and each alternative ancestral range at a node was weighted by the frequency of occurrence for the node.

Molecular phylogenetic analyses

The aligned mtDNA gene fragments from Rhacophorus consisted of 1,935 bp nucleotide positions before trimming (Dataset S1). After trimming, 1,851 nucleotide positions were retained for genealogical reconstructions. The fragments contained 934 constant and 917 potentially phylogenetically informative characters. The ML or BI phylogenetic approaches based on Dataset S1 resulted in virtually identical topology, and all terminal clades obtained relatively high-supporting values (Fig. S1). The genus Rhacophorus was supported as monophyletic containing four major clades (Fig. S1). For further probing of the dispersal process and diversification of the Asian tree frog, the molecular dating and ancestral area reconstructions were carried out. The phylogenetic tree, collected from the molecular dating, showed three distinct clades (A, B and C) in the genus of Rhacophorus (Fig. 2). There were some difference in the species distribution areas among the three clades. Species in clade A were mostly distributed in Southeast Asia and East Asia, species in Clade B were distributed in Southeast Asia and India, and species in lineage C only found in Southeast Asia. Clade A contained six groups, A1 to A6 (Fig. 2). The phylogenetic tree, based on Dataset S2, was largely consistent with the results from Dataset S1 (Fig. S2). However, there were some minor differences between them, such as the polyphyletic of clade B and C in Fig. S2 . But, generally, it did not affect the results of ancestral area reconstructions of Rhacophorus.

Molecular dating, ancestral area reconstructions and lineage through time

Dating analyses based on Dataset S1 suggested that the most recent common ancestor (MRCA) of Rhacophorus dates back to 29.51 Ma (median value; 95% of the highest posterior density [HPD] = 25.00–34.07 Ma) (Table 1 and Fig. 3).The MRCA of Clade A and Clade B was estimated at 27.38Ma (95% HPD = 22.44–32.17 Ma). The MRCA of Clade A was 21.56 Ma (95% HPD = 17.92–25.22 Ma) and the MRCA of Clade B was 26.73 Ma (95% HPD = 21.56–31.83 Ma).

Ancestral area reconstructions from S-DIVA and BBM analyses were largely similar with some minor differences (Fig. 2). All analyses supported Southeast Asia (Area W, Fig. 2) as the ancestral area of Rhacophorus and most speciation events were attributed to dispersal. The empirical LTT plot of Rhacophorus showed that, after a lengthy period of constant diversification, the diversification rate of the genus had increased during the middle Pliocene. The cumulative curve of species-birth per time interval showed that the diversification of Rhacophorus fluctuated through time, especially during 14–12 Ma and 10–4 Ma (Fig. 4).

The dispersal process of Rhacophorus and its spread toward East Asia and India

Previous studies have indicated that the diversification of Rhacophoridae was closely linked to the India-Asia collision (57 Ma–35 Ma) (Li et al., 2013). Southeastern Asia houses three globally significant hot spots divided by sharp, yet porous biogeographic boundaries (Evans et al., 2003; Favre et al., 2015; Schmitt, Kitchener & How, 1995; Wallace, 1860). Studies have shown that the dynamics of the formation of biodiversity in Southeastern Asia is assumed to be interrelated with many geological events and a unique climatic history. Events such as the continuing processes of volcanic uplift and the emergence of many new islands in Indo-Australian Archipelago during the Miocene-Pliocene (Fig. 3C) (Esselstyn, Timm & Brown, 2009; Hall, 1996; Hall, 1998; Hall, 2002; Lohman et al., 2011), the rapid uplifts of the Tibetan Plateau (Shi et al., 1999), repeated sea level fluctuations during the Pleistocene (Bird, Taylor & Hunt, 2005; Esselstyn, Timm & Brown, 2009; Hall, 1998; Heaney, 1985; Heaney, 1986; Jansa, Barker & Heaney, 2006; Voris, 2000) and the onset of the Asian monsoon system (An et al., 2001; Qiang et al., 2001; Sun & Wang, 2005; Zhisheng et al., 2001). Many phylogeographical studies of plants and animals support this assumption (Deng et al., 2011; Klaus et al., 2016; Shi et al., 1999), such as those on Lilium (Gao et al., 2013), Delphinieae (Jabbour & Renner, 2012), Hyoscyameae (Tu et al., 2010), Mandragoreae (Tu et al., 2010), Saussurea (Wang et al., 2009), birds (Lei, Qu & Song, 2014; Tietze & Borthakur, 2012; Tietze et al., 2013; Yang, Dong & Lei, 2009), Hynobiidae (Zhang et al., 2006), lizards (Guo et al., 2011) and Spiny Frogs (Che et al., 2010), so the diversification and speciation in Rhacophorus may also be related to the special geological formations and the climatic history.

The phylogenetic analysis shows that Rhacophorus is composed of multiple lineages. In the phylogenetic tree with timescale, calculated by BEAST, Rhacophorus is composed of three major clades, A, B and C (Fig. 2). Among these clades, Clade C was the basal branch of Rhacophorus, which contained ten species from Southeast Asia, and the age of the MRCA of Rhacophorus was estimated at 29.51 Ma (i.e., 95% CI [25–34.07 Ma], Fig. 2 and Table 1). The MRCA of Clades B and Clade A was 27.38 Ma (95% CI [22.44 Ma–32.17 Ma]) during the Oligocene (Fig. 2, Table 1). The members of Clades A and B are mainly distributed in the south of the Tibetan Plateau margin, India and Eastern Asia (Fig. 1). Clade A contained six groups which were distributed in three areas: Southeast Asia (group A1), the south of the Tibetan Plateau margin (group A2) and an Eastern Asia (group A3 to A6) (Figs. 1 and 2). The MRCA of Clade A occurred 21.56 Ma ago (95% CI [17.92–25.22 Ma]; Fig. 2) and the time of the split of different groups was estimated at 14.09 Ma (A2 vs A3∼A6, 95% CI [10.96–17.41 Ma]), 11.39 Ma (A3 vs A4∼A6, 95% CI [8.89–14.16 Ma]), 8.56 Ma (A4 vs A5∼A6, 95% CI [6.43–10.88 Ma]) and 5.33 Ma (A5 vs A6, 95% CI [3.92–6.99 Ma]) respectively (Table 1). In addition, the LTT plot analysis indicated an increased diversification rate during two periods (14–12 Ma and 10–4 Ma) (Fig. 4). Basically, the above mentioned phylogeographical information reflected the trend of diversification and the speciation process. Obviously, the distribution of these species expanded continuously from Southern Asia to India and Eastern Asia, reaching as far as Japan (Fig. 2).

During the Oligocene and Miocene the uplift progressed, causing the extension of the Tibetan Plateau (Harrison et al., 1992; Mulch & Chamberlain, 2006). The start of the uplift of the northern Tibetan Plateau occurred at about 30 Ma BP (Sun & Wang, 2005) or slightly earlier (Wang et al., 2012b). Then, the eastern parts of the Tibetan Plateau likely reached an elevation comparable to the present-day elevation in the Mid to Late Miocene (from 15 to 5 Ma) (Axelrod, 1997; Currie, Rowley & Tabor, 2005; Jacques et al., 2011; Spicer et al., 2003; Tapponnier et al., 2001; Valdiya, 1999; Zhang et al., 2013). The southeastern edge of the Tibetan Plateau, the Hengduan mountain range, experienced rapid uplift only after the Miocene (5.33 Ma), reaching a peak elevation shortly before the Late Pliocene (5.33–2.66 Ma) (Li & Fang, 1999; Mulch & Chamberlain, 2006; Sun et al., 2011; Zheng et al., 2000), which separated several major rivers that ran in parallel (the Yangtze, Mekong, and Salween valleys) (Clark et al., 2004). This series of rapid Tibetan Plateau uplifts dramatically changed the terrain and landform in this area, which resulted in speciation, especially in animal groups (Che et al., 2010; Deng et al., 2011; Gao et al., 2013; Jabbour & Renner, 2012; Lei, Qu & Song, 2014; Li et al., 2013; Shi et al., 1999; Tietze & Borthakur, 2012; Tietze et al., 2013; Tu et al., 2010; Wang et al., 2009; Zhang et al., 2006). Zhang et al. (2006) found that the origin and phylogenetic divergence of the Hynobiidae had a correlation to the uplift of the Tibetan Plateau (Zhang et al., 2006). The phylogenetic history of Paini (Anura: Dicroglossidae) illuminates a critical aspect of the timing of geological events, especially for the uplift of the Tibetan Plateau (Che et al., 2010). On the other hand, the Tibetan Plateau and its adjacent mountain ranges acted as an orographic barrier to atmospheric circulation in Asia and consequently contributed to the formation of the Asian monsoon system, which was one of the major climatic changes in this region (Early Miocene, 24 Ma) due to the Tibetan Plateau’s considerable size and altitude (Guo et al., 2008; Kutzbach, Prell & Ruddiman, 1993; Liu & Yin, 2002; Ruddiman & Kutzbach, 1991; Song et al., 2010; Sun & Wang, 2005; Tang et al., 2013; Zhisheng et al., 2001). In the following millions of years, the East Asian monsoon intensified three times (∼15 Ma, ∼8 Ma and 4–3 Ma) (An et al., 2001; Jacques et al., 2011; Molnar, Boos & Battisti, 2010; Song et al., 2010; Sun & Wang, 2005; Valdiya, 1999; Wan et al., 2007; Zhisheng et al., 2001). The development of the Asian monsoon system directly gave birth to the warm and humid climate in the south of China (Sun & Wang, 2005), which was maybe favorable for the geographical spread and speciation of amphibians (Che et al., 2010; Thorn & Raffaelli, 2001; Wu et al., 2013; Zhang et al., 2006). In addition, the climate oscillations that began about 2.8 million years ago, in the Late Pliocene (Deng et al., 2011), also provided the chance for diversification and speciation of many species (Zhang, Fengquan & Jianmin, 2000), such as birds (Lei, Qu & Song, 2014), the Tibetan woolly rhino (Coelodonta thibetana) (Deng et al., 2011) and stream-dwelling frog (Feirana quadranus) (Wang et al., 2012a). Molecular dating suggested that the TMRC of Clade A and Clade B was during the Oligocene (22.44 Ma–32.17 Ma) (Fig. 3, Table 1). At same time, ancestral area reconstructions supported Southeast Asia (W) as the ancestral area of Rhacophorus and the dispersal events happened from ancestral area of Clade A and Clade B (Fig. 2, node b). In addition, the land and sea in the Indo-Australian Archipelago changed greatly during this period (Lohman et al., 2011), which may promote the dispersal events from Southeast Asia. In Clade A, the time of the split of subgroups was estimated from 14.09 to 5.33 Ma (Table 1). In addition, the time of most nodes in Clade B also occurred during this period (Fig. 3). Based on the LTT plot analysis, there were two increased diversification rate periods (14–12 Ma, 10–4 Ma) in Rhacophorus (Fig. 4). The series of Tibetan Plateau rapid uplifts (from 15 to 2.66 Ma) dramatically changed the landscape, which resulted in the diversification of species or speciation in this area (Che et al., 2010; Deng et al., 2011; Gao et al., 2013; Jabbour & Renner, 2012; Lei, Qu & Song, 2014; Li et al., 2013; Shi et al., 1999; Tietze & Borthakur, 2012; Tietze et al., 2013; Tu et al., 2010; Wang et al., 2009; Zhang et al., 2006) and the biotic interchange between the Indian subcontinent and mainland Asia (Klaus et al., 2016). In addition, the intensified East Asian monsoon (∼15 Ma, ∼8 Ma and 4–3 Ma) directly gave birth to the warm and humid climate in the south of China, which was favorable for the geographical spread and speciation of amphibians (Che et al., 2010; Thorn & Raffaelli, 2001; Wu et al., 2013; Zhang et al., 2006). Obviously, the diversification events in Rhacophorus were in line with the time frame of the orogenic movement and climatic histories, especially the staged rapid uplift of the Tibetan Plateau and the enhanced Asian monsoon system (Figs. 2 and 3). Therefore, we think that the diversification and speciation events in Clade A and Clade B, are related to the staged uplift of the Tibetan Plateau and the subsequent chain-reaction events, such as the establishment of the Asian monsoon system, which facilitated the radiations and speciation of amphibians (Che et al., 2010; Thorn & Raffaelli, 2001; Wu et al., 2013; Zhang et al., 2006).

Overall, the evolutionary history of Rhacophorus originated approx 30 Ma Bp (Oligocene). Basically, it is the dispersal process from its ancestral area, Southeast Asia, toward India and East Asia. During the process, Rhacophorus diversified by multiple factors, such as geographic shifts, the staged rapid rises of the Tibetan Plateau with parallel climatic changes, the reinforcement of the Asian monsoons (15 Ma, 8 Ma and 4–3 Ma) and alternating glacial-interglacial oscillations.