A Phytogeographic Divide Along the 500 mm Isohyet in the Qinghai-Tibet Plateau: Insights From the Phylogeographic Evidence of Chinese Alliums (Amaryllidaceae)

The Qinghai-Tibet Plateau (QTP) has been biogeographically divided into the eastern monsoonal and the western continental climatic zones along the 500 mm isohyet. However, this biogeographic hypothesis has been rarely tested using a phylogeographic approach. The members of the genus Allium subgenus Cyathophora coincidentally distribute across this biogeographical divide. Intriguingly, Allium fasciculatum of subgenus Amerallium co-occurs in the distribution range of subgenus Cyathophora. To illuminate the role of this biogeographic divide on the genetic divergence, we genotyped 466 individuals of 52 populations of subgenus Cyathophora and 110 individuals of 19 populations of A. fasciculatum using three chloroplast DNA fragments, whole nrITS and nine nuclear microsatellite loci, supplemented with the present environmental space and paleo-distribution modeling. Our phylogeographical evidence recovered the concordant east–west genetic breaks both for subgenus Cyathophora and A. fasciculatum along the 500 mm isohyet. The divergence time estimations and environmental niche differentiations suggested this east–west genetic breaks could have been triggered by the climatic-induced vicariance during the early Pleistocene. Noticeably, this split within subgenus Cyathophora could have been deepened by the morphological vicariance from the eastern umbel to the western spicate, while that within A. fasciculatum could have been obscured due to the pollen flows from the east to west caused by the postglacial expansion. The genetic structures and ecological niche modelings (ENMs) recovered the distinct responses to the Quaternary climatic oscillations for species constricted to different climatic zones, further highlighting the profound effect of the climatic differences and tectonic uplifts on the genetic diversification. Overall, our findings offer strong evidence for the existence of a biogeographic divide between the eastern monsoonal and the west continental climatic zones of the QTP nearly along the 500 mm isohyet.


INTRODUCTION
The Asian monsoon system, as a by-product of the Qinghai-Tibet Plateau (QTP) uprising, has become an increasing attraction not only because of its controversial onset, also its profound effect on the reorganization of the Asian climate system (Sun and Wang, 2005;Liu and Dong, 2013). With the establishment of Asian monsoon regimes around the Oligocene/Miocene boundary, the broad arid belt which widely stretched across China from west to east in the late Cretaceous and Paleocene, has been strongly restricted to the northwest China, and the arid climate conditions have been replaced by the humid ones in East China (Sun, 1979;Guo et al., 2008;Miao et al., 2012). Undoubtedly, the episodes of the dramatic tectonic uplifts in the QTP around 15-13 Ma, 8-7 Ma and 3.5-1.6 Ma, could have strengthened the Asian monsoon at different topographic units of the QTP (Liu and Dong, 2013;An et al., 2014). This accordingly resulted in the heterogeneity of climatic conditions in the QTP Wan et al., 2007;Zhang et al., 2012), as indicated by the occurrence of the 500 mm isohyet during these processes (Sun and Wang, 2005). It has been hypothesized that the 500 mm isohyet subdivides the QTP into the eastern monsoonal climatic zone (majorly controlled by the Indian Asian monsoon) and the western continental climatic zone (controlled by the mid-latitude westerlies). Consequently, the extreme continental glaciers with higher the actual equilibrium line altitudes (ELAs) could have developed in the western platform during the Pleistocene glaciations, whereas the dispersed mountain glaciers with lower ELAs in the eastern regions (Figure 1). This is because the ELAs determined by the balance of ablation and accumulation of glaciers, is higher in arid areas and lower in humid areas under the same temperature conditions (Shi et al., 1998;Shi, 2002;Liu and Dong, 2013). Thus, the ELAs can be used as a predictor to explore how species responded to the Quaternary climatic oscillations across the 500 mm isohyet.
The previous studies have borne out that the climatic-induced lines might have played an important role to plant divergence. For example, the 'Tanaka-Kaiyong Line' in Southwest China, acting as a geographic boundary of East Asian Sino-Himalayan and Sino-Japanese Floras, promoted the intraspecific divergence of Sophora davidii (Fan et al., 2013). Likewise, the broad arid belt now existing in the west-most part of China, worked as a climatic barrier to cause the north-south split within the Asian butternut (Bai et al., 2016). However, whether the biogeographic divide along 500 mm isohyet was also valid in the genetics and acted as a natural line of vicariant event to partition the plants in the QTP into west and east genealogies, has not been examined. Thus, pursing phylogeographic studies for the endemic plant components in the QTP regarding their genetic breaks along 500 mm isohyet during the regional monsoon intensification episodes, will shed light on the role of this biogeographic divide.
Among those endemic plant components in the QTP, the genus Allium L. (Amaryllidaceae) subgenus Cyathophora species extend their distribution range across the eastern monsoon and western continent climatic zones along the 500 mm isohyet. Strikingly, the clear morphology differences occur between populations restricted to these two zones: the eastern species comprising A. mairei (AMA), A. rhynchogynum, A. cyathophorum (ACY), A. farreri (AFR), and A. tetraploidum (ATE, unpublished new species, Li et al., 2016) have umbel inflorescence, while the western species only including A. spicatum (ASP) presents spicate inflorescence (Figure 1). Thus, this morphological differentiation concomitant with speciation in subgenus Cyathophora could be a perfect proxy to reflect the role of the biogeographic divide along 500 mm isohyet on genetic divergence. Interestingly, the field investigations found that A. fasciculatum (AFS) belonging to subgenus Amerallium, sympatrically distributes with A. cyathophorum in the eastern QTP, and with the A. spicatum in the western QTP (Figure 1 and Supplementary Table S1), however, no obvious morphology variations are presented between the east and west groups of A. fasciculatum (Supplementary Figure S1). In this circumstance, a comparative model can be constructed to test whether a concordant east-west genetic beak occurred both within subgenus Cyathophora and A. fasciculatum which are genetically distant. In addition, the comparisons of the genetic structures in the same/distinct climatic zone(s) can test whether the climatic differences could have an influence on the species' responses to the Quaternary climatic oscillations. Therefore, combined with the present environmental space and the ecological niche model (ENM) of species distributions, we conducted the phylogeographical analyses for subgenus Cyathophora and A. fasciculatum using three chloroplast fragments (trnD-trnT, trnL-trnF, and rps16), whole nrITS, and nine pairs of microsatellites, aiming to (i) examine the common genetic breaks along the 500 mm isohyet; (ii) test the relevance between climatic differences and the genetic breaks; (iii) address the influence of climatic differences on species' response to Quaternary climatic oscillations.  (Li et al., 2016) (Supplementary  Table S1). Allium macranthum was used as outgroup to examine the monophyly of A. fasciculatum (Supplementary Figure S2). Voucher specimens were deposited in Sichuan University Herbarium (SZ).

MATERIALS AND METHODS
Total genomic DNA was isolated from the well silica-dried foliar tissue using a protocol of plant genomic DNA kit (Tiangen FIGURE 1 | The spatial variation of the equilibrium line altitudes (ELAs) during the LGM in the QTP [referred from Shi (2002) to in honor of his contributions] and localities of sampled populations (p1-71) for subgenus Cyathophora species (p1-52) and A. fasciculatum (p53-71) across the 500 mm isohyet (blue line), which separates the QTP into the eastern monsoonal (humid-semihumid) and western continental (arid-semiarid) climatic zones [referred from Sun and Wang (2005)]. The species abbreviation holds same in the present paper: AFS, A. fasciculatum; AMA, A. mairei; ACY, A. cyathophorum; AFR, A. farreri; ATE, A. tetraploidum; ASP, A. spicatum.
Biotech, Beijing, China). Three chloroplast DNA fragments (trnD GUC -trnT GGU , trnL UAA -trnF GAA , and rps16) and the whole nrITS were amplified using the common primer pairs for the individuals of A. mairei and A. fasciculatum. One individual of Allium macranthum was amplified to confirm the monophyly of A. fasciculatum in the cpDNA-based phylogenetic analyses. For trnD GUC -trnT GGU and trnL UAA -trnF GAA , the PCR began with 94 • C for 4 min, followed by 35 cycles of 94 • C for 45 s, 60 • C for 1 min, and 72 • C for 1 min with a final-extension for 10 min at 72 • C. For rps16 and nrITS, the most conditions followed the above one, except the different denaturation temperatures, respectively, at 52 • C and 55 • C (Supplementary Table S2). PCR for all primers was conducted in a 30 µL volume comprising 3 µL plant total DNA, 0.8 µL forward primer, 0.8 µL reverse primer and 15 µL 2× Taq MasterMix (cwbio, Beijing, China). The purified PCR products of DNA sequences were bidirectionally sequenced with the primers in ABI Genetic Analyzer (Applied Biosystems Inc., Foster City, CA, United States). Using a double-suppression PCR method, nine pairs of microsatellite primers from 25 ones (Lee et al., 2011;Hyoun-Joung et al., 2012) ( Supplementary Table S2), were cross-amplified and genotyped for these 576 individuals in this study. PCR products were separated on 3.5% of agarose gel followed by staining with ethidium bromide. Alleles were sized using Peak Scanner software v.1.0 (Applied Biosystems).

Genetic Dataset Analysis
DNA sequence was edited to get consensus sequence using SeqMan (DNAstar, Burland, 2000). Sequence alignment was conducted in Clustal X version 2.0 (Larkin et al., 2007), followed by fine manual adjustment. DNAsp v 5.0 (Librado and Rozas, 2009) was used to generate the haplotype files and to estimate haplotype diversity (Hd) and nucleotide diversity (π). Genealogical relationships among haplotypes were inferred from a median-joining method optimized by maximum parsimony in Network 4.6.1.3 1 . The neutrality tests for DNA sequences using the Tajima's D test and Fu' Fs statistics were calculated by Arlequin version 3.5 (Excoffier and Lischer, 2010). The newly haplotype sequences were submitted to GenBank with KY744819-KY744939. The haplotypes of DNA sequences for A. cyathophorum, A. farreri, A. tetraploidum, and A. spicatum were cited from Li et al. (2016) with GenBank numbers KP164313-KP164411, KP114563-KP114599, KP127671-KP127673, and KT369742-KT369758. MICRO-CHECKER v. 2.2.3 (Van et al., 2004) and Cervus v. 3 (Marshall et al., 2010) were used to test for evidence of deviations from the Hardy-Weinberg equilibrium, genotyping errors and null alleles. FSTAT 2.9.3 (Goudet, 2001) was used to estimate the genetic diversity statistics of nSSRs dataset. For each microsatellite locus, the genetic diversity was estimated using the observed number of alleles (A O ), the observed heterozygosity (H O ), the expected heterozygosity (H E ), the average genetic diversity within the population (H S ) and fixation index F is . For each population, the genetic diversity was assessed across all loci using A O , H O , H S , F is , and allele richness (R S ). Linkage disequilibrium (LD) among loci was tested in FSTAT using Delta' and R 2 . Genetic groups in the nSSR dataset were identified by a Bayesian model-based clustering approach using STRUCTURE v2.3.4 (Pritchard et al., 2010), following the admixture model and assuming independent allele frequencies among populations. Ten independent runs were performed for each K (=1-20) with 100,000 MCMC replications and 10,000 used as burn-in steps. To determine the true number of gene pools indicated by the optimal K-value, we monitored the mean of log e P(D), with independent runs for each K (Pritchard et al., 2000). Alternatively, plots of ad hoc posterior probability models of Delta K were also observed to predict the most appropriate K (Evanno et al., 2005). Principal coordinate analysis (PCoA) was conducted on the microsatellite data using GENALEX 6.5 (Peakall and Smouse, 2012).
Using DNA sequences, the phylogeographic structure for each taxon was estimated in PermutCpSSR1.2.1 (Pons and Petit, 1996) by calculating interpopulation differentiation (G ST ) and the number of population subdivision (N ST ), with running 1,000 random permutations for the significance of N ST > G ST . The average genetic diversity within populations (H S ), total genetic diversity (H T ) were also computed. An analysis of molecular variance (AMOVA) was conducted using ARLEQUIN version 3.5 to quantify variation in DNA sequences (cpDNA and nrITS) and nSSRs using -and R-statistics, with 10,000 permutations to test the significance (Excoffier et al., 1992).

The Divergence Time Estimations
In the absence of known palaeo-botanical information for Allium or related genera to calibrate tree nodes, a constant mutation rate 1.52 × 10 −9 s/s/y within the cpDNA dataset (Yamane et al., 2006) was used to estimate the divergence time in BEAST v.2.5.0 (Bouckaert et al., 2014). The best-fitting substitution model (GTR+I+G) was calculated in jModelTest 2.1.10 (Posada, 2008). An uncorrelated lognormal relaxed clock and Yule process tree prior was specified with a random starting tree. A total 1 http://www.bio-soft.net/tree/Network.htm of 50,000,000 generations are set in Markov chain, both echo state to screen and log parameters every 5,000 generations. The first 10% trees were discarded as burin-in and the remaining were used to get a maximum clade credibility (MCC) tree with TreeAnnotator v1.10. Finally, MCC tree with mean ages for each node and their 95% credible intervals was displayed in Fig Tree 1 (Hijmans et al., 2005) for each georeferenced location. The determinant bioclimatic variables which exhibited pairwise Pearson correlation coefficients r < 0.7 for each taxon of interest were finally selected. Of which, Bio1 Annual mean temperature, Bio2 Mean diurnal range, Bio3 Isothermality, Bio4 Temperature seasonality, Bio12 Annual precipitation, Bio14 Precipitation of driest month, Bio15 Precipitation seasonality, Bio17 Precipitation of driest quarter and Altitude were included (Supplementary Table S5). A principle component analysis (PCA) was first performed to identify the ecological niche differentiation based on the determined bioclimatic variables using prcomp function (Dray and Dufour, 2007). The PCA loading scores for axes one (PC1) and two (PC2) of each determined bioclimatic variable were extracted using get_pca_ind function in R factoextra package. ANOVA permutation analysis (PERMANOVA) was then conducted to estimate the variation of principle components (PC1 and PC2) and each determinant bioclimatic variable within the taxa of interest using R lmPerm package to identify the bioclimatic variables which determined species' divergence (Wheeler, 2010). The principle components and bioclimatic variables with significant differences indicated by the ANOVA were further confirmed by the Tukey's honestly significant difference (HSD) test using the R Stats package.
To assess the effect of geographic distances and environmental conditions (i.e., IBD and IBE) to the observed genetic differentiation, the correlations between pairwise F ST estimated from the microsatellite markers and geographic distance were investigated, as well as between pair F ST and environmental distance. A matrix of pairwise F ST between populations was calculated using ARLEQUIN version 3.5. A matrix of geographic distance was calculated by the geographic coordinates (Supplementary Table S1) using GENALEX 6.5, and matrices of the environmental distance were computed using the significantly different bioclimatic variables, PC1 and PC2 on a basis of Euclid method. Mantel tests were conducted on the matrices of pairwise F ST /(1-F ST ) and geographical/environmental distances by running 10,000 permutations in R vegan package.

The Ecological Niche Modeling
The ecological niche modeling (ENM) was performed in MAXENT v3.3.1 (Phillips and Dudik, 2008) to predict potential distribution range shifts for involved species at the present, the Last Glacial Maximum [LGM, c. 21 kyr before present (BP)] and the Last inter-glacial (LIG, c. 130 kyr BP), respectively. These determinant bioclimatic variables exhibiting pairwise Pearson correlation coefficients r < 0.7 of each group (Supplementary Table S5) were used to perform ENM for current distributions. The ecological layers for the last glacial maximum (LGM) and Last Inter-Glacial were also retrieved from the WorldClim website at 2.5 arc-min and 30 arc-seconds spatial resolution, respectively. For LGM prediction, we used data from the Model for Interdisciplinary Research on Climate (MIROC). 75% of the distribution sites were randomly selected as a prediction program with the commission error taken into consideration, and the remaining were conducted to estimate the goodness of the prediction. The maximum number of iterations and the convergence threshold were set at 1,000 and 1 × 10 −5 , respectively. The accuracy of each model prediction was evaluated by using the area under the receiver operating characteristic curve (AUC; Fawcett, 2006).
The cpDNA-derived AMOVA revealed a stronger population structure for each species and group ( ST ≥ 0.59, all P < 0.001), except for A. spicatum with a weaker one ( ST = 0.32, P < 0.001). For the east and west groups of A. fasciculatum (east/west, ST = 0.59/0.65, both P < 0.001), 25% ( CT = 0.25, P < 0.001) of the total variation was partitioned to between groups, and 65% explained by variation among populations within group ( Table 3).

The Divergence Time Estimations for Each Lineage
The crown node of east-west split within subgenus Cyathophora was dated around 2.11 Ma (95% HPD, 1.35-3.05 Ma; Figure 2A), and that within A. fasciculatum around 1.37 Ma (95% HPD, 0.54-2.92 Ma; Figure 3A). The divergence among the clades of A. mairei was dated back to Pliocene-Pleistocene (3.15-1.34 Ma), as well as the divergence among the western clades of A. fasciculatum (4.34-1.74 Ma). The intraspecific divergence and intra-clade divergence for these species mainly were dated back to the largest Quaternary glaciation in the QTP (c. 1.2-0.6 Ma) (Yi et al., 2005). For example, the intraspecific divergence of A. cyathophorum was estimated around 1.22 Ma, that of A. spicatum around 0.97 Ma, that of A. farreri around 0.69 Ma and that of A. tetraploidum around 0.57 Ma (Figure 2A). For the intra-clade divergence of A. mairei was estimated around 1.06 to 0.53 Ma (Figure 2A), and that of the west group of A. fasciculatum around 0.93 to 0.76 Ma, and that of the east group of A. fasciculatum around 0.68 Ma ( Figure 3A).

nrITS Diversity and Network Topology
A total of 46 and 6 ribotypes based on 133 and 4 variable sites were, respectively, identified for subgenus Cyathophora species (Supplementary Table S9) and A. fasciculatum (Supplementary Table S10). The nrITS-based H T was higher than H S , but lower than the cpDNA-based H T (Table 2). However, no   well-defined phylogeographic structures were detected in each species (all P > 0.05), except in A. mairei (N ST > G ST , P < 0.05) ( Table 2). Neutrality test statistics based on Tajima's D and Fu's Fs were generally non-significant for all taxa (Supplementary Table S8).

Nuclear Microsatellite Diversity and Population Structure
A few null alleles were detected in subgenus Cyathophora and A. fasciculatum (Supplementary Table S11). When the null alleles were filtered, little effects on the genetic structures of subgenus Cyathophora and A. fasciculatum were showed. Thus all the nSSRs were used in the following analyses. In total, 251 and 2 | Estimates of the nucleotide diversity (π), average genetic diversity within populations (H S ), total genetic diversity (H T ), interpopulation differentiation (G ST ), and the number of population subdivision (N ST ) (mean ± SE in parentheses) on cpDNA and nrITS datasets for subgenus Cyathophora species and A. fasciculatum.  Table 1). No significant genotypic linkage disequilibrium for pairwise loci was detected within any population (P > 0.05). The nSSR-based AMOVA demonstrated significant population structure for each species (R ST ranged from 0.16 to 0.59, all P < 0.001). For the east and west groups of A. fasciculatum, 24% ( CT = 0.24, P < 0.05) of the total variation was partitioned between groups, and 40% explained among populations within group (Table 3). STRUCTURE results showed that the individuals of subgenus Cyathophora and A. fasciculatum can be, respectively, clustered into seven (K = 7) ( Figure 5B) and three groups (K = 3) ( Figure 5C) according to the highest lnP(D) and delta K values (Supplementary Figure S3). The geographic distribution of seven clusters of subgenus Cyathophora was highly congruent with that of DNA sequence (cpDNA and nrITS) in representing species (Figure 4), whereby individuals of A. cyathophorum (ACY1 and ACY2) and A. spicatum were further subdivided (ASP1 and ASP2) (Figures 5A,B). The geographic distribution of three clusters of A. fasciculatum was greatly consistent with that of cpDNA data in representing the east and west groups of A. fasciculatum (Figure 4A), whereby the west individuals of A. fasciculatum were further subdivided into West1 and West2 (Figures 5A,C). Observationally, the genetic admixture line near to the 500 mm isohyet comprising three western populations (ND65, LS60, and NML61) was also found in the STRUCTURE result (Figure 5C), as that presented in the nrITS data ( Figure 4B). Likewise, this genetic admixture was also supported by the PCoA for A. fasciculatum, in despite of the absence of genetic admixture in LS60 and NML61, being thought probably due to the limited individuals for each population ( Figure 5D and Table 1).

PCAs for Ecological Differentiation, IBDs and IBEs
More than ten georeferenced localities were sampled for each group of interest, with exception of A. tetraploidum having nine ones due to its narrow distributions (Supplementary Table S3). Results of the PCAs are exhibited in Figure 6. For all six species, the PCA recovered two components that cumulatively explained 77.3% of variation ( Figure 6A). The scatter plot showed when these six species were split into east and west groups, the greater niche overlap was presented along PC1 and PC2 ( Figure 6B). East group including A. mairei, A. cyathophorum, A. farreri, A. tetraplodium and the east populations of A. fasciculatum, occupies an ecological niche with higher annual mean temperature (Bio1), isothermality (Bio3) and precipitation of driest month (Bio14), while lower temperature seasonality (Bio4) and precipitation seasonality (Bio15), whereas west group comprising A. spicatum and the west populations of A. fasciculatum, occupies a contrast ecological niche compared to the east group. The niche overlap was also found in the east and west groups of subgenus Cyathophora (totally 79.1% of variation, Figure 6C) and A. fasciculatum (totally 66.6% of variation, Figure 6D) along PC1 and PC2. The scatter plot presented that the east group of subgenus Cyathophora has a distinct ecological niche with higher annual mean temperature (Bio1), isothermality (Bio3) and precipitation of driest quarter (Bio17), while the west group has an ecological niche with higher mean diurnal range (Bio2), precipitation seasonality (Bio15) and altitude ( Figure 6C). For A. fasciculatum, the scatter plot showed that the east group has an ecological niche with higher annual mean temperature (Bio1) and precipitation of driest month (Bio14), while the west group possesses an ecological niche with higher mean diurnal range (Bio2) and altitude ( Figure 6D). The large overlap of ecological niche was revealed in the sympatric species both for A. cyathophorum and the east group of A. fasciculatum (Figure 6E), as well as A. spicatum and the west group of A. fasciculatum ( Figure 6G). Additionally, the geographically parapatric species in subgenus Cyathophora (A. mairei vs. A. cyathophorum, and A. cyathophorum vs. A. farreri) exhibited slightly overlap, while the allopatric species A. mairei vs. A. farreri showed clearly ecological niche differentiation. By contrast, the allotetraploid species A. tetraploidum showed great niche overlap with its potential diploid parents A. cyathophorum and A. farreri, and a larger niche breadth compared to the diploid parental progenitors (Figure 6F).
The significant difference of determinant bioclimatic variables in the specific group was approximately consistent between the AVOVA tests and Tukey's HSD tests (Supplementary Table  S12). For the east-west split among all species, Bio1, Bio2, Bio4, Bio14, Bio15 and altitude were significantly different in both tests (P < 0.05), while PC1 was significant difference in the Tukey's HSD tests (P = 0.0137), but not in the AVOVA test (P > 0.05). For such east-west split in subgenus Cyathophora, Bio1, Bio2, Bio15, Bio17 and altitude were significantly different in both tests (P < 0.001). For such split within A. fasciculatum, Bio1, Bio12, Bio14, Bio15 and altitude were significantly different in both tests (P < 0.02). There were non-significant differences of the determinant bioclimatic variables in the sympatric species including the west group of A. fasciculatum -A. spicatum, and the east group of A. fasciculatum -A. cyathophorum (both P > 0.05) (Figure 6). The Mantel test showed no significant IBDs throughout subgenus Cyathophora (r = −0.063 P = 0.456), and A. fasciculatum (r = 0.057, P = 0.08). Likewise, no significant IBEs were detected within A. fasciculatum (all P > 0.05, Supplementary Table S13), however, the significant IBEs were identified within subgenus Cyathophora, whereby Bio2 (r = 0.132, P = 0.007) and Bio3 (r = 0.156, P = 0.012) could be the most relevant variables explaining the population genetic structures of subgenus Cyathophora (Supplementary Table S13).

The Ecological Niche Modeling
The ENM simulations based on Maxent model showed a perfect predictive power for A. fasciculatum, A. spicatum, and the eastern subgenus Cyathophora species, with the AUC = 0.967 ± 0.009, 0.986 ± 0.003 and 0.981 ± 0.005 (mean ± SD), respectively. The model predicted suitable habitats for these units were highly coincident with their actual geographic range (Figures 7A,D,G), except that of A. spicatum expanding to the eastern QTP, however, with a clear gap between the western and eastern distribution, which may attribute to these determinant climate variables resulting in their niche overlap ( Figure 6C). Compared with the present distributions, the suitable habitats of A. spicatum and the west group of A. fasciculatum were largely contracted toward the eastern Himalayas with the lower ELAs during the LGM period, whereas this tendency was weakened during the LIG period (Figures 7A-F). The suitable habitats of the east group of A. fasciculatum, A. farreri, and A. tetraploidum, respectively, were largely contracted in situ during the LGM, however, no significant changes for that of A. mairei and A. cyathophorum (Figures 7A,B,G,H).

Climatic-Induced Phytogeographic Divide Along the 500 mm Isohyet in the QTP
The morphology clearly demonstrates that there are two distinctive groups for the genus Allium subgenus Cyathophora, which well correspond to the eastern monsoonal and the western continental climatic zones in the QTP across the 500 mm isohyet. Of which, the eastern group with umbel inflorescence comprises A. mairei, A. cyathophorum, A. farreri and A. tetraploidum, while the western group with spicate inflorescence only includes A. spicatum (Figure 1). However, A. fasciculatum which co-exists with A. cyathophorum in the eastern QTP, and with A. spicatum in the western QTP, shows no obvious morphological variation (Supplementary Figure S1). Strikingly, the genetic breaks both in subgenus Cyathophora and in A. fasciculatum strongly supported this climatically east-west split across the 500 mm isohyet (Figures 4, 5). Nonetheless, this split in the nuclear genealogy of A. fasciculatum has been obscured most likely due to the genetic introgression indicated by the narrow genetically admix cline (Barton and Gale, 1993), which geographically was adjoining to the 500 mm isohyet (Figures 4B, 5A). The significant phylogeographical structures in cpDNA data but not in nrITS data for most species strongly supported the limited seed flows but frequent pollen flows among these Allium populations in the QTP ( Table 2). In combination with the expansion of the derived R47 ( Figure 3C) from the east to the west within A. fasciculatum, it is suggested the pollen flows from the east group to the west group in A. fasciculatum could have resulted in this narrow genetic admixture line (Figure 3), which reasonably explains why the genetic admixture only was detected in the nuclear genealogy but not in the cpDNA genealogy (Rieseberg and Soltis, 1991).
The divergence time estimations for this east-west break were traced back to the early Pleistocene (c. 2.11 Ma for subgenus Cyathophora; c. 1.37 Ma for A. fasciculatum) (Figures 2, 3), which predates the largest glacial age in the QTP (c. 1.2-0.6 Ma), but coincidently falls into the climatic condition differentiation in the QTP around 3.6-1.2 Ma (Liu and Dong, 2013), suggesting the east-west genetic breaks within subgenus Cyathophora and A. fasciculatum could have been caused by the climatic differences. This hypothesis is also supported by the distinct climatic niches between the east and west groups (Figure 6B) determined by the significantly different bioclimatic variables (Supplementary Table S12) both in subgenus Cyathophora and A. fasciculatum (Figure 6). Noticeably, the larger niche overlap between the east and west groups (Figure 6), probably suggests the shallow or intermediate climatic-derived divergence (Antonelli, 2017). However, no significant correlations between the genetic divergence and bioclimatic variables were found within A. fasciculatum, but the significantly ones within subgenus Cyathophora (Supplementary Table S13). Considering the striking morphological differentiations between the east (umbel inflorescence, larger seeds) and west (spicate inflorescence, smaller seeds) groups of subgenus Cyathophora (Figure 1) (Li et al., 2016), it is assumed this climatic-induced genetic beak could have been deepened in subgenus Cyathophora due to the prezygotic sterility resulted from the morphological vicariance (e.g., Landmann and Winding, 1993;Kennedy et al., 2012;Kou et al., 2014;Wen et al., 2014). However, no such case occurs between the east and west groups of A. fasciculatum (Supplementary Figure S1), it is asserted that the frequent pollen flows within A. fasciculatum could have obscured this climatic-induced genetic break (Figures 4B, 5A).
To sum up, it seems reasonable to support that the climatic-derived 500 mm isohyet acts as a natural line of vicariant event to partition subgenus Cyathophora and A. fasciculatum into western and eastern genealogies during the climatic difference episode.

Divergence Caused by Tectonic Uplift and Quaternary Climatic Changes
Our phylogenetic analyses indicated that A. fasciculatum could have ancestrally located in the western QTP (Figure 3), and subgenus Cyathophora ancestrally settled in the eastern QTP (Figure 2) (Li et al., 2016). In terms of 'center-periphery' hypothesis, it is predicted a higher genetic diversity, lower genetic differentiation and more constant population size in the ancestral areas, but not in the recolonized areas (Eckert et al., 2008). The DNA sequence-based AMOVA recovered a higher genetic differentiation ( ST ) and total genetic diversity (H T ) in the ancestral areas than those in the recolonized areas both for subgenus Cyathophora (eastern ST ≥ 0.77, western ST = 0.32 for cpDNA; eastern ST ≥ 0.39, western ST = 0.14 for nrITS) and A. fasciculatum (western ST = 0.65, eastern ST = 0.59 for cpDNA; western ST = 0.89, eastern ST = 0.65 for nrITS) (Tables 2, 3). Based on the cpDNA sequences, a higher average genetic diversity within populations (H S = 0.479) and a larger effect population size inferred from the nucleotide diversity (π = 0.93 × 10 −3 ) were found in the western ancestral group of A. fasciculatum, when comparing to the eastern derived group (H S = 0.287 and π = 0.14 × 10 −3 ) ( Table 2). Conversely, it was revealed that all the eastern species of subgenus Cyathophora have lower average genetic diversity within populations (H S ≤ 0.381) and smaller effect population size (π ≤ 0.24 × 10 −3 ) when comparing to the western species A. spicatum (H S = 0.541 and π = 0.50 × 10 −3 ) ( Table 2). The BEAST-based divergence estimations showed that the crown divergence at the ancestral areas both for subgenus Cyathophora (5.19-3.26 Ma; Figure 2A) and A. fasciculatum (4.34-1.74 Ma; Figure 3A) was estimated to occur around the Pliocene-Pleistocene, coinciding to the tectonic motions of the QTP . It is thus inferred that the tectonic uplift could have facilitated greatly for the genetic differentiation among populations in the ancestral areas, and the ancestral effect population of subgenus Cyathophora could have been subdivided into different small parts. This tectonic uplift-driven diversification has been increasingly recognized as an important process that can have diverse outcomes with respect to the biodiversity in the QTP (e.g., Favre et al., 2015;Zhang et al., 2016;Xing and Ree, 2017).
The nSSRs-based AMOVA revealed a lower genetic differentiation both for subgenus Cyathophora (eastern R ST ≤ 0.43, western R ST = 0.49) and A. fasciculatum (western R ST ≤ 0.59, eastern R ST = 0.11). Compared to cpDNA and nrITS, microsatellite markers usually have faster mutation rates and thus are more related to more recent population history. The intraspecific divergence within subgenus Cyathophora (1.06-0.69 Ma; Figure 2A) and the intra-clade divergence of A. fasciculatum (0.69-0.33 Ma; Figure 3A) were mainly dated back to the largest Quaternary glaciation in the QTP (c. 1.2-0.6 Ma), probably highlighting the importance of Quaternary glaciations to their recent diversification in the QTP Wen et al., 2014;Du et al., 2017).
Of the east species, the ENMs showed larger contractions of the suitable habitats, respectively, for A. farreri, A. tetraploidum, and A. fasciculatum during the LGM (Figure 7), suggesting a contraction-expansion pattern for these species when they responded the Quaternary glacial oscillations. The single dominant haplotype among populations, respectively, within A. farreri (C19/R24), A. tetraploidum (C1/R1) and A. fasciculatum (C56/R47) (Figures 2, 3) most probably indicates the species specific glacial refugium, and the population diversifications could have been resulted in by this contraction-expansion pattern. For A. mairei, no clear contraction of suitable habitat during the LGM (Figure 7), the greatly dominant population-specific haplotypes (Figure 2) and significantly phylogeographical structure (Table 2), suggest that it could have survived from the dispersed mountain glacial refugia  resulted from tectonic uplift in its early divergence (Figure 2A). For A. cyathophorum, two genetic lineages (Figures 2, 5B) possibly indicate two isolated glacial refugia (Figure 2). However, the striking lack of phylogeographic structure within A. cyathophorum (G ST > N ST , P > 0.05) ( Table 2), possibly suggests that the frequent gene flows could have blurred the imprints of refugial isolation. However, this is challenged by the range displacement most likely due to the invasion of A. tetraploidum to A. cyathophorum (the genetic information from A. tetraploidum was detected in p9 of A. cyathophorum, Li et al., 2016). The STRUCTURE-based result showed that the p7-11 of the ACY1, were uniquely sub-clustered, but not for the p12-16 of ACY1 and p17-21 of ACY2, which contained genetic information from both lineages, presumably indicating the formation of a moving hybrid zone due to the genetic admixture from ACY1 to ACY2 ( Figure 5A) (Godbout et al., 2012;Cinget et al., 2015). The wider niche for A. tetraploidum than the parental progenitors ( Figure 6F), meaning the stronger ability to inhabit different conditions, and to use different resources (Carrillo-Angeles et al., 2016), probably provides a prerequisite for this invasion. In contrast, a significantly phylogeographical structure was recovered in the east population of A. fasciculatum (Table 2), nevertheless the sympatric distribution and similar niches to A. cyathophorum ( Figure 6E). Considering their approximately same divergence time (Figures 2A, 3A), it is suggested that the range displacement could have a deep influence on the genetic structure of A. cyathophorum.
Of the west species, the ENMs presented the same contraction-expansion pattern, respectively, for the suitable habitats of A. spicatum and A. fasciculatum when responding FIGURE 6 | Scatter plots of PC1 and PC2 (A-D) showing ecological differentiation between the east and west groups, respectively, for all species, subgenus Cyathophora and A. fasciculatum, as well as (E-G) the niche differentiation among all east species and all west species based on the bioclimatic variables exhibiting pairwise Pearson correlation coefficients r < 0.7 (Supplementary Table S5). * * P < 0.001, * P < 0.05 indicating significant difference of the bioclimatic variables between the east and west groups estimated by PERMANOVA and Tukey's HSD tests (Supplementary Table S12). the glacial crisis (Figure 7), which could have contributed much to the divergence of A. spicatum (Figure 2A) and the intra-clade divergence of A. fasciculatum ( Figure 3A). However, the phylogeographic structures ( Table 2) and AMOVA results indicate a distinct population genetic structures between A. fasciculatum and A. spicatum (Table 3), notwithstanding their similar niches (Figure 6G). It deems this pattern can be much attributed to the key role of the tectonic motions during the Pliocene-Pleistocene in the early divergence of A. fasciculatum ( Figure 3A). This proposal is also supported by the similar genetic structures between A. mairei and the west group of A. fasciculatum (Table 3), most likely due to the tectonic motions for their early inter-clade divergence (Figures 2A, 3A). The contracted ranges of A. spicatum and A. fasciculatum during the LGM indicate their refugia would locate at eastern Himalayas with lower ELAs (Figure 7). This strongly supports the local survival of A. spicatum and A. fasciculatum from the Quaternary freezing damages in the western QTP, agreeing with some tolerant cold organisms, such as Aconitum gymnandrum (Wang et al., 2009), Juniperus tibetica complex (Opgenoorth et al., 2010), Hippophae thibetana , Orinus thoroldii (Liu et al., 2015) and Eospalax baileyi (Tang et al., 2010). It is implied that the aridification in the western QTP possibly could have reduced the prevalence of large ice sheets, characterized . The dots signifying the species' present georeferenced localities, and the dash line dividing the eastern species of subgenus Cyathophora into the northern species A. tetraploidum and A. farreri, and the southern species A. cyathophorum and A. mairei.
by the higher ELAs (Figure 1) (Shi et al., 1998;Shi, 2002), which in turn provides a promise for the local survival from the glacial crisis.
Taken together, it is inferred that the eastern species have the species specific glacial refugia and consequently resulted in the diversifications in situ, most probably due to the development of the dispersed mountain glaciations. By contrast, the western species mostly could have largely contracted to the eastern Himalayas where the larger ice sheet could have limited due to the drought in the western QTP. In this way, the distinct responses to the Quaternary climatic oscillations could have been recovered for species restricted to different climatic zones. Additionally, the early divergence trigged by the tectonic motions and genetic introgression could have a deep effect on the species' genetic structures.

CONCLUSION
Our chloroplast and nuclear DNA results for subgenus Cyathophora and A. fasciculatum, along with the present environmental space, robustly support the existence of a biogeographic divide between the eastern and western QTP along the 500 mm isohyet. The causative factor for this genetic divide is suggested to be the climatic differentiation in the eastern and western QTP. Interestingly, this climate-mediated split within subgenus Cyathophora could have been deepened due to the prezygotic sterility resulted from the morphological vicariance from the east umbel inflorescence to the west spicate inflorescence, while that within A. fasciculatum could have been obscured due to the pollen flows from the east to west caused by the postglacial expansion. It is recovered that the tectonic motions could be responsible for the ancestral divergence, and the Quaternary climatic fluctuations for the recent divergence. Strikingly, distinct responses to the Quaternary climatic fluctuations have been revealed in the components of the eastern and western QTP, which could be closely related with the climatic differences. Overall, our findings highlight the importance of orogeny and climatic changes in the QTP in shaping species diversification during the Pliocene-Pleistocene.