Phylogeny of spiny frogs Nanorana (Anura: Dicroglossidae) supports a Tibetan origin of a Himalayan species group

Abstract Recent advances in the understanding of the evolution of the Asian continent challenge the long‐held belief of a faunal immigration into the Himalaya. Spiny frogs of the genus Nanorana are a characteristic faunal group of the Himalaya–Tibet orogen (HTO). We examine the phylogeny of these frogs to explore alternative biogeographic scenarios for their origin in the Greater Himalaya, namely, immigration, South Tibetan origin, strict vicariance. We sequenced 150 Nanorana samples from 62 localities for three mitochondrial (1,524 bp) and three nuclear markers (2,043 bp) and complemented the data with sequence data available from GenBank. We reconstructed a gene tree, phylogenetic networks, and ancestral areas. Based on the nuDNA, we also generated a time‐calibrated species tree. The results revealed two major clades (Nanorana and Quasipaa), which originated in the Lower Miocene from eastern China and subsequently spread into the HTO (Nanorana). Five well‐supported subclades are found within Nanorana: from the East, Central, and Northwest Himalaya, the Tibetan Plateau, and the southeastern Plateau margin. The latter subclade represents the most basal group (subgenus Chaparana), the Plateau group (Nanorana) represents the sister clade to all species of the Greater Himalaya (Paa). We found no evidence for an east–west range expansion of Paa along the Himalaya, nor clear support for a strict vicariance model. Diversification in each of the three Himalayan subclades has probably occurred in distinct areas. Specimens from the NW Himalaya are placed basally relative to the highly diverse Central Himalayan group, while the lineage from the Tibetan Plateau is placed within a more terminal clade. Our data indicate a Tibetan origin of Himalayan Nanorana and support a previous hypothesis, which implies that a significant part of the Himalayan biodiversity results from primary diversification of the species groups in South Tibet before this part of the HTO was uplifted to its recent heights.

Consequently, our understanding of the origin and historic biogeography of the terrestrial faunas inhabiting the HTO is far from conclusive, and has been hindered by a lack of and potential misinterpretation of data (Favre et al., 2015;Renner, 2016;Spicer, 2017;Su et al., 2019).
The key to disentangling the paleoclimatology and biogeography of the HTO lies particularly in the Himalaya (Molnar, 1986;Searle, 2013) for the following reasons: i) Most of the modern geological models differ with respect to the uplift history of this part of the mountain system (Schmidt, Opgenoorth, & Miehe, 2016, and refs. therein), ii) the uplift of the Greater Himalaya markedly influences the climate in the interior of High Asia and on a global scale (Boos & Kuang, 2010;Hodges, Hurtado, & Whipple, 2001;Sanwal et al., 2013), and iii) South Tibet may have been an important evolutionary center during the Cenozoic, which impacted the modern faunas of Central and East Asia (Schmidt, Opgenoorth, Holl, & Bastrop, 2012;Weigold, 2005).
So far, at least three different scenarios exist in the literature on the origin of the Himalayan wildlife: first, the long-held belief of a faunal immigration scenario. Under this scenario, Himalayan taxa are assumed to have migrated along the Himalaya, that means, across preexisting barriers (deep valleys and high mountain ridges) (Martens, 2015). Many phylogenetic studies with a primary focus on the Greater Himalaya have been conducted in groups with higher dispersal potential, such as birds, butterflies, and plants with windor bird-dispersed seeds (e.g., Deodati, Cesaroni, & Sbordoni, 2009;Favre et al., 2015;Mani, 1986;Martens, Tietze, & Päckert, 2011;Rajbhandary, Hughes, Phutthai, Thomas, & Shrestha, 2011;Xie, Ash, Linde, Cunningham, & Nicotra, 2014;Zhang, Kang, Zhong, & Sanderson, 2012). For the majority of these organisms, it has been reported that they originated through long-distance dispersal from the mountains of China-Indochina along the southern slope of the Himalayan chain, associated with very little in situ speciation (Deodati et al., 2009;Johansson et al., 2007;Liu et al., 2013;Martens & Eck, 1995;Martens et al., 2011;Rajbhandary et al., 2011;Tabata, 2004;Xie et al., 2014;Zhang, Meng, Zhang, Vyacheslav, & Sanderson, 2014). For several Palearctic species groups, on the other hand, it has been shown that the influx of organisms came from the West along a climatically temperate corridor that enabled dispersal from Central Asia and the Pamiro-Alai region into the Himalaya (Alcaide, Scordato, Price, & Irwin, 2014; Martens, 2015).
Second, recent studies of Asian forest-dwelling beetles and anurans suggest a Tibetan-origin scenario for at least some Himalayan faunal components (Hofmann et al., 2017;Schmidt et al., 2012). A similar explanation has been proposed for the evolution of highly isolated non-wind-dispersed Himalayan alpine plants (Miehe et al., 2015). This Tibetan-origin scenario postulates that during the early phase of mountain uplift, South Tibet was an independent center of evolution that was well separated from other mountainous regions (Schmidt, 2006;Schmidt et al., 2012). According to this model, the final Plateau uplift, associated with the desiccation of the Tibetan highlands, "squeezed out" ancestral lineages that subsequently diversified by vicariance (Schmidt et al., 2012;Yin & Harrison, 2000).
If so, adaptation to high altitudes and primary diversification of local species groups would have happened or at least been initiated in the high mountains of Paleo-Tibet, potentially long before the final uplift of the Greater Himalaya. Colonization of the Greater Himalaya would have taken place in the course of its growth by ancestral species, originating in the immediately adjacent mountains to its north (i.e., at the southern edge of what is now Tibet). This area of origin, however, had been lost due to aridification, leading to faunal extinction or turnover there (Schmidt et al., 2012).
The third scenario is based on a previous study in Nanorana spiny frogs (Che et al., 2010). It implies a setting of explicit vicariance, driven by geographical isolation, climatic conditions, and an assumed low dispersal ability of spiny frogs. Accordingly, species formation among major lineages occurred as the species were "trapped" in the mountain mass and became separated when it uplifted (Che et al., 2010). However, in that study only three samples from the Himalaya were included, but none from the vast areas of the central and western parts of this mountain range.
Spiny frogs of the genus Nanorana, subfamily Dicroglossinae, are a characteristic species group of the HTO, living mostly in swift boulder-strewn streams (Dubois, 1975). These frogs are found across the Himalayan arc from northern Pakistan and northern India, through Nepal, Sikkim, and Bhutan, to western China (Hengduan Shan), Myanmar, Thailand, Laos, northern Vietnam, and to montane southern China (Frost, 2019). The genus comprises 30 accepted species (Frost, 2019; Appendix 1) and is subdivided into three subgenera (Nanorana Günther, 1896, Paa Dubois, 1975, and Chaparana Bourret 1939Frost, 2019). However, the phylogenetic and taxonomic relationships among spiny frogs, especially within the subgenera, but also within the subfamily, are still not completely resolved, as shown by taxonomic refinement during the last decade including the description of new species (Che et al., 2009;Frost, 2019;Huang et al., 2016;Jiang et al., 2005;Pyron & Wiens, 2011).
Herein, we explore whether the phylogeny of Himalayan Nanorana is better explained by an east-to-west immigration into the Himalaya, by a Tibetan-origin scenario, or by a strict vicariance model. We generated a gene genealogy using mitochondrial and nuclear DNA fragments, reconstructed phylogenetic networks, and assessed the ancestral areas and the relative role of dispersal, vicariance, and extinction to examine the history of Himalayan spiny frogs. Based on the nuDNA, we also generated a time-calibrated species tree.
In the case of an immigration scenario, the gene tree is expected to show nested phylogeographic distribution. The dating of the divergence between the lineages from the western part of the Greater Himalaya to a time younger than those at the eastern part would likewise argue for an immigration scenario, as would the occurrence of close relatives of Himalayan taxa in adjacent mountains. Gene flow among the regions at the southern slope of the Himalayan mountain chain from east to west or vice versa would result in a directional or intermingled nature of the haplotype patterns, due to (repeated) immigration into the Himalaya.
Under the Tibetan-origin scenario, endemic Himalayan lineages are expected to show disjunct distribution with no close relationships to lineages occurring in adjacent parts of the mountain system, and which are of greater phylogenetic age. Himalayan lineages might be the same age as lineages from the eastern HTO (or might be even older than these), while lineages from the Tibetan Plateau should be closely related to any of the lineages from the HTO margin. Findings of shared haplotypes in adjacent drainage or mountain systems would provide evidence of dispersal in more recent times.
Under the strict vicariance scenario, dispersal of the lineages into the different parts of the HTO should have occurred more or less simultaneously during an early phase of the uplift. Lineages from the Tibetan Plateau are predicted to be at least as old as those occurring in the different parts of the Greater Himalaya. Consequently, similar to the Tibetan-origin scenario, lineages are expected to be represented by deep phylogenetic branches and to correspond strongly to certain parts of the HTO. However, geographical ranges of sister lineages are not expected to overlap. We would expect no evidence for contemporary dispersal paralleling the Greater Himalaya, because barriers to dispersal in this part of the HTO were never as effective as today.
Our results may help to better understand the evolution and Cenozoic history of Himalayan biodiversity. They may also provide insights into how molecular phylogenies of poorly dispersing extant species groups can be integrated into advanced reconstructions of the highly complex geomorphological and paleoecological evolution of the HTO.

| Sampling and DNA extraction
A total of 150 individual samples of Nanorana spiny frogs covering major parts of the southern slope of the Himalaya chain were included in this study (Figure 1), coming from scientific collections (Chinese Academy of Science, CAS; Natural Museum of Erfurt, NME; Russian Academy of Science, RAS; Appendix 2). Samples were supplemented by NCBI sequence data from further Nanorana and outgroup species (Figure 1 and Appendix 2).

| Sequence alignment
Ribosomal RNA (rRNA) 12S and 16S sequences were aligned based on their secondary structures using RNAsalsa 0.8.1 (Stocsits, Letsch, Hertel, Misof, & Stadler, 2009) and the ribosomal structure model of Bos taurus downloaded from www.zfmk.de/en/ resea rch/resea rch-centr es-and-group s/rnasalsa. We generated a first structure-based multiple sequence alignment using complete rRNA sequences of several Nanorana species (N. yunnanensis, N. taihangnica, N. parkeri, N. pleskei) and further dicroglossid species (Fejervarya cancrivora, Hoplobatrachus rugulosus, Limnonectes fragilis, Quasipaa shini) available from NCBI (Appendix 2). Remaining sequences were then aligned with this initial file using the MUSCLE algorithm (Edgar, 2004) in MEGA X (Kumar, Stecher, Li, Knyaz, & Tamura, 2018). To exclude ambiguously aligned sites, we used trimAL (Capella-Gutierrez, Silla-Martinez, & Gabaldon, 2009) under the automated1 option, which implements a heuristic selection to decide the most appropriate mode depending on the alignment characteristics. The sequences of the protein-coding genes were also aligned with MUSCLE using default settings in MEGA X.
Alignment on the basis of nucleotides and amino acids produced similar results, since no ambiguities, such as deletions, insertions, or stop codons, were found.
We did not phase our nuclear data as many outgroup and ingroup taxa had only single representative individuals. The data set was partitioned a priori by gene fragments, and PartitionFinder 1.1.1 (Lanfear, Calcott, Ho, & Guindon, 2012) was applied to optimize partitions using linked branch lengths, the Bayesian information criterion (BIC), the "greedy" search algorithm, and the BEAST or RAxML option (Lanfear, Calcott, Kainer, Mayer, & Stamatakis, 2014). We ran RAxML with the GTRGAMMA model and 1,000 bootstrap replicates. BEAST2 analyses were based upon five independent runs with a chain length of 250 million each, thinning interval of 25,000, a lognormal relaxed clock model, a Yule tree prior, a random tree as starting tree, and the site models selected using bModelTest version 1.1.2 (Bouckaert & Drummond, 2017).

| Time-calibrated species tree estimation
A time-calibrated species tree was estimated using the uncorrelated relaxed clock method in *BEAST. Because mtDNA sequence data are expected to be systematically biased toward the calibration point and, thus, to overestimate divergence times (Zheng, Peng, Kuro-o, & Zeng, 2011), we used only the unphased nuclear gene data for our dating analysis. This data set contained 182 taxa (of these, 32 were obtained from GenBank). We used the same run parameters as for the gene tree (i.e., five independent runs, 250 Mio generations, 25,000 sampling frequency, Yule speciation process, random starting tree, bModelTest on gene partitions), combined the runs after checking for convergence of modeled parameters, and annotated a maximum credibility clade tree using a burn-in of the first 10% of the sampled trees.
Three carefully chosen fossil calibrations were selected to obtain divergence dates on the species tree and applied using an offset lognormal distribution: i) a minimum age of 33.9 Ma for the most recent common ancestor (MRCA) of Ranoidea based on the fossil Thaumastosaurus gezei (Rage & Roček, 2007) Roelants, Haas, & Bossuyt, 2011). Since anuran fossil records are notoriously rare and do not exist for dicroglossid frogs, we therefore relied on using previously published fossil-calibrated divergence estimates to place a broad prior on the divergence of Nanorana and Quasipaa as further, internal calibration point. Based on Bossuyt et al. (Bossuyt, Brown, Hillis, Cannatella, & Milinkovitch, 2006), we used a calibration of 38.1 ± 9.4 Ma (28.7-47.5 Ma) for that node with a normal distribution prior. In all runs, we constrained Dicroglossinae to be monophyletic, as well as Microhylidae, Nanorana, and Natatanura, respectively; we also constrained Sooglossus to be the sister clade to Ranoidea since these relationships have been well established (Roelants et al., 2011;Wiens, Sukumaran, Pyron, & Brown, 2009;Zhang et al., 2013).

| Haplotype network
We used SplitsTree4 v.14.8 (Huson & Bryant, 2006) and the NeighborNet algorithm based on the uncorrected p-distances to construct a haplotype network of the concatenated mtDNA and (unphased) nuDNA sequence data set, respectively. To increase resolution, we reduced the two data sets to sequences that had less than 10% missing data (and less than 20% when analyzing mtDNA and nuDNA together).

| Ancestral area reconstruction
To infer ancestral distributions and to elucidate the biogeographic history of Himalayan spiny frogs, we used two different approaches, namely Statistical Dispersal-Extinction Cladogenesis (S-DEC) and Statistical Dispersal-Vicariance Analysis (S-DIVA) implemented in RASP v.4.1beta (Yu, Harris, Blair, & He, 2015). Assignment of the samples to biogeographic units within the HTO proved to be difficult because current models of the spatiotemporal evolution of the HTO and of its biogeographic history are partly very different and no widely accepted model exists for the region. These differences arise particularly from biogeographic conclusions depending on whether they are based on distributional patterns of species groups with high dispersal ability (e.g., wind-dispersed plants, butterflies, birds) or with very low dispersal ability (e.g., ground beetles; for an overview, see Martens, 2015). Because anurans disperse actively "on foot," they are considered poor dispersers. They show, moreover, remarkable stasis in ecological niches, suggesting that dispersal will have been historically constrained between similar climatic conditions (Wiens, 2011). Due to the potentially restricted nature of their movements and dependence on relatively stable local environments, anuran species distributions should reflect a high level of paleoenvironmental history of a region. For the HTO, however, this history is still largely unknown (Renner, 2016). Therefore, we provisionally assigned the samples of dicroglossid species to certain parts of the HTO based on the species' contemporary distri-  Figure 1a), locality numbers refer to samples/ sequences listed in Appendix 2 (Clark et al., 2005;Schmidt et al., 2015). Species with wider distribution were assigned to the unit where most of the distributional area of the respective species is located. For example, N. liebigi has occasionally been reported from the westernmost part of the East Himalaya but is widely distributed in the Central Himalaya and is therefore assigned to the latter. Nanorana quadranus also occurs on the northeastern foothills of the Tibetan Plateau but is widely distributed along the Sichuan Basin (see Appendix 1) and consequently assigned to the latter. The northernmost records of N. yunannensis are located in Sichuan along the eastern Plateau margin; however, the distribution of this species is centered in the southeastern part of the HTO (unit F). Further, distributional areas of some species, for example, Quasipaa spinosa and Q. verrucospinoa (unit G), appear to overlap with species assigned to unit F, though these G-assigned species occur at significantly lower elevations than those of the montane regions of the Plateau margin.

| Phylogenetic analyses
The ML and BI gene trees (mtDNA + nuDNA) were well resolved and yielded almost identical tree topologies, except for the place-  (Kass & Raftery, 1995). Optimal partitions and substitution models were assessed in PartitionFinder 1.1.1 (Lanfear et al., 2012) with branch lengths linked, a "greedy" search algorithm, the BIC, and the BEAST option (we run topology tests based on both, best scheme of gene partition and of genes + codon partition for protein-coding genes). To reduce computation time, we included only data of dicroglossid frogs and excluded outgroups which had been used for the calibration analysis. The BF model selection clearly preferred model m1 based on both the stepping stone and path sampling method (Appendix 5). In fact, our split trees (Appendix 4) also showed a NW-Central Himalaya sister group, especially based on nuclear genes.
The geographical distribution of the Himalayan lineages, which are strongly supported in the combined mtDNA and nuDNA data analyses, is presented in Appendix 6a-d. These results show two apparently contradictory aspects. On the one hand, strong support is evident for local endemism of lineages to different valley systems within the Greater Himalayan mountain arc, even at intraspecific level. On the other hand, contemporary dispersal within, but also across, drainage boundaries is shown by some of the spe-

| Divergence times in spiny frogs and biogeographic reconstructions
Topology based on nuDNA was consistent with that based on the concatenated mtDNA and nuDNA sequence data (Figures 2 and   3). Dating analysis suggests an origin of the spiny frogs (Quasipaa, Nanorana) in the Lower Miocene (23 Ma, 95% HPD 11-35 Ma; F I G U R E 2 Maximum-likelihood (left) and Bayesian inference tree (right) based on concatenated mtDNA and nuDNA sequence data. Numbers at branch nodes refer to bootstrap values > 70% (ML tree, 1,000 replicates) and or posterior probabilities ≥ 0.9 (Bayesian inference tree). The subgenus (Chaparana, Nanorana, or Paa) is also indicated at the respective clade. A map below the trees shows the localities of the samples within the phylogenetic clades; the color code corresponds to the color of the respective clade Figure 3 and Appendix 7). This date is consistent with a transcriptome-based study (Sun et al., 2018) and only slightly younger

| Origin and evolutionary history of Himalayan spiny frogs
The commonly believed immigration scenario for faunal elements in the Himalaya assumes they originated via long-distance dispersal, for example, from the mountains of China-Indochina along the Himalayan chain (see Hofmann et al., 2017;Schmidt et al., 2012 and references therein). Accordingly, as observed in birds, butterflies, and plants (e.g., Li et al., 2019;Manish & Pandit, 2018;Martens et al., 2011;Xie et al., 2014), all lineages of Himalayan spiny frogs recorded at present would have dispersed from east to west or would have appeared during the potential range expansion along the Greater Himalaya. Thus, they must have been able to cross the many north-south stretching mountain ranges and the epigenetic transverse valleys, by which the Himalaya has always been intersected and that might form effective barriers to dispersal of amphibians (Sánchez-Montes, Wang, Ariño, & Martínez-Solano, 2018;Zhou et al., 2017). Relating thereto, the phylogenetic placement of the clade from the NW Himalaya is of particular interest. Based on the immigration model, it can be expected that the NW Himalaya clade represents a terminal lineage of the highly diverse Central Himalaya clade. However, in all our analyses, the NW Himalaya clade was placed basally relative to the Central Himalaya clade (model m1) or as sister clade to the East Himalaya clade (model m2), though m1 had a significant higher likelihood than m2 based on the topology tests (Appendix 5). This not only indicates the presence of ancestral lineages in the NW Himalaya but also provides a strong argument against dispersal of spiny frogs paralleling the Greater Himalayan mountain arc.
The results of our analyses demonstrate that the ancestral Asian spiny frogs descended, during the Lower Miocene, from a tropical ancestor in eastern China that gave rise to the genera Quasipaa and Nanorana, with Nanorana comprising three major clades (subgenera), namely the lower montane Chaparana, the (high-)montane Paa, and the (sub)alpine, nominal subgenus Nanorana. Up to this point, our results confirm previous findings (Che et al., 2010). However, we recovered a highly supported distinct Central Himalayan Paa clade, which rapidly diversified during the Pliocene and Pleistocene. It includes the high-montane species Nanorana liebigii, N. cf. blanfordii, N. cf. ercepeae, N. cf. polunini, N. cf. rarica, N. cf. rostandi, and several undescribed lineages from the southern slopes of the Himalaya. This clade shows strong phylogeographic structure and constitutes-most parsimoniously together with the NW Himalaya clade-the sister clade to high-montane Paa species occurring in the East Himalaya and F I G U R E 3 Time-calibrated species tree of Asian spiny frogs using *BEAST and ancestral areas reconstructed with S-DEC and S-DIVA (above and below branches). Outgroups were excluded for readability. Bayesian posterior probability values ≥ 95% are indicated by an asterisk at the respective node. Codes in square brackets next to undescribed lineages (N. sp.) specify their region of origin (HimPr = Himachal Pradesh; Chainpur = Chainpur Himal) or simply a working label (A, B, C). Nanorana subgenera are indicated to the right of the clades. Pie charts at nodes show probabilities of alternative ancestral range for MRCAs obtained from S-DEC and S-DIVA, corresponding to the color key legend. Colored boxes around clades indicate their areas of distribution: A = NW Himalaya, B = Central Himalaya, C = East Himalaya; D = Transhimalaya and adjacent parts of the Tibetan Plateau; E = (sub)alpine parts of the eastern margin of the Tibetan Plateau; F = high-montane regions of the southeastern margin of the Tibetan Plateau including the high mountains of northern Vietnam; G = subtropical and meridional eastern China; H = Sichuan Basin and mountains of Northeast China the adjacent high-montane regions of the southeastern margin of the Tibetan Plateau. According to our dating results, the MRCA of these sister clades (Paa) and the Plateau clade (Nanorana) occurred in the Upper Miocene, ca. 9 Ma. The MRCA of Chaparana and Nanorana + Paa occurred almost ten million years earlier, during the Lower Miocene (ca. 18 Ma). A similar age for that split has been found in other studies (Che et al., 2010;Chen et al., 2017;Wiens et al., 2009). Between these two nodes, the ancestral lineages of Nanorana + Paa apparently did not diversify for several million years before major diversification occurred. There are two alternative explanations for such a pattern.
Either these ancestral lineages persisted in their area of origin, which had remained stable over a long period until it facilitated new habitats (or niches), or this ancestral area successively disappeared in the course of the stepwise uplift and drying of Tibet led to extensive species extinction and range shifts into new suitable areas. The first scenario is unlikely because mountain building is accompanied by both gradual and abrupt environmental changes, which ultimately drive biodiversity dynamics (Huang, Meijers, Exres, Mulch, & Fritz, 2019).  (Wang, Shi, & Zhou, 1982;Wang et al., 2008). Due to the rising Himalayan mountain belt and the associated continuous aridization of southern Tibet, many ancestral lineages might have successively been lost to extinction or been forced to track the displaced suitable environment along the transverse valleys of the Himalaya, such as the Arun, Brahmaputra, Kali Gandaki, or the Indus catchment (see Figure 1b and Appendix 1 for geographical reference), which otherwise formed an effective barrier against the north-south dispersal of these species (Schmidt et al., 2012). Shifting ranges in response to substantial climatic changes is a common phenomenon that defines phylogeography especially in temperate taxa and mountain biota (Giezendanner, Bertuzzo, Pasetto, Guisan, & Rinaldo, 2019;Hoorn, Perrigo, & Antonelli, 2018). The strong re-  (Schmidt et al., 2012). The Kashmir Himalaya may be considered an "exile area," meaning it is the distributional area of the Paa descendants which originated from the western slope of Paleo-South Tibet.
Similarly, from the eastern macroslope of the HTO, ancestral species followed the main river gorges into the southeastern Plateau margin (see modern distribution of N. chayuensis, N. maculosa in Appendix 1 and Figure 1a; for evolution of the great river systems see Brookfield, 1998;Lang & Huntington, 2014).
Further support for a South Tibetan-origin scenario is provided by the phylogeographic pattern of the Tibetan Plateau frog (N. parkeri), which is the only known Nanorana taxon that occurs in the alpine zone while all other members of this genus are high-montane to subalpine taxa. This points to a stepwise high altitude adaptation. The relatively recent divergence time between the alpine N. parkeri and the subalpine N. pleskei/N. ventripunctata from the eastern Plateau margin (< 5 Mya) fits the concept of a final uplift of the Tibetan Plateau rather late in the Neogene (Su et al., 2019). An almost identical biogeographic history has been demonstrated in Scutiger boulongeri, the Tibetan alpine toad (Hofmann et al., 2017).
In contrast to the Tibetan-origin scenario, the strict vicariance model suggests that spiny frogs colonized the uplifting HTO and were then "trapped" in their habitats when the mountains evolved (Che et al., 2010). Considering this scenario, we must assume that, for example, the spiny frog lineage that occurs on the central Tibetan Plateau (N. parkeri) reflects its ancestral spatial distribution and that the uplifting Tibetan Plateau triggered further diversification within that lineage. Since modern geological models consistently suggest that at least the southern central part of the Tibetan Plateau had reached significant heights before the Greater Himalaya was uplifted (Deng & Ding, 2015;Mulch & Chamberlain, 2006;Wang et al., 2008), the Plateau lineage is expected to branch basally relative to the other clades occurring in the Greater Himalaya. However, N. parkeri clusters together with species from the eastern Plateau margin, representing a relatively young, terminal branch as discussed above (Figure 2). For the current distribution of spiny frogs in the Greater Himalaya, two alternative explanations emerge under the vicariance model. On the one hand, the ancestral Paa lineages may have colonized the different parts of the Greater Himalaya before this mountain arc was uplifted.
Consequently, these lineages must have been elements of a tropical fauna, while the Nanorana lineage of the Tibetan Plateau must have been already adapted to high altitudes. This implies that the ancestral species were able to adapt to new conditions under a dramatically changing environment in their current locations. Given the existence of subtropical vegetation (cloud forests) in the Tibetan Himalaya north of the Central Himalaya during the Miocene and even during the Late Pliocene, as demonstrated by plant fossil data (Su et al., 2019;Xu, 1981Xu, , 1982, this environment must have been rapidly transformed from subtropical montane to arid alpine habitats. However, since most amphibian populations can only accommodate very gradual environmental changes, their phenotypic plasticity and genetic variance might be insufficient to generate the phenotypic changes necessary to cope with such massively transformed conditions (Urban, Richardson, & Freidenfelds, 2014).
On the other hand, like the immigration scenario, the ancestral Paa lineages may have dispersed parallel to the growing Greater Himalayan mountain arc-a concept that is not supported by our data as we have discussed above.
A central argument of the strict vicariance scenario is the geographical and ecological separation of Nanorana lineages induced by the uplift of the Himalayan and Tibetan region. This argument, however, is compatible neither with a distribution that stretches over several mountains and deep valleys (Appendix 6a-d), nor with the overlapping altitudinal ranges (Hu, Xie, Li, & Jiang, 2011) and sympatry of some spiny frog species, e.g., Chaparana yunnanensis, C. unculuanus, and Paa maculosa; Appendix 1). Given that barriers to dispersal were never as high as in the present HTO, the modern distributional patterns raise questions about a strict vicariance scenario.

| Conclusions
Phylogeography of the HTO encompasses approximately 50 million years of evolution and one of the most complex mountain systems on the planet. Given the resulting biodiversity, generalizations are difficult, but some patterns, particularly those that are found congruent among different taxonomic groups, are especially informative. Our results show strikingly similar phylogeographic patterns as observed in other Himalayan faunal elements that are assumed to have originated in South Tibet, for example, forest-dwelling Pterostichus ground beetles (Schmidt et al., 2012) and Scutiger lazy toads (Hofmann et al., 2017). These patterns are indicated by i) the presence of several distinct lineages of a monophyletic clade each occurring in restricted areas along the Greater Himalaya, ii) the presence of at least one ancient lineage endemic to the central portion of the Himalaya, and iii) lineages which are geographically adjacent are not necessarily closely related. Indeed, genetic separation, diversity, and geographical distribution of Himalayan spiny frogs seem to be best explained by a South Tibetan origin rather than by the alternative immigration or strict vicariance scenario. A Tibetan origin of high-montane faunal elements of the Himalaya is coherent with isotope records as well as with fossil evidence of Late Paleogene and Miocene subtropical to temperate vegetation north of the Greater Himalaya (Su et al., 2019;Wang, Deng, & Biasatti, 2006;Xu et al., 2012), suggesting the existence of cloud forests in Paleo-Tibet.
Although our molecular data markedly support a Tibetan origin for Himalayan spiny frogs, some ambiguities remain, for example, in the phylogenetic placement of some Himalayan lineages. Denser sampling in additional mountain systems of the Greater Himalaya and the use of a more extensive set of nuclear markers, for example, through the use of next-or third-generation sequencing, will be required to improve the resolution for relationships of Himalayan spiny frogs and allow the tracking of the true pattern of their cladogenesis. Potentially, this would also reveal numerous additional, undescribed lineages.
We hope that this work will encourage further studies aiming to analyze the biogeographic history in a range of additional taxonomic groups in the Himalaya, to assess the evolutionary significance of the HTO and to understand modern distributional patterns against the background of existing geological scenarios.

ACK N OWLED G EM ENTS
We are very grateful to Matthias Hartmann and Leo J. Borkin for providing tissue samples, Catarina Ferreira for comments on the manuscript, and Santos Tamang and Markus Sandner for assistance during fieldwork. We would also like to express our gratitude to the two anonymous reviewers,

DATA AVA I L A B I L I T Y S TAT E M E N T
All sequences were uploaded to GenBank; details on individual samples and accession numbers are available in Appendix 2.