Shared response to changes in drainage basin: Phylogeography of the Yunnan small narrow‐mouthed frog, Glyphoglossus yunnanensis (Anura: Microhylidae)

Abstract Aim With the late Cenozoic uplift of the Qinghai–Tibetan Plateau (QTP), drainage of the southeastern edge of the QTP changed significantly. However, the impact of this dramatic change on the geographical distribution and genetic diversity of endemic organisms is still poorly understood. Here, we examined the geographical patterns of genetic variation in the Yunnan small narrow‐mouthed frog, Glyphoglossus yunnanensis (Microhylidae), and two alternative hypotheses were tested: That is, the geographical distribution of genetic variation was determined by either the contemporary drainage basin or historical drainage basins. Location The Mountains of southwest China. Materials and methods Analyses were based on 417 specimens collected from across the distribution of the species. We reconstructed the genealogy (Bayesian and maximum parsimony methods) and assessed demographic history based on DNA sequencing data from mitochondrial and nuclear markers. We also mapped the genetic diversity and estimated the divergence times by a relaxed clock model. Results The species has maintained a relatively stable population size without recent population expansion. Four major maternal lineages were identified with good support, one representing a possible cryptic species and the other three showing further subdivision. The distribution of these deeply differentiated lineages/sublineages corresponded well to geographical regions. The secondary contact zones and phylogeographic breaks in distinct lineages of G. yunnanensis were almost concordant with those of Nanorana yunnanensis. Main conclusions Lineage division conformed to the hypothesis of drainage system evolution, that is, the phylogeographic pattern of G. yunnanensis was shaped by historical drainage patterns. Concordance in phylogeographic patterns may suggest a shared response to common hydrogeological history and also might indicate that there was more contribution of the drainage history than ecological or life‐history traits in structuring genetic variation between these two disparate codistributed taxa G. yunnanensis and N. yunnanensis.


| INTRODUC TI ON
The southeastern edge of the Qinghai-Tibetan Plateau (QTP) is an orogenetically young region that encompasses two geologically dynamic and biodiverse areas, namely the Hengduan Mountain and Indo-Burma hot spots (Myers, 2003;Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000;Tang, Wang, Zheng, & Fang, 2006). These mountainous regions of southwest China and the rearrangement of their major river drainage systems are considered evolutionarily important (Stewart, Lister, Barnes, & Dalen, 2010). The Yunnan-Guizhou Plateau (YGP) formed a planation surface (or peneplain) with smooth relief during the middle-late Pliocene (Cheng, Liu, Gao, Tang, & Yue, 2001;He & He, 1993;He, He, & Zhu, 1985). With the strong uplift of the QTP, the Yunnan planation surface began to collapse from the late Pliocene to early Pleistocene (He et al., 1985). Along with the obvious change in topography, the previous lake-centered drainage network was broken. In the late early Pleistocene, the most significant drainage system changes occurred on the Central Yunnan Plateau (CYP) and in northwest Yunnan (Cheng et al., 2001;Clark et al., 2004).
Increasing evidence has shown that the rearrangement of drainage systems has played a key role in shaping current geographical patterns of genetic and species diversity of aquatic organisms in many parts of the world (Adamson, Hurwood, & Mather, 2012;Burridge, Craw, & Waters, 2006;de Bruyn et al., 2013;Hurwood & Hughes, 1998;Kozak, Blaine, & Larson, 2006;Thomaz, Malabarba, Bonatto, & Knowles, 2015). The historical drainage systems in the southwest mountains of China are markedly different from present systems (Brookfield, 1998;Clark et al., 2004). Prior to the most recent uplift of the QTP, the initial drainage system along the plateau margin primarily consisted of tributaries of a single, southward flowing paleo-Red River. With uplift of the QTP, disruption of the paleo-drainage system occurred by river capture and reversal. The Jinsha, Yalong, and Dadu rivers connected to form the modern Jinsha River, which redirects drainage away from the southward Red River eastward into the East China Sea, and the Upper Mekong and Upper Salween rivers drained into their modern drainage positions. The impact of paleo-drainage rearrangement on species diversification and geographical distribution has been evaluated via molecular phylogenetic studies among freshwater fish (Guo, He, & Zhang, 2005;He & Chen, 2006;Peng, Ho, Zhang, & He, 2006;Rüber, Britz, Kullander, & Zardoya, 2004) and a few riparian plant species and amphibians (Wang, Mao, Zhao, & Wang, 2013;Yan et al., 2013;Zhang, Chen, et al., 2010;Zhang, Comes, & Sun, 2011).
The potential for diversification triggered by environmental changes or climatic oscillations likely varies among taxa. Organisms with different ecological preferences may have differing responses to the same event, as well as contrasting evolutionary histories.
Previous research on freshwater fish has provided evidence that changes in ancient drainage systems contributed to allopatric speciation and population range expansion (Guo et al., 2005;He & Chen, 2006;Peng et al., 2006;Rüber et al., 2004). Furthermore, phylogeographic studies on riparian plant species endemic to the hot-dry river valleys of the eastern Sino-Himalayan region have indicated that current genetic structures were historically sculpted by paleo-drainage patterns (Yue, Chen, Sun, & Sun, 2012;. In contrast, however, recent research on pine species endemic to the major river valleys of southwest China has indicated that spatial genetic structure is a reflection of current geography and environmental factors rather than hydrogeological history (Wang, Mao, Zhao, & Wang, 2013).
The two patterns mentioned above have also been observed in amphibians from southwest China. For example, the current population structure of the Yunnan spiny frog, Nanorana yunnanensis, was primarily shaped by the historical drainage system (Zhang, Chen, et al., 2010). Studies on other amphibians from the same region have revealed that present-day genetic structures were shaped predominantly by Pleistocene climatic oscillations (Li, Yu, Rao, & Yang, 2012;Yu, Zhang, Rao, & Yang, 2013), with deep river valleys identified as strong geographical barriers to dispersal (Li, Chen, Tu, & Fu, 2009;Yuan et al., 2016;Zhang, Rao, Rao, Yang, Yu, & Wilkinson, 2010). Differentiation in ecological adaptations among taxa has also been proposed to explain variation in the geographical patterns of genetic diversity. For example, the genetic structures of species living in lentic environments, for example, Yunnan pond frog (Babina pleuraden) (Li et al., 2012) and Red Knobby Newt (Tylototriton shanjing) (Yu et al., 2013), differ from those inhabiting lotic environments, for example, N. yunnanensis (Zhang, Chen, et al., 2010). Drainage system evolution may have significant influence on the genetic structure of species preferring lotic water (Yu et al., 2013). However, Leptobrachium ailaonicum (Zhang, Rao, et al., 2010) occupies montane streams similar to N.
yunnanensis (Zhang, Chen, et al., 2010), but both have highly distinct phylogeographic patterns. Thus, the responses of organisms to changes in drainage patterns are complex. Further research is needed to explore the impact of hydrogeological history of drainage basins on the geographical distribution and genetic structure of endemic taxa in southwest China. traits in structuring genetic variation between these two disparate codistributed taxa G. yunnanensis and N. yunnanensis.

K E Y W O R D S
drainage history, genetic structure, Glyphoglossus yunnanensis, phylogeography, secondary contact, shared response, southwest China Glyphoglossus yunnanensis occurs in southwest China (Yunnan, southern Sichuan, and western Guizhou), adjacent northern Vietnam, and presumably in adjacent Laos (Frost, 2019). It inhabits moist and soft soil, and breeds in puddles of rainwater and ditches at elevations of 1,700-3,600 m (personal observation). Its distribution overlaps considerably with N. yunnanensis and B. pleuraden. Vicariance theory assumes that common genealogical distribution patterns and geographical position of lineage breaks across multiple organisms in the same geographical range may originate from shared biogeographical history (Rissler & Smith, 2010). Thus, G. yunnanensis is an ideal species for exploring the mechanisms driving the phylogeographic patterns of species in this region from a comparative phylogeographic perspective.
Our objective was to examine the impacts of drainage system rearrangement on the genetic patterns of G. yunnanensis. We tested two alternative hypotheses: That is, geographical distribution of the genetic variation of G. yunnanensis was determined primarily by the contemporary drainage basin or by the historical drainage basins. We collected samples from throughout the distribution range of G. yunnanensis. Using DNA sequence data from both mitochondrial and nuclear genes, we analyzed genealogical and population genetics of the nucleotide sequence data. If the geographical structure of the genetic variation was best explained by the contemporary drainage basin, we would expect haplotypes from the same drainage basins have closer phylogeographic relationships and cluster together. On the other hand, phylogeographic breaks would be expected among haplotypes sampled from the same drainage basins if historical drainage basins were responsible for the phylogeographic structure.

| Sampling
From 2013 to 2016, a total of 417 individuals were collected from 48 localities throughout the distribution range of G. yunnanensis, except Vietnam (Supporting Information Table S1 in Appendix S1;  Table S1 in Appendix S1. Solid triangles, circles, diamonds, and fivepointed star correspond to major matrilines (lineages W, C, E, and LD, respectively) in Figure 2. Overlap of different symbols denotes sympatric occurrence of distinct clades. Grayish dotted line designates geographical distributions of four main haplotype subclades of clade C in Figure 2. Both (b) and (c) are adapted from T. C. . Paleo-drainage pattern prior to major river reversal/capture (b) and modern river pattern after putative capture and reversal events (c) are shown, respectively

| Laboratory procedures
Using the proteinase K/sodium dodecyl sulfate (SDS) method (Sambrook & Russell, 2001), genomic DNA was extracted from tissue samples. Partial segments of two mitochondrial genes, cytochrome oxidase I (COI) and cytochrome b (CYTB), were amplified for all individuals. Mitochondrial ribosomal subunit 16S rRNA and nuclear gene tyrosinase (Tyr) were amplified only for a small subset of individuals (176 and 58 samples, respectively) belonging to different lineages identified by preliminary analysis. References for primers are given in Supporting Information Table S3 in Appendix S1. For the mitochondrial gene COI, due to difficulty in polymerase chain reaction (PCR) amplification and sequencing for some samples, primer H-t COI (Stuart & Parham, 2004) was used to pair with primer Chmf4 in this study.
All amplifications were conducted in 25 μl volume reactions and were initiated at 95°C for 5 min followed by 35 cycles at 94°C for 1 min, 46°C-55°C for 1 min, 72°C for 1 min, and a single final extension at 72°C for 10 min. Standard annealing temperatures are given in Supporting Information Table S3 in Appendix S1. Negative controls were run for all amplifications. The PCR products were purified with a Gel Extraction Kit (Tsingke Co., Ltd., Beijing, China). Cycle sequencing reactions were performed using a BigDye Terminator Cycle Sequencing Kit (v2.0, Applied Biosystems, USA), and sequencing was conducted on an ABI PRISM 3730XL automatic DNA sequencer (Applied Biosystems, USA) with both forward and reverse primers.
Nuclear gene clonal sequencing was carried out for individuals that contained more than one ambiguous site and failed to infer alleles using PHASE v2.1.1 (Stephens & Scheet, 2005;Stephens, Smith, & Donnelly, 2001). First, the nuclear gene Tyr fragment was re-amplified following the above procedures. The purified PCR products were then cloned into a pClone007 Versatile Simple Vector Kit (Tsingke Co., Ltd., Beijing, China) and transferred into Trelief TM 5α chemically competent cells (Tsingke Co., Ltd., Beijing, China). Plasmids carrying the PCR fragment were then extracted and sequenced. For each PCR product, at least four clones were randomly selected, with one then sequenced.
All newly derived sequences have been deposited on GenBank under the accession numbers MN851304-MN852232 and MN860291-MN860466.

| DNA sequence alignment
DNA sequences were assembled and edited using DNASTAR 5.0 in Lasergene v.7.1.0, aligned using ClustalX 1.81 (Thompson, Gibson, Plewniak, Jeanmougin, & Higgins, 1997) with default parameters, and then examined and revised by eye in MEGA 7.0 (Kumar, Stecher, & Tamura, 2016). Protein-coding nucleotide sequences were translated to amino acids to confirm alignment. Allele sequences of individuals with heterozygous sites were inferred using PHASE 2.1.1, with the algorithm applied five times with different seeds and the assumption of the stepwise mutation mechanism for multiallelic loci relaxed. The input files for PHASE were generated using the web tool SEQPHASE (Flot, 2010). Individuals whose haplotype analysis failed underwent clonal sequencing. The final nuclear gene dataset was the combination of both haplotype estimates and clonal sequencing.

| Datasets
Three datasets were available under study. Two mtDNA datasets were prepared for different analyses. DATA I, for genealogical reconstruction, were the concatenation of 16S rRNA, COI, and CYTB, and contained 195 individuals plus seven outgroup taxa. DATA II, which excluded individuals representing the cryptic species based on genealogical analysis, were the combination of COI and CYTB of 406 samples for genetic analysis. Identical haplotypes were generated using DnaSP 5.10 (Librado & Rozas, 2009). Nuclear gene was the third dataset.

| Genealogical reconstruction and divergence dating estimation
To estimate phylogenetic relationships among mitochondrial haplotypes, Bayesian inference (BI) and maximum parsimony (MP) were implemented using DATA I. The best-fit models of sequence evolution were determined using ModelTest 3.7 (Posada & Crandall, 1998). The BI analysis was performed using MrBayes 3.2.1 (Ronquist & Huelsenbeck, 2003) with six million generations, and two independent runs starting from different random trees were performed with four Markov chains. The chains were sampled every 100 generations. Convergence between the two runs was checked in Tracer 1.7.1 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018). The first 25% of trees were discarded as burn-in, and the remaining sampled trees were used to construct majority rule consensus trees and estimate Bayesian posterior probabilities (BPP) of the tree nodes. The MP analyses were conducted in PAUP4.0b10 (Swofford, 2002) using a heuristic search with 1 000 random-addition sequence replicates, with support for nodes of the resulting MP tree assessed based on 1 000 bootstrap replicates.
Lineage divergence times were estimated in BEAST v1.8.4 (Drummond, Suchard, Xie, & Rambaut, 2012) implemented using the CIPRES Science Gateway portal (Miller et al., 2015). No fossil data were available to serve as internal calibration points within the genus Glyphoglossus or within the closely related genus Microhyla.
Therefore, we assumed a substitution rate of 0.65%-1.00% (mean = 0.8%) per million years (Ma) for CYTB based on evolutionary rates commonly suggested for anurans (Macey et al., 1998(Macey et al., , 2001Monsen & Blouin, 2003). The mutation rate was multiplied by the ratio of the average distance for combined sequences versus that for CYTB alone, after which we deduced the substitution rate of the concatenated fragment (Qu et al., 2011). The net average distance (Da) was estimated to be 0.059 for combined gene fragments, 0.855 times greater than that of CYTB alone. The evolutionary rate of 0.65%-1.00% (mean = 0.80%) per Ma for CYTB was multiplied by a factor of 0.855 to deduce the substitution rate of all fragments combined, 0.556%-0.855% (mean = 0.684%). We employed an uncorrelated log-normal relaxed molecular clock with constant-size tree prior. Analyses were implemented for 100 million generations using the GTR model of nucleotide substitution with gamma-distributed rate variation among sites. Effective sample size for each parameter and convergence were checked using Tracer v1.7.1.
The allele network of the nuclear gene Tyr was constructed.
After generation of a neighbor-joining (NJ) tree based on uncorrected p-distances in MEGA 7.0, the network of haplotypes was visualized in Haploview (Salzburger, Ewing, & Von Haeseler, 2011).

| Molecular diversity and genetic structure
Analysis of DATA II was conducted in ARLEQUIN 3.5 (Excoffier & Lischer, 2010). Haplotype diversity (h) and nucleotide diversity (π) were estimated for overall and for each population with more than five samples, respectively. Geographical partitioning of genetic diversity was inferred via two approaches. First, the G. yunnanensis populations were grouped using spatial analysis of molecular variation (SAMOVA) (Dupanloup, Schneider, & Excoffier, 2002) and were separately analyzed based on mtDNA and nuDNA data. Second, the genetic structure of the populations was inferred by an analysis of variance framework (analysis of molecular variance, AMOVA) (Excoffier, Smouse, & Quattro, 1992 (Kimura, 1980) genetic distance.

| Population demography
The demographic history of each main phylogenetic clade was investigated using Extended Bayesian Skyline Plot (EBSP) analyses (Drummond et al., 2012;Heled & Drummond, 2008) implemented in BEAST2 (Bouckaert et al., 2014), based on mtDNA DATA II (COI + CYTB, 1,140 bp) and nuDNA data (546 bp). The substitution rate (mean = 0.00738) of the mtDNA gene fragment was deduced as in previous BEAST analyses and used as a reference to guide rate estimates for the nuDNA. The analysis was set up following the recommendations of Heled (2015) for BEAST2. Chains were run for 100 million generations sampling every 10,000 steps, and three independent runs were conducted. Final graphs were generated following Heled (2010) with a burn-in cutoff of 20%.

| Sequence characteristics
No premature stop codons were observed in the mitochondrial protein-coding genes, indicating that the sequences were obtained from functional genes rather than nuclear mitochondrial pseudogenes. For the ingroup, the 16S rRNA, COI, and CYTB alignments For Tyr, a total of 58 individuals were sequenced. One to seven heterozygous sites were observed in some individuals. After allele estimates and clonal sequencing, the nuclear gene was 546 bp in length, and 71 haplotypes were generated from 101 allele sequences that included 23 parsimony informative sites.

| Genealogical reconstruction and divergence time estimation
GTR + I + G was selected as the best evolutionary substitution model. The BI, MP, and BEAST trees consistently revealed that the ingroup contained four highly supported major lineages and that these independent evolutionary lineages showed strong geographical structure: that is, LD, W (western), E (eastern), and C (central) (Figures 1 and 2). Lineage LD included haplotypes only from populations 1-2, and 48, and was sister to the group of the other three lineages. Lineage W contained three sublineages (i.e., W1, shared by W/C2, and H2 and H24 by E/C. No alleles were shared between lineages W and E.

| Population genetic structure
The SAMOVA results revealed similar significant F CT values for different numbers (K) of groups (from three to nine) (Table 1). For the mtDNA data, the F CT value reached a plateau at K = 6, but the six suggested groups did not completely correspond to mtDNA lineages ( Figure 2). Localities 3-4, 6, 16-17, and 42-47, geographically corresponding to sublineages C2 and C4, were grouped together (Supporting Information Figure S1 in Appendix S1). Furthermore, the eastern region included two groups: that is, 5, 19, and 24-25 grouped to form E1, and 20-23 and 26-28 grouped to form E2. With the increase in K value, F CT also increased, but to a lesser extent, and new subdivision emerged (Supporting Information Figure S1 in Appendix S1), which were roughly consistent with mtDNA lineages.
The AMOVA results revealed significant geographical structure in the mtDNA genetic variation at all hierarchical levels examined (p < .01; Table 1). When populations were grouped according to the mtDNA lineages (six strongly supported clades W, E, and C1-C4 by three methods) (Figure 2), AMOVA produced the maximum F CT value (0.7186). In contrast, grouping of contemporary drainage basins did not adequately explain the distribution of genetic diversity across the study (F CT = 0.3033).
Pairwise F ST estimates ranged from 0.00 to 1.00, and analysis of genetic differentiations showed strong evidence for correspondence between geographical area and lineage as well as between geographical subdivision and sublineage. Genetic differentiations between geographical zones and between subdivisions were significantly high (Table 2). In contrast, F ST estimates within geographical zones and subdivisions were lower. Therefore, the mtDNA of G. yunnanensis demonstrated considerable geographical structure.

| Genetic structure and historical drainage pattern
The matrilineal genealogy ( Figure 2) revealed clear genetic divergence (W, C, and E) in G. yunnanensis, and further division occurred in each main clade. The analysis of mitochondrial variation by SAMOVA and AMOVA supported this division (Table 1; Supporting Information Figure S1 in Appendix S1). However, our nuDNA data did not find the pattern identified in mtDNA, possibly due to incomplete lineage sorting or implied gene flow. So, we largely restricted this section discussion to the present mtDNA data as they provided a finer geographical resolution compared with the nuDNA analysis.
The geographical distribution of G. yunnanensis exhibits greatest similarity to the montane stream amphibian species N. yunnanensis, and the Yunnan pond frog B. pleuraden. In contrast to montane stream species, however, G. yunnanensis breeds and lays eggs in still water, mainly in the Jinsha River basin. This geographical distribution highlights its use as an ideal model organism for testing scenarios regarding the role of historical drainage rearrangement on species genetic diversity. We may expect its genetic structure to be unaffected by drainage basin history, such as that found for the Yunnan pond frog (B. pleuraden) (Li et al., 2012) and red knobby newt (T. shanjing) (Yu et al., 2013), two other standing-water taxa. To prove this conjecture, we hypothesized that a strong association would exist between the geographical distribution of genetic variation and contemporary drainage basins (H0) and that historical drainage basins would be accountable for the phylogeographic structure (H1). H0 would be rejected by evidence suggesting a phylogeographic break in the genetic relationships of haplotypes from the same drainage basin. Individuals from Jinsha River basin were deeply divergent and scattered in different mtDNA lineages (Figures 1 and 2; see the details below). The AMOVA results (Table 1) showed that differentiation among populations within the same drainage was almost twice that between different drainages (48.93% and 30.33%, respectively) when groupings were based on the current drainage system. These lines of evidence showed there was a phylogeographic discontinuity in the genetic variation within the same drainage basin. Thus, we rejected H0 and supported that contemporary drainage basins are connected by historically independent drainage systems. Groupings based on the mtDNA lineages, which conformed to the hypothesis of drainage system evolution, explained the greatest amount of variation (Table 1) (Wang & Wang, 2005), which is proposed to be associated with the rapid uplift of the QTP in the same period (Cheng et al., 2001;Clark et al., 2004;Yang, Li, Huang, & Ge, 2010). This scenario is further supported by phylogeographic studies on endemic amphibian species and riparian plant species in this region (Zhang, Chen, et al., 2010;.
The above hypothesis is also supported by the presence of three geographically distinct haplogroups (sublineage C1, C2, and C3) in the Jinsha River drainage basin (Figures 1 and 2).  (Clark et al., 2004;Jiang, Wu, Xiao, & Zhao, 1999;Li & Ming, 2011), which resulted in significant changes from lacustrine to fluviatile sediment (Jiang et al., 1999). Sublineage C2 includes four populations from the Jinsha/ Yalong River drainage basin (3-5 and 16) and one population (6) from the Red River drainage basin, which agrees with the hypoth-  These genetic patterns have also been revealed in other endemic species (Li et al., 2012;Yue, Chen, Sun, & Sun, 2012;Zhang, Chen, et al., 2010;. The C3 sublineage mostly occurs in the Lower Jinsha River, but also shows a scattered distribution in the Middle Jinsha River (population 1) (Figure 1). The sharing of some haplotypes among these disjunctive populations (populations 1 vs. 7-14) is difficult to explain, but similar patterns have also been found in riparian plant species such as Terminalia franchetii . Sublineage C3 also consists of populations from Niulan River, and their close genetic relationship is concordant with the hypothesis that the Niulan River-Yongshan-Xinshizhen paleo-drainage was a north-flowing tributary (Li & Ming, 2011;Li, Yang, et al., 2009), different from the south-flowing Yalong and Dadu rivers.
Since the Pliocene, with the exception of the upper reaches in northwest Yunnan, the drainage basins on both sides of the Red River are of antiquity and have experienced little change relative to other regions (Cheng, Chen, Luo, & Peng, 1994;Cheng et al., 2001;Clark et al., 2004), thus resulting in higher genetic homogeneity (Li et al., 2012;Zhang, Chen, et al., 2010;Zhang, Rao, et al., 2010). That is, we would expect populations within the same drainage system are more closely related. Populations from east of Yunnan and from west of the Red River respectively clustered together and formed Lineage E and lineage W (Figure 1). Lineage C4 includes haplotypes from both the Ruili and Salween River basins (Figure 1). The patterns are in accordance with the hypothesis.
C4 is distributed geographically far from the other sublineages of group C, but genetically closer to C2 ( Table 2). The close genetic relationship between these geographically isolated sublineages (C1-C4) indicates a probable connection or no obvious dispersal barriers among these regions during the Pliocene, which coincides with the hypothesis that the YGP formed a planation surface with smooth relief during the middle-late Pliocene (Cheng et al., 2001;He & He, 1993;He et al., 1985). Several phylogeographic studies of endemic species also confirm this inference Zhang, Chen, et al., 2010;. However, more samples from areas between the sublineages are needed to clarify the issue.
The phylogeographical structure of G. yunnanensis substantially corresponds with historical drainage patterns and was largely in concordance with codistributed but disparate Yunnan spiny frog (N. yunnanensis, Zhang, Chen, et al., 2010). Expectation of the concordance of genetic structure among taxa with disparate traits is almost null in most comparative phylogeographic researches (Papadopoulou & Knowles, 2016). Thus, the congruence here might indicate that there was more contribution of the drainage history than ecological or life-history traits in structuring genetic variation between these two disparate codistributed taxa.
However, there are discordant phylogeographic patterns between

| Persistence during Pleistocene climatic oscillations
Cyclical Pleistocene glaciations have had profound influence on population dynamics (Hewitt, 2000(Hewitt, , 2004, whose spatial effects depend on latitude and topography. Population demographic dynamics vary with life history and geography (Hewitt, 2004 (Kozak et al., 2006). Much effort has been made to explore the pattern and extent of Pleistocene climatic oscillations on the geographical distribution and genetic diversity of organisms on the QTP and adjacent areas (Qiu, Fu, & Comes, 2011;Wan et al., 2018;Yan et al., 2013;Yang, Dong, & Lei, 2009;Yu et al., 2013;Yuan et al., 2016).
yunnanensis (Zhang, Chen, et al., 2010), exhibit completely distinct responses to Pleistocene climate fluctuations. In our opinion, ecological opportunity (Wagner, Harmon, & Seehausen, 2012) and microecological niches should also be included to further investigate the responses of organisms to Pleistocene climatic oscillations.

| Secondary contact
Sympatric distribution of distinct lineages was demonstrated in a few populations (Figure 1; populations 3-5, 42, and 45-46), which may result from incomplete lineage sorting, admixture events, or both (Wendel & Doyle, 1998 Figure 1). Furthermore, the higher π (Supporting Information Table S1 in Appendix S1) and absence of intermediate divergent haplotypes are suggestive of secondary contact.
Although nuclear data from Tyr do not well-resolve relationships among G. yunnanensis, they show some noteworthy features. Like the matriline patterns of N. yunnanensis (Zhang, Chen, et al., 2010), sympatric distribution of distinct lineages also occurred in areas adjacent to major drainage divides (Figure 1). For example, lineages C and E co-occurred along the CYP, which forms the Jinsha, Red, and Pearl River watershed (Wang & Wang, 2005) ( Figure 1). This further indicates that the evolution of drainage basins in southwest China has had significant impact on the contemporary geographical distribution of G. yunnanensis genetic variation.

| IMPLI C ATI ON S FOR CON S ERVATI ON
Sympatric distribution of distinct lineages in G. yunnanensis and N.
yunnanensis occurs in many areas. For example, in the CYP, localities 4, 5, and 17 from this study largely correspond to N. yunnanensis sample sites 32, 55-56, and 70 (Zhang, Chen, et al., 2010), which all harbor a particularly high level of nucleotide diversity. The same pattern also appears in western Yunnan for both species. These areas may be major genetic reservoirs for organisms and should be considered for priority protection, pending further monitoring and assessment of additional species.

| CON CLUS IONS
Our study documents the genetic structure within G. yunnanensis.
We found that historical drainage systems are responsible for the divergence between lineages/sublineages. The frog diverged prior to Pleistocene climatic change and experienced a stable demographic history thereafter. Sympatric distributions of distinct lineages resulted from secondary contact, and these areas are major genetic reservoirs and require priority protection. The phylogeographic patterns of G. yunnanensis are almost in concordance with those of N.
yunnanensis, thus showing a shared response to hydrogeological history. Comparative phylogeographic studies of other endemic taxa in southwest China are now required to estimate whether their population genetic structures still preserve evidence of these late Cenozoic changes in palaeohydrology.

ACK N OWLED G M ENTS
We thank Tan

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
DRZ, DQR, and HX conceived the ideas; HH and DRZ carried out fieldwork, with active assistance from XQS, SL, JSW, and SQY. DRZ performed laboratory and data analyses, and drafted the manuscript. GHY helped to revise the manuscript and gave constructive advice. All authors contributed to and approved the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
All sequences were uploaded to GenBank. Additional supporting information may be found online in the Supporting Information sec- Guo-Hua Yu https://orcid.org/0000-0002-0220-6550