Examining the interglacial high‐elevation refugia scenario in East Asian subtropical mountain systems with the frog species Leptobrachium liui

Abstract The effects of Quaternary climatic oscillations on the distributions of organisms in different parts of the world are not equally well understood, limiting the ability to understand the determinants of biodiversity. Compared with the mountain regions in southern Europe and southwestern North America, such effects on high‐elevation species in the East Asian subtropical mountain systems located in southern and southeastern China have seldom been addressed. In this study, using Leptobrachium liui (Megophryidae), we made one of the earliest attempts to examine the interglacial high‐elevation refugia scenario in these Asian mountains. Based on our current understanding of the study system, we formulated a hypothesis that these frogs of western origin were distributed more widely and continuously during glacial phases, allowing eastward dispersal, and that they are currently isolated in interglacial refugia at higher elevations. Microsatellite data and mitochondrial and nuclear sequence data were obtained with extensive sampling followed by the synthesis of phylogeographic and population genetic analyses and modeling of the species distribution. The analyses revealed a sequential eastward divergence of microsatellite clusters and gene lineages accompanied by a decline in genetic diversity. Molecular dating estimates revealed divergence events in the Pleistocene, and a reduction in local populations was inferred to have occurred at a time comparable to the end of the last glacial. Strong genetic isolation by distance reflecting a more continuous historical distribution was detected. Furthermore, environmental niche models inferred a wide planar distribution during the last glacial maximum, providing further support for the hypothesis.

Relevant studies in this region are rare and have primarily focused on southwestern China and adjacent areas with their complex topography and high mountains, such as the Himalayas (e.g., Gao, Zhang, Gao, & Zhu, 2015;He & Jiang, 2014).
East of these areas, numerous subtropical mountains stretch across southern and southeastern China. Their most recent rapid uplifts have been dated to before the Quaternary or during the early part [2.6-0.8 million years ago (Ma)] of the Quaternary (Guo, 1998;Zhou & Ren, 1984). With the highest peaks reaching ~2,100 m above sea level, these mountains are surrounded by lowlands ~200-400 m in elevation and exhibit strong temperature and habitat gradients ( Figure 1). The limited available evidence suggests that certain coldadapted species ranged widely during glacial stages but became isolated at high elevations in these mountains during the present interglacial. First, pollen fossils of some current upland taxa were common at lower elevations during the last glacial maximum (LGM) (Zheng, 2000). Second, using mitochondrial DNA data, Tian, López-Pujol, Wang, Ge, and Zhang (2010) identified several interglacial montane refugia for the pine species Pinus kwangtungensis and noted another three plant species with similar distributions and genetic patterns (but see Qiu, Fu, & Comes, 2011). Last, using molecular dating, Wu, Wang, Jiang, and Hanken (2013) found that all species of the endemic salamander genus Pachytriton had diverged before the Quaternary, and they reported that these cold-adapted animals are currently restricted to high-elevation refugia. Many other taxa that currently live at high elevations in these mountains may have experienced similar evolutionary histories.
Pope's spiny toad, Leptobrachium liui (Pope, 1947) (Figure 2), is ideal for evaluating this possibility. Leptobrachium liui is restricted F I G U R E 1 Sampling sites and geographic pattern of expected heterozygosity (H e ) inferred from the 8-loci microsatellite data. Sampling sites were grouped into five distribution areas, which were then grouped into two ranges, west (W) and east (E). W W : the most western area in which L. leishanense occurred; W M : middle W; W E : eastern W; E W : western E; E E : eastern E to habitats at middle to high elevations (700-1,600 m) and is found in several major southern and southeastern China mountain ranges that are separated by low elevations or connected only by middleelevation passes and that form the northeastern range limit of the genus. This distribution may have evolved via west-to-east dispersal followed by long-term occupation, as suggested by the deep eastward divergences among regional mitochondrial lineages (Kropf, Kadereit, & Comes, 2002;Waters, Fraser, & Hewitt, 2013;Zheng, Li, & Fu, 2008). In a much smaller range west of and adjacent to the distribution of L. liui, its sister species L. leishanense is restricted to similar habitats at 1,050-1,800 m in elevation (Zheng et al., 2008). These frogs breed in cold permanent streams from October to December, the hatching period extends over the entire winter (Li et al., 2009), and the tadpoles usually require 3 years to reach metamorphosis.
Two factors suggest that the current disjunct distributions of these frogs were not caused by human activities (Sanderson et al., 2002), which are more intensive at low elevations. First, certain distribution localities were identified by systematic herptile surveys from low-to-high elevations conducted no later than the early 1960s (e.g., Liu & Hu, 1962) during a time when pollution was not a significant issue in the region. Second, no low-elevation populations have been recorded in reserves providing natural habitats for L. liui or L. leishanense.
There is an ambiguous species boundary between the two taxa.
The only morphological difference between the species is an internal subgular vocal sac in male L. liui that is not found in L. leishanense; however, both species call underwater (Zheng, Rao, Murphy, & Zeng, 2011), and intraspecific variation in the vocal sac condition has been reported in certain anurans (Che et al., 2009;Yu, Zhang, & Yang, 2010). Currently, adequate molecular evidence distinguishing these two species is lacking. In a nuclear gene tree created by Zheng et al. (2008), most haplotypes of L. liui are regionally endemic, while L. leishanense shares its haplotypes with western L. liui. In a mitochondrial gene tree created by Matsui et al. (2010), L. leishanense is nested within lineages recognized as L. liui in Zheng et al. (2008) andFrost (2018). Therefore, L. leishanense may be L. liui populations occupying the western edge and a possibly ancestral range of the generally east-west distribution, and including L. leishanense is imperative when testing a hypothesis on the historical distribution of L. liui.
The interglacial high-elevation refugia scenario includes a contraction after the end of the last glacial ~0.01 Ma, which is expected to have caused a reduction in population size. Rapidly evolving and hence highly polymorphic microsatellites may provide adequate data to assess whether there has been a recent decline in population size (Peery et al., 2012). Approximate Bayesian computation (Beaumont, 2010) estimates the parameters and relative posterior probabilities for a set of user-specified demographic models, which can be complex (e.g., >1 population size changes). Range expansions of limited extent might have occurred in some cold-adapted species after the warmest period of the present interglacial, the Holocene climatic optimum (Bell et al., 2010), and resulted in detectable demographic expansions. A model defining an increase after a decline in population size can be used to account for this possibility. In southern China, the optimum was dated to approximately 10,000-6,000 years ago, and it was followed by a decrease in the extent of forests . Another possible demographic change that may complicate the inference is a more recent decline, within hundreds of years, caused by human impact. This can be considered at first by a simple model of one change/decline in population size at a time to be estimated. With immediate effects on microsatellite variation and limited recovery time, the sign of a newly occurred reduction would be readily detected (Cornuet & Luikart, 1996;Garza & Williamson, 2001;Hoban, Gaggiotti, & Bertorelle, 2013).
Under the interglacial refugia scenario, L. liui may exhibit an isolation by distance (IBD) pattern of genetic structure across the range.
If these frogs ranged more widely and continuously during the glacial phases, their distribution over the long term should have been more continuous, with gene flow interrupted only by short interglacial disjunction into high elevations (Schmitt, 2007). Change in the genetic differentiation between populations often requires long time periods (Ibrahim, Nichols, & Hewitt, 1996;Leng & Zhang, 2011;Palo et al., 2004). It is unlikely that an IBD pattern established over long periods of a continuous distribution (Slatkin, 1993;Wright, 1943) would have been erased by the increased genetic differentiation during the current interglacial of ~0.01 million years (Mastretta-Yanes et al., 2018). In theory, the effect of such increased genetic differentiation on the Mantel correlation test for IBD may be analogous to adding a constant to the genetic-distance variable, which will not affect the detection of a correlation. If an IBD pattern is detected, the genetic distances are expected to be better predicted by log-transformed geographic distances (two-dimensional IBD) than by raw geographic distances (one-dimensional IBD) when the habitat is more planar (Rousset, 1997) over the long term.
F I G U R E 2 Male Leptobrachium liui with maxillary nuptial spines in Huaping Reserve, Guangxi, China. Photograph credit: Yuchi Zheng The predicted planar distribution during glacial stages can also be examined by correlative species distribution modeling. Specifically, ecological niche models are used to generate the distribution projection based on reconstructed LGM climatic variables. Different variables reconstructed using different models can result in dissimilar and even contradictory projections with respect to a given question of interest. In such situations, molecular evidence can be useful in interpreting the projections (Gavin et al., 2014).
In this study, using extensive sampling of both species with 907 individuals from 25 sites, 14 microsatellite loci, one mitochondrial locus, and one nuclear locus, we tested a hypothesis integrating eastward dispersal with the interglacial high-elevation refugia scenario. We hypothesised that L. liui originally occurred in the west, and spread from high-elevation refugia to a more continuous and planar habitat during Quaternary glacial periods that allowed dis- (e) two-dimensional IBD fitting the genetic structure better than one-dimensional IBD if an IBD pattern was found; and (f) a wider and more continuous distribution during the LGM inferred by correlative species distribution modeling. We examined these predictions based on the combined data of L. liui and L. leishanense because we detected lower genetic differentiation between L. leishanense and nearby L. liui than between the latter and the other populations of L. liui, suggesting that L. leishanense is a synonym of L. liui. The microsatellite loci included VIB-B4, VIB-C10, VIB-D5, CHA1, CHA6, CHA9, LSMT6, LSMT8, CHA4, CHA7, CHA8, CHA11, LSMT4, and LSMT9 (Bi, Deng, Crosby, & Fu, 2010;Hu, Xia, Zheng, & Zeng, 2012;Wang et al., 2011). The PCR fragments containing primers fluorescently labeled with FAM, TAMRA, or HEX were evaluated along with the GS500 size standard using an ABI 3730 analyzer (Applied Biosystems, Foster City, CA, USA). The raw data were visualized, and the alleles were manually scored using GeneMarker 1.95 (SoftGenetics, State College, PA, USA). Tests for the presence of stuttering, allele dropout, and null alleles were performed using Micro-Checker 2.2.3 (van Oosterhout, Hutchinson, Wills, & Shipley, 2004). An exact test of Hardy-Weinberg equilibrium and a likelihood ratio test for linkage disequilibrium were conducted in GENEPOP 4.2 (Rousset, 2008). The 745-bp mitochondrial nad1 locus consisted of regions of the tRNA-Leu (16 bp) and NADH dehydrogenase subunit 1 genes; the nuclear rag1 locus contained an 800-bp partial sequence of the recombinase activating protein 1 gene. The primers used in the PCR and Sanger sequencing (ABI 3730) experiments included Leu-1L, Leu-5L, Rag1-3L, Rag1-2H, Rag1-6H (Zheng et al., 2008), and a newly designed ND1-6H (5′-TgACCCTAAgAACAgAATT

| Sampling and molecular data
ACTgATAAggT-3′). The rag1 PCR products were directly sequenced, and when two or more dimorphic sites resulting from heterozygosity were detected, 3 to 16 (mean 9.2) clones were sequenced to distinguish among the haplotypes.
Three congeners were used as outgroups for reconstructing gene trees. In the reported mitochondrial gene trees, L. boringii belonged to the clade (or one of the clades) most closely related to L. liui and L. leishanense (Matsui et al., 2010;Zheng et al., 2008).
Thus, L. boringii was selected as the outgroup for the nad1 data. In the reported rag1 tree, the lineages of L. liui and L. leishanense did not coalesce until the most recent common ancestor of L. ailaonicum, L. boringii, L. chapaense, L. leishanense, L. liui, and L. promustache because of incomplete lineage sorting (Zheng et al., 2008). The sister group of this clade contained L. ngoclinhense and L. xanthospilum, which were used as outgroups for the rag1 data. The nad1 sequences of 29 individuals and the rag1 sequences of 27 individuals were obtained from Zheng et al. (2008). Detailed specimen information, GenBank accession numbers, and haplotype designations are presented in Supporting Information Table S1.
Sequence alignments were performed with ClustalX 2.1 (Larkin et al., 2007) and checked against the amino acid sequences. The Bayesian information criterion implemented in jModelTest 0.1.1 (Posada, 2008) was used to select the nucleotide substitution model that best fit the data.

| Microsatellite clustering, gene tree, and divergence date
The clustering method implemented in STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000) was applied to the microsatellite data. The admixture ancestry model was adopted, with the value of the Dirichlet parameter α for the degree of admixture inferred from the data. The prior population information was not used. The correlated allele frequencies model was used, with the value of the allele frequencies parameter λ fixed to 1. The burnin length and the run length after burnin were set to 100,000 and 200,000 steps, respectively. The number of clusters (K) evaluated ranged from one to the number of sampling sites included plus one. For each K, 20 independent runs were conducted, and the 10 runs with the highest likelihoods were included in the follow-up analyses. In STRUCTURE HARVESTER 0.6.94 (Earl & vonHoldt, 2012), the results were summarized and compared using the Evanno method and the Delta K (Evanno, Regnaut, & Goudet, 2005) was calculated. Combining runs with the same mode of clustering, a graphical representation of the Q-matrix was generated using CLUMPAK (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015). The following four criteria were considered in detecting values of K that captured the major structures in the data: Delta K, a high rate of change of the likelihood, the most individuals from a sampling site being assigned to a population with high possibility (say ≥0.85), and a reasonable biological interpretation (Meirmans, 2015). When hierarchical clusters were obtained using the selected K values, they were summarized as a tree-like population history.
Gene trees of the unlinked nad1 and rag1 loci were separately reconstructed because the relationships among haplotypes were of primary interest. Using the nad1 data, divergence times were estimated simultaneously with Bayesian topology in BEAST 2.1.2 (Bouckaert et al., 2014). The monophyly of the ingroup was constrained. The uncorrelated lognormal model was used to describe the relaxed clock. The calibrated Yule model was specified for the tree prior. Suitable fossil calibrations were not available. An evolutionary rate of 0.69% ± 0.3% divergence per million years per lineage that is commonly used for the anuran nad1 gene was used to calibrate the tree (e.g., Hoffman & Blouin, 2004;Macey et al., 1998;Sanguila, Siler, Diesmos, Nuneza, & Brown, 2011). The uncorrelated lognormal relaxed clock mean (ucldMean) was constrained between 0.0039 and 0.0099 substitutions per site per million years, and a normal distribution was applied (mean 0.0069, sigma 0.0015). For rag1, Bayesian inference was performed with MrBayes 3.2.1 (Ronquist & Huelsenbeck, 2003), using four chains. In both BEAST and MrBayes analyses, a run length of 20 million generations and a sampling frequency of 2,000 generations were applied, and four independent runs were performed to avoid local optima. The performance of the runs was assessed using Tracer 1.6 (Rambaut, Suchard, Xie, & Drummond, 2014) to ensure convergence and adequate effective sample sizes (>200; Drummond, Ho, Phillips, & Rambaut, 2006).
Then, the last 90% of the samples for each of the four runs were combined to produce a maximum clade credibility tree.

| Genetic diversity and structure
Diversity parameters were estimated for each local population, and the inverse distance-weighted method incorporated in ArcGIS 9.2 (ESRI, Redlands, CA, USA) was used to present spatial patterns in diversities. For the microsatellite data, the expected heterozygosity was estimated using Arlequin 3.5.1.2 (Excoffier & Lischer, 2010), and allelic richness and private allelic richness were calculated with rarefaction to the smallest sample size of 8 using HP-Rare 1.1 (Kalinowski, 2005). The haplotype diversity and average number of pairwise nucleotide differences were obtained from the nad1 and rag1 data using DnaSP 5.10.01 (Librado & Rozas, 2009).
Various measures of population differentiation were calculated and used for further analysis. For the microsatellite data, the F ST value was estimated according to Weir (1996) using the software FreeNA (Chapuis & Estoup, 2007), applying the ENA method for correcting genotype data with null alleles, and the D est value (Jost, 2008;Leng & Zhang, 2011) was estimated using SMOGD 1.2.5 (Crawford, 2010). For the nad1 and rag1 data, analysis of molecular variance (AMOVA) was conducted using Arlequin. To assess the major pattern of differentiation between regional population groups, principal coordinate analysis (PCoA) was conducted on individual interpopulation difference matrices using GenAlEx 6.502 (Peakall & Smouse, 2012), with negative values set to zero. The Isolation by Distance Web Service 3.23 (Jensen, Bohonak, & Kelley, 2005) was used to perform the Mantel test for correlations between the genetic (F ST , D est , or AMOVA Φ ST ) and geographic distances, and 10,000 randomizations were applied, with the negative genetic distances set to zero. The geographic distances were calculated using coordinates and analyzed with and without logarithmic transformation.

| Population size change
Addressing the predicted recent decline ~0.01 Ma in each local population, the approximate Bayesian computations were applied to the microsatellite data using DIYABC 2.1.0 (Cornuet et al., 2014). Five demographic scenarios were compared. In the one-change scenario, a change in population size occurred sometime within 10-10,000 generations ago. Because the tadpoles of these frogs require 3 years to reach metamorphosis, a generation time of 4 or more years could be reasonably assumed. Consequently, the prior time extended to no <40,000 years ago. The following three scenarios were special cases of the one-change scenario: Population growth occurred in the onegrowth scenario, no changes (rate = 0) in population size occurred in the constant-size scenario, and a decline in population size occurred in the one-decline scenario. Lastly, in the growth-after-decline scenario, two changes in population size occurred 10-10,000 generations ago; a decline first occurred from a population size larger than the current one, and then, an increase in population size occurred.
This scenario was designed to model a decline after the LGM followed by growth after the Holocene climatic optimum, with the resulting current population remaining smaller than that at the LGM.
For all scenarios, the prior effective population size was set to 10-10,000, and uniform distributions were applied to the time and population size priors. All loci were included in a single locus group.
With default priors valid for many species (Cornuet et al., 2015), a generalized stepwise mutation model allowing single-nucleotide indels was applied to the data. All four summary statistics for one sample of microsatellite data were used. For each scenario, 10,000,000 simulations were performed to build the reference table. Among them, the 10,000 (0.1%) with summary statistics closest to the observed ones were used to estimate the historic demographic parameters. The relative posterior probabilities of different scenarios were computed using the direct estimate and logistic regression approaches with 500 and 10,000 simulated datasets closest to the observed dataset, respectively (Beaumont, 2010). For analyses evaluating confidence in the scenario choice, the posterior-based error computation was adopted. Computations were processed on 1000 simulated pseudo-observed test datasets, each of which drew (with replacement) the scenario ID and parameter values from the 500 simulated datasets closest to the observed one. Then, the above procedure for computing posterior probabilities of scenarios was conducted on each pseudo-observed test dataset, and the right scenario and the scenario having the highest probability were reported.

| Species distribution modeling
The LGM and current species distributions were modeled for comparison. The analyses were conducted at a spatial extent of 21.25-31°N and 106.5-122.25°E. This area did not extend substantially from the ranges of the study species and the mountains (Figure 1).
The current range of the species was expected to be indicative of its LGM distribution, which was fragmented into various present-day local refugia according to our hypothesis. Ecological niche modeling was conducted using Maxent 3.3.3k (Phillips, Anderson, & Schapire, 2006). The Maxent model, which works by optimizing a set of constraints representing the incomplete information about distribution and evaluating the environmental suitability of each grid cell within the study area (Elith et al., 2011;Phillips et al., 2006), has been shown to have good predictive performance across various applications (Elith et al., 2006). Nineteen bioclimatic variables were obtained from the WorldClim database (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), with a resolution of 2.5 arc-minutes. Georeferenced occurrence records were obtained from our field surveys, and only one occurrence per grid cell was ensured. The final compilation included 29 records from both taxa (Supporting Information Table   S2). There is a need for variables to be as proximal as possible, and over-fitting could occur if there are too many strongly collinear variables, especially with the small sample size of the present records (Qiao, Peterson, Ji, & Hu, 2017;Synes & Osborne, 2011). We hence conducted a Pearson's correlation test between variables (Dormann et al., 2013) and a jackknife analysis of variable importance to the Maxent model . Two criteria were considered for choosing between highly correlated temperature/precipitation variable pairs (|r| > 0.8), namely, larger contribution to model development and putative greater biological importance. As a result, 10 variables were retained as follows: annual mean temperature, mean monthly temperature range, annual temperature range, mean temperature of the wettest quarter, mean temperature of the warmest quarter, mean temperature of the coldest quarter, annual precipitation, precipitation of the driest month, precipitation seasonality, and precipitation of the warmest quarter.
In the Maxent analysis, we mainly used the default settings and ran models with 10 bootstrap replicates by randomly assigning localities to the training (75%) and test (25%) data sets. The default settings have been validated over a wide range of species using different sets of environmental variables in various regions of the world and have been shown to achieve good performance (Phillips & Dudík, 2008). The replicates from Maxent were used as proxy models to develop the consensus-based ensemble forecasting (Araújo & New, 2007), and the mean suitability of the replicate outputs was calculated using ArcGIS 9.2 (ESRI). The easily interpretable logistic output format was selected with suitability values ranging from 0 (lowest) to 1 (highest) (Phillips & Dudík, 2008). Model performance was assessed using the average AUC (area under the receiver operating curve) score. The most important variables to the model were identified by their percent contribution values and the results of the jackknife tests (Phillips, 2017).
How they affected the distribution was interpreted with the aid of the response curves made using only the corresponding variable.
The paleodistributions in the LGM were inferred with the same resolution of 2.5 arc-minutes (Hijmans et al., 2005). The climatic variables of the LGM reconstructed by the community climate system model (CCSM; Collins et al., 2006) and the model for interdisciplinary research on climate (MIROC; Hasumi & Emori, 2004) were used. For each important variable, the LGM reconstructions and current estimates were compared to characterize the distribution and extent of possible novel climate conditions, a potential issue of model projecting (Phillips, 2017).

| Molecular data
Eight microsatellite loci, VIB-B4, VIB-C10, VIB-D5, CHA1, CHA6, CHA9, LSMT6, and LSMT8, were successfully genotyped for all sampling sites included, forming the 8-loci data used for interpopulation analyses and comparisons. From the 14-loci data with missing genotypes, Micro-Checker did not return evidence of stuttering or allele dropout and detected 13 cases (4.8%) of potential null alleles (all p < 0.05), with one to three cases per locus.
After sequential Bonferroni corrections, four cases (1.5%) of significant (0.05 level) deviation from Hardy-Weinberg equilibrium were identified in three populations, and 31 departures (2.7%) from linkage equilibrium (one or two cases per locus pair) were significant at the 0.05 level.
The nad1 dataset had 122 haplotypes and 745 nucleotide sites of which 184 sites were variable and 138 were parsimony-informative among the ingroup members. Among the 138 sites, 0, 19, and 2 sites belonged to the tRNA-Leu gene, the 1st codon position, and the 2nd codon position, respectively. The rag1 dataset had 43 haplotypes and 800 sites, of which 33 were variable and 16 were parsimonyinformative among the ingroup members. Of those 16 sites, 2 and 3 were from the 1st and 2nd codon positions, respectively. To avoid potential overparameterization of the intraspecific analysis, we did not divide the nad1 and rag1 datasets into partitions. For gene tree reconstruction, the TrN+G model was selected for the nad1 data, and the K80 + I model was selected for the rag1 data.

| Microsatellite clusters, gene trees, and divergence dates
The STRUCTURE clustering analysis detected sequential eastward divergences and grouped L. leishanense (W W ) with nearby L. liui (W M ) (Figures 1 and 3; Supporting Information Table S3). Using the 8-loci dataset containing all sampling sites, K = 2 was selected (mean α = 0.031). The samples from areas W W and W M and those from W E and E were separately assigned to two groups, which were then ana-  (Pritchard, Wen, & Falush, 2010). Within the W E +E group, a tree-like population history was inferred. For K = 2, the individuals from E were separated from those from W E . For K = 3, the individuals from E E were further separated from those of E W . When K = 4, within E E , the individuals of sites 1 and 2 were separated from those of sites 3 and 4. Sequential eastward divergences that occurred approximately 3.2-0.4 Ma and, more generally, strong and long-standing geographic structures were revealed by the nad1 time-calibrated tree (Figure 4).
The nad1 tree exhibited mostly well-supported major lineages consistent with the mitochondrial topology in Zheng et al. (2008) with fewer localities but additional genes (2,577 sites). It could be summarized as (W W ,(W M ,(W E ,E))), with a 1.5-0.6 Ma divergence date between W E and E. In E, the E E lineage was nested within the E W lineages and was estimated to have diverged from the latter 0.6-0.2 Ma. The distribution could be divided into eight areas/subareas (W W , four W M subareas, W E , E W , and E E ) that correspond exclusively to various genetic lineages that diverged approximately 4.2-0.4 Ma (Figures 1 and 4). Within the areas/ subareas, each main subclade was formed by all haplotypes from one or a few adjacent local populations, and only populations close to one another shared haplotypes.
The rag1 tree showed eastward divergences and close relationships between W W and W M haplotypes (Supporting Information Figure S1). Four haplotypes were observed in the eastern range E, one in E E and four in E W . They were most closely related to those found in the nearby western populations 8-11 (W E ) and 12, forming two separate lineages. These and another two W E lineages were nested within several W M and W W lineages.

| Genetic diversity and structure
A general pattern of a decline in genetic diversity was detected moving eastward along the distribution, with low levels of genetic diversity found in both the most southern (22) and northern (20) sites (Supporting Information Table S4). The spatial pattern of the expected heterozygosity for the microsatellite data is presented in Figure 1 as an example.
The close relationship between W W and W M populations was also revealed by PCoA and AMOVA of the microsatellite and rag1 data. The samples from W M + W W , W E , and E were easily distinguished in the PCoA plots, and the W W points were mixed with or located very close to the W M points ( Figure 5 and Supporting Information Figure S2). The rag1 AMOVA Φ CT estimates were 0.370 (p < 0.05) for the group pair W M and W E + E, 0.047 (p = 0.198) for W M and W W , and 0.422 (p < 0.05) for W W and W E + E. These lines of comparison had not been applied to the nad1 data, in which the gene tree showed a sister-group relationship between W W and the others.

F I G U R E 3
Bayesian clustering for all samples (K = 2) and for the individuals from sites 1-11 (K = 2 to K = 5) based on the 8-loci data. High Delta K values were obtained for these K priors. Each individual is represented as a vertical bar colored in proportion to its estimated membership in each cluster. The top phylogenetic hypothesis of clusters corresponding with individual distribution areas was summarized from these results 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 1 3 1 5 1 6 1 7 1 8 1 9 2 0 2 1 2 2 2 3 2 4 2 5 A strong two-dimensional IBD model fitted the genetic variation better than the one-dimensional IBD model across the range.
Genetic variation was explicitly spatially structured and low levels of differentiation were limited to local scales (Supporting Information   Table S5). For the microsatellite data, the mean F ST and D est estimates were 0.31 and 0.63, respectively. For rag1 and nad1, as   revealed by AMOVA, the among-population variation accounted for 54.7% and 95.2% (both p < 0.05) of the total variation, respectively.
These genetic variations could be primarily explained by logarithmic geographic distances, which consistently performed better than the linear geographic distances (Table 1).

| Population size change
The DIYABC analysis of the microsatellite data suggested that a slight reduction in population occurred thousands of generations ago in local populations (Supporting Information Table S6). Four sites were combined into two, 1 + 2 and 3 + 4, given the nonsignificant F ST and short geographic distances. For each population, all available loci (8 to 14, mean 11) were used. Fairly consistent results of scenario choice were produced by the direct and logistic regression approaches. The one-decline scenario resulted in the highest probabilities in 21 (95.5%) and 19 (86.4%) of the 22 populations in the direct and logistic regression approaches, respectively. In the other cases, the growth-after-decline scenario resulted in the highest probabilities for population 3 + 4 for both approaches and populations 5 and 22 for the logistic regression approach. In population 3 + 4, if the growth-after-decline scenario was not considered, the most likely scenario would have been distinct growth ~3.6 thousand generations ago. It was not surprising that the one-change scenario was not the most likely in any comparisons because all the special cases of it were covered in other scenarios. However, demographic parameter estimates under the one-change scenario and the most likely scenario were highly consistent, suggesting that the data had sufficient information to distinguish between these simple scenarios (Supporting Information Table S6). On the other hand, the estimated confidence in scenario choice among the one-growth, constant-size, one-decline, and growth-after-decline scenarios was not high. The mean error rate of the direct approach was 0.429, with a mean of 46.3% cases of the best scenario for the observed data being falsely selected for the pseudo-observed test data and a mean of 30.7% cases of the best scenario being falsely considered as suboptimal.
For the logistic regression approach, the corresponding values were 0.427%, 40.8%, and 34.4%, respectively. These evaluations, although they might be influenced by the number of scenarios compared (Cabrera & Palsbøll, 2017), implied that the information was not sufficient for determinable scenario-comparison and precise parameter estimation in individual populations. Interpreting the results from each population jointly, a decline of 1/4-fold and three thousand generations ago fitted the data best on average, and the scenarios of a more recent decline or a commonly occurring subsequent growth after decline were not supported.

| Species distribution modeling
Compared with the current distribution, the inferred distribution of these frogs during the LGM was more continuous and planar. The Maxent model provided reasonable discrimination with high AUC values, with an AUC training of 0.975 ± 0.007 (mean ± SD) and an AUC test of 0.971 ± 0.014. Two extreme variables, T war (mean temperature of the warmest quarter; 66.8%) and Prec dry (precipitation of the driest month; 17.0%), provided the greatest contributions to model development (Supporting Information Table S7). Additionally, the jackknife tests revealed that T war was most informative and that Prec dry contained the greatest amount of information that the other variables could not provide. Their response curves (Supporting Information Figure S3) indicated that high maximum temperature and low minimum precipitation severely limited the distribution of  these cold-adapted frogs, the tadpoles of which live in permanent streams, which was biologically reasonable. As a result, the present disjunct montane distribution was accurately inferred (Figure 6).
These performances suggested that the reliability of the LGM projection would be unlikely to suffer from low reliability of the Maxent model itself. Although the LGM potential distributions based on the CCSM and MIROC reconstructions did not match one another closely, they both overlapped with the major extent of the present distribution and were largely located within the present latitudinal limits (Figure 6). Such a comparison could be obtained using projections in the mainland portion of the study area. In the present mainland, for T war and Prec dry , the only LGM reconstruction outside the training range was a marginally lower minimum precipitation according to MIROC in the most southwestern corner (Supporting Information Figure S3). The response curve showed that the probability of the presence had already declined to near zero at the dry limit of the minimum precipitation range. Therefore, in this case, little effect of the Maxent clamping strategy, treating novel conditions as the same as the limit of the training range, was expected (Phillips, 2017). As T war and Prec dry contributed 83.8% to model development, we also expected limited, if any, complications caused by novel conditions of the other variables in examining the prediction of a wider and more continuous LGM projection.

| D ISCUSS I ON
Our hypothesis that L. liui of western origin responded to Q uaternary climatic cycles by spreading from high-elevation interglacial refugia to more continuous and planar habitats during glacial periods, which allowed eastward dispersals that extended the distribution range to the current one, was not rejected. A closer relationship between L. leishanense and nearby L. liui than between the latter and the other L. liui populations was revealed by PCoA and AMOVA as well as clustering and gene tree analyses of the nuclear microsatellite and rag1 data (Figures 3 and 5). Given this finding and a lack of morphological diagnostic characters, the predictions of the hypothesis were examined using the results from the combined data of the two putative species. The results are consistent with the predictions and provide details for the hypothesized scenario.
Concordant patterns of sequential eastward divergence were separately observed in the microsatellite clusters and gene trees.
The microsatellite clusters and nad1 lineages corresponding exclusively to regions W E , E W , and E E diverged eastward from the western W W and W M regions (Figures 1, 3, and 4). A lower level of turnover was presented by the rag1 tree than by the nad1 tree (Supporting Information Figure S1). Similar differences between nuclear and mitochondrial loci have frequently been reported and have been explained, for instance, by the smaller effective population size of the mitochondrial marker and male-mediated gene flow (Knopp & Merila, 2009;Moritz, Dowling, & Brown, 1987). Despite this difference, the same order of divergence was supported by the nesting of rag1 lineages. The single E E rag1 haplotype was identical to one of the four E W haplotypes, most of which were nested within one of the three W E lineages, which in turn were nested within a number of W M and W W lineages. The other E W haplotype was nested with a haplotype from site 12, which was located close to W E and showed some admixture of eastern microsatellite alleles as well (Figure 3). The predicted eastward decline in variation during dispersal was detected in all types of genetic data (Figure 1; Supporting Information Table   S4 The eastward divergences were estimated to have occurred mainly in the Pleistocene (2.588-0.012 Ma) (Figure 4). we believe that the result agrees with the prediction.
As predicted, the two-dimensional IBD model fitted the genetic variation better than the one-dimensional IBD model across the range (Table 1). In favor of the prediction, with a coefficient of determination (R 2 ) up to 0.630, the two-dimensional IBD pattern detected was fairly strong. Various measures of population differentiation also showed generally high levels of genetic turnover (Supporting Information Table S5). This finding is consistent with the long-term habitation suggested by the nad1 time tree and the hypothesized interglacial disjunction (Wielstra et al., 2013).
In general, the inferred distributions during the LGM were more extensive and continuous than the current distribution and were located well within the present latitudinal range. Between the CCSM and MIROC projections, a notable difference was observed around site 22. This local population exhibited exclusively unique sequence variation, suggesting long-term persistence (Figure 4 and Supporting Information Figure S1). Such independent evidence was consistent with the CCSM projection but inconsistent with the MIROC one, which showed that site 22 and a nearby region would be unsuitable in the LGM ( Figure 6). Consequently, the colder and wetter CCSM reconstructions around site 22 might provide a better fit to the two variables that contributed the most to the Maxtent model (Supporting Information Figure S3). In addition to the caveat that these projections should be viewed with caution, this finding again emphasized the various capabilities of different reconstructions when applied to specific parts of the world (e.g., Arenas, Ray, Currat, & Excoffier, 2012;Igea et al., 2013). In recent years, an increasing number of studies have simultaneously used species distribution modeling and molecular and/or fossil data (Gavin et al., 2014). When meta-analyses of such studies reveal that a climatic variable reconstructed using a particular model is too high/low in a specific region, the finding may be used to help interpret projections and determine model modifications (McGuire & Davis, 2013).
The findings of this study are among the first pieces of evidence that certain cold-adapted organisms in the subtropical mountains of the East Asian mainland responded to Quaternary climatic oscillations by expanding their ranges during glacials and retreating to high elevations during interglacials. Given this response, the current ranges of these frogs are refugia used for various lengths of time.
These ranges overlap extensively with the sampling of previous studies suggesting the same scenario, with Tian et al. (2010) suggesting refugia for pine trees in six mountains (sites 12, 14, 17 + 18, 19, 21, and 23) in the west and Wu et al. (2013) suggesting refugia of salamanders in eight mountains (sites 1 + 2, 5 + 6, 8 + 9, 10, 12, 17 + 18, 21, and 25) across the range. The eastern range (E) was inferred to be established through dispersals that occurred in different Quaternary glacial periods, forming the northeastern limit of the range of this primarily tropical genus of 35 described species (Frost, 2018). In shaping diversity patterns, the contribution of the glacial climate through facilitating the establishment of new distributions for cold-adapted organisms remains unclear. Modern climate change, on the contrary, is likely a threat to many cold-adapted species (e.g., Gottfried et al., 2012). Our results support that L. liui is isolated in sky-island refugia by warm low-elevation habitats. Subtropical East Asia is expected to be warmer in the 21st century (Collins et al., 2013). This trend may lead to habitat reduction and local extinctions of this species, for which the conservation efforts can benefit from the distribution of genetic diversity and evolutionary processes inferred in the present study (Moritz, 2002).
Since this work was submitted for publication, Li et al. (2018) reported a phylogeographic study including 13 L. liui and one L. leishanense sampling sites and using mitochondrial cytb and nad4 genes and eight microsatellite loci. Based on the time-calibrated gene tree, the authors proposed that the late Miocene intensification of East Asian monsoons Xia et al., 2009), the warmer climate of the Pliocene (Jiang, 2009;Salzmann et al., 2011), and an early Pleistocene uplift of the Wuyi Mountains (Liu, 1984) shaped the genetic structure of L. liui, providing a hypothesis for the causes of its early dispersal and divergence. In addition to molecular dating, they also assessed regional genetic diversity, microsatellite clustering with STRUCTURE, and IBD. The mitochondrial time-calibrated tree of Li et al. (2018) (Rousset, 1997;Whitlock & Mccauley, 1999) and the smaller sample size (n = 13) might have been relevant to the failure to detect a strong IBD pattern across the range of L. liui. More generally, these differences again demonstrated the substantial effects of the geographic sampling scheme on spatial genetic structure analyses.

ACK N OWLED G M ENTS
We thank John Wiens, Jinzhong Fu, Qiang Dai, Jeffrey Streicher, Martha Muñoz, Shea Lambert, and Tereza Jezkova for their suggestions and comments. We also thank Aurélie Bonin and five anonymous reviewers for their helpful comments on the manuscripts. This work was supported by the National Natural Science Foundation of China (30900134, 31270568, 31272282, 31372181, 31572243) and the Chinese Academy of Sciences (09C3011100, Y1C3051100).

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
Y.Z. and X.Z. conceived the study. Y.Z. conducted the field and molecular lab work. Y.Z. and J.H. performed the analyses. Y.Z. led the writing with assistance from J.H. and X.Z.