Diverse range dynamics and dispersal routes of plants on the Tibetan Plateau during the late Quaternary

Phylogeographical studies have suggested that several plant species on the Tibetan Plateau (TP) underwent recolonization during the Quaternary and may have had distinct range dynamics in response to the last glacial. To further test this hypothesis and locate the possible historical dispersal routes, we selected 20 plant species from different parts of the TP and modeled their geographical distributions over four time periods using species distribution models (SDMs). Furthermore, we applied the least-cost path method together with SDMs and shared haplotypes to estimate their historical dispersal corridors. We identified three general scenarios of species distribution change during the late Quaternary: the ‘contraction-expansion’ scenario for species in the northeastern TP, the ‘expansion-contraction’ scenario for species in the southeast and the ‘stable’ scenario for widespread species. During the Quaternary, we identified that these species were likely to recolonize along the low-elevation valleys, huge mountain ranges and flat plateau platform (e.g. the Yarlung Zangbo Valley and the Himalaya). We inferred that Quaternary cyclic glaciations along with the various topographic and climatic conditions of the TP could have resulted in the diverse patterns of range shift and dispersal of Tibetan plant species. Finally, we believe that this study would provide valuable insights for the conservation of alpine species under future climate change.


Introduction
Phylogeography has provided a way to understand the demographic history of species, which was primarily affected by periodic climate oscillations during the Quaternary, generating varying degrees of range contraction and expansion [1]. In general, species retreated to refugia during the glacial period, preserving multiple and original haplotypes within these regions [2]. During the postglacial period, several species recolonized in other areas, producing nearly complete haplotype uniformity in the colonized populations [3]. Hence, the range dynamics of

Species selection and determination of shared haplotypes
To test whether species in different parts of the TP have distinct range dynamics, we chose 20 plant species from three regions of the TP (Table 1), which experienced demographic expansion during the Quaternary. Based on the phylogeographical study of each species, we collected the haplotype composition and geographical coordinate of each population. To ensure the congruence in DNA markers, we mainly selected maternally inherited markers (e.g. cpDNA for angiosperms and mtDNA for gymnosperms). From the spatial haplotype structures of each plant (S1 Fig), we obtained which populations had shared haplotypes.

Modelling species distribution
We applied the available bioclimatic and occurrence data to model the species distribution. Based on the WorldClim database (http://www.worldclim.org/) [49], we obtained 19 bioclimatic variables (S1 Table) Table) were mainly collected from published literatures and further supplemented with data from field surveys and online database, including the Global Biodiversity Information Facility (GBIF, http://data.gbif.org) and the Chinese Virtual Herbarium (CVH, www.cvh.org.cn). As spatial clusters of localities in the gathered occurrence data may result in model's inflation and over-fit towards environmental biases [51][52][53], we applied a spatially rarefy method to filter the occurrence data [53]. Based on the climatic heterogeneity, localities of occurrence were spatially filtered at 5, 15 and 25 km 2 in areas of high, medium and low climatic heterogeneity, respectively. To calculate climatic heterogeneity, we conducted a principal component (PC) analysis of 19 bioclimatic variables and then calculated the mean standard deviation of the first three climate PCs. These analyses were implemented using the SDMtoolbox v1.1c [50]. Finally, the remaining localities combined with low-correlated bioclimatic layers were used to estimate species distribution.
In this study, we used the maximum entropy algorithm in MAXENT v3.3 [9] to model the current and past distributions of 20 species with a convergence threshold of 10 −5 , 500 iterations, and the localities of occurrence divided into testing data sets and training data sets (25% and 75%, respectively). We then conducted simple statistics on data distribution of range dynamics for each species and different regions in R 3.2.5 [54].

Visualizing the dispersal routes
We identified the dispersal routes of plant species on the TP under the least-cost path (LCP) method using the SDMtoolbox v1.1c [50] in ArcGIS 10.1 (ESRI). (i) We first obtained a dispersal cost layer (resistance layer) by inverting the species distribution layer (1-SDM) and then created a cost distance raster for each sample locality based on the resistance layer. (ii) Using the cost distance raster, we obtained the corridor layers between two localities with shared haplotype. (iii) In contrast to previous studies [12,55], which generally obtained linear dispersal routes, we applied the categorical LCP method to better represent environmental heterogeneity in dispersal. We classified the value of each corridor layer into four intervals (three cutoff values); if the lowest value of corridor layer was hypothesized to be 1.00, the four intervals were: 1.00~1.01 (1.00 + 1.00 Ã 1%), 1.01~1.02 (1.00 + 1.00 Ã 2%), 1.02~1.05 (1.00 + 1.00 Ã 5%) and 1.05~highest value. These four intervals were then reclassified as new values (5, 2, 1, 0, respectively). (iv) We summed all of the pairwise reclassified corridor layers, and the spatial dispersal routes for each plant were obtained. Finally, based on the routes of 20 species, we identified the possible recolonization pathways of species on the TP during the Quaternary.

Results
During the late Quaternary, plant species in three regions of the TP had distinct range dynamics (Fig 2, S2 Fig). For the species in the northeastern part, they mainly presented contractionexpansion scenario before and after the LGM. In contrast, the species in the southeastern TP underwent surprising range expansion during the LGM and then contracted to some extent. For the widespread species across the TP, their distribution exhibited almost no change since the LIG. Phylogeographical studies have suggested that at least 20 plant species on the TP experienced recolonization during the Quaternary. To date, we have identified the dispersal routes of each plant (S3 Fig) and predicted several major recolonization paths for species on the TP (Fig 3). These dispersal routes are also supported by fossil data and other evidence (summarized in Table 2). We found that at least 3 dispersal corridors existed, and the Hengduan Mountains (HM) seemed to act as the main dispersal source (i.e. glaciation refugia). These paths included: i) the northward path from the northern part of the HM to the Qilian Mountains; ii) the westward path, in which species started in the HM region and recolonized along the Yarlung Zangbo River or Himalaya Mountains; and iii) the southward path, which went from the northeastern part to the southern part throughout the majority of TP's platform, converged with the Yarlung Zangbo River, and moved further into the vast western region of the TP.

Range dynamics of plants on the TP during the late Quaternary
Within the European and North American continents, contraction-expansion in species distribution was generally predominant [66,67], owing to the existence of a large and complete ice sheet during the LGM. However, there was no such formation of an ice sheet on the TP [21], and our analysis suggests that the range dynamics of species on the TP were complex during the late Quaternary. It might be argued that these distinct shifts in distribution were associated with the heterogeneous climate and topography across the TP.
In the northeastern TP, the tabula rasa scenario seemed to be dominant, which is consistent with previous phylogeographical studies and pollen fossil records [25,26,64]. In this region, the climate is colder and drier than that of the southern TP, and the platform is basically flat; only the eastern part has a series of mountain ranges (e.g. the Qilian Mountains and the Minshan Mountains in the northern HM). During the LGM, according to the rare palaeovegetation on the platform and the climatic characteristics [68,69], we assumed that platform populations of species became extinct and only marginal populations survived within the eastern mountains. Subsequently, when the warmer and wetter Holocene period came, forest and alpine meadow recolonized on the platform [64]. Therefore, the alternation of vegetation on the platform suggests that the LGM might have had greater effects on species distribution in this region.
In the southeastern TP, the species exhibited an expansion-contraction response during the late Quaternary. We assumed that this response was closely related to regional topographic and climatic features. This region is mainly located within the East Himalaya and Hengduan Mountains, which provide ample altitudinal range (e.g. niche) and stable ecological environment [70]. On the other hand, the Indian monsoon produces a warmer and wetter climate compared to other regions of the QTP [31,71]. The pollen data indicate that the late Quaternary was characterized by a warm temperate climate [72], and phylogeographical studies suggest that this region acted as the major glaciation refugia for many alpine species on the TP [18,73]. Hence, during the LGM, suitable habitats could have facilitated the persistence of Phylogeographic analysis [11,60] The great canyon of Yarlung Zangbo River has abundant plants, which is a young floristic region originated from the East Asia flora since the Tertiary; Floristic analysis [61] The species Pedicularis longiflora and Tsuga dumosa underwent westward dispersal along this path; Least-cost path analysis [16,62] The seed dispersal of Sophora moorcroftiana relied on the wind and waterflow in the midstream of Yarlung Zangbo River.
Seed dispersal experiment [63] Northward path (including the Hehuang Valley and the Qilian Mountains) The alternation of grassland and forest in the northeastern TP during the late Quaternary; Vegetation reconstruction based on pollen [64,65] The species Juniperus przewalskii and Picea crassifolia underwent invasion of platform after the last glaciation.
Phylogeographic analysis [25,26] Southward path The species Potentilla fruticosa and Allium przewalskianum underwent extensive recolonization events from the northeastern to the southern TP. populations in low-elevation areas, and populations on the top of mountains could have moved downward and expanded into other areas because of the abundant niche in this region. In contrast, during the LIG and Holocene period, the warmer and wetter climate would have driven alpine species upward, resulting in their limited distribution in an isolated area. For the widespread species on the TP, as discussed above, during the late Quaternary, the populations in the northeastern platform may have contracted to refugia, whereas the populations in the southeastern part may have expanded in range. For the populations in the central and southern TP's platform, a lot of evidence suggests that the LGM had a limited influence on species distribution of the platform [22,56]. In addition, there might be several microrefugia (e.g. the Three Rivers' Headstream Regions) within these platform areas [74,75]. Therefore, as a whole, widely distributed species maintained stable populations throughout the LGM.

Dispersal routes of plants on the TP
Phylogeographical studies have been able to test which species experienced rapid and extensive range expansion. In this study, using the LCP method, in combination with SDMs and shared haplotypes, we have identified the possible dispersal routes of Tibetan plant species during the Quaternary. Undoubtedly, determining the historical dispersal routes is important to supplement the demographic history of species on the TP.
In this study, we identified three common pathways across the TP. Among the initial locations of these routes, the East Himalaya-Hengduan Mountains seemed to be the dispersal source, which may be associated with this region's role as the major glaciation refugia [18,73]. According to historical biogeographic studies, this place was also an important center of origin and radiation of alpine plants, triggering numerous species to migrate into other areas from this region in the Neogene or earlier period [34,58]. Surprisingly, the migration routes estimated in those studies partially overlapped with the paths during the Quaternary. Thus, it might be argued that species on the TP tend to disperse along several constant routes. We assumed that this could be related to the topographic and monsoon climatic features of the TP. The northward and westward routes identified in this study were basically along the Hehuang Valley and Yarlung Zangbo Valley. Why did species migrate along these valleys? These valleys are not only lower elevation areas (average elevation < 4000 m), but also important vapor channels; in particular, the Yarlung Zangbo Valley is the largest vapor channel of the TP [30]. Hence, with lower topography and strong airflow, these valleys are likely to facilitate the dispersal of seeds and pollen. Additionally, many fossils found in these valleys along with phylogenetic and floristic analysis also support the existence of these corridors (see Table 2). Apart from these two valleys, we predicted that the Yangtze River, Mekong River and Salween River might also be important dispersal pathways. Even not obviously identified in our study, other evidence implied that alpine species would migrate along these valleys during the glacial period [10,11].
Mountain ridges are another important dispersal corridors, like the Himalaya Mountains in this study, which is a very long-narrow mountainous area. Evidence has shown that the Quaternary glaciations likely had a great impact on plants in this region, resulting in many species becoming extinct [60]. However, because the Himalaya is linked with the Hengduan Mountains, some species could have easily moved to the Himalaya during the postglacial period. In other montane areas around the world, huge mountains (e.g. the Andes) are also regarded as important dispersal corridors [76].
We found that the TP's platform (average elevation > 4500 m), which has flatter topography, could facilitate long distance dispersal of alpine species, like the southward path in this study. Despite the existence of a series of east-westward mountains across the platform, their relative elevation is actually not so high. Through Barrier analysis, we have also demonstrated that these mountains are not acted as geographical isolations (unpublished data). Thus, several drought-tolerant species (e.g. Potentilla fruticosa and Allium przewalskianum), originating in the northern TP or adjacent desert (e.g. the Tengger Desert), could migrate into the southern Tibet across the entire TP [22,23].

Conclusions
The plant species on the TP exhibited diverse demographic histories during the late Quaternary. There were at least three types of changes in range. Despite the detailed range shifts are distinct across different parts of the TP, they actually had no significant difference, and as a whole, they remained relatively stable status in response to the last glacial, suggesting that LGM had a limited impact on species distribution. In fact, an earlier glaciation (the Naynayxungla glaciation, 0.72-0.50 Ma) might have had a significant influence. After this glaciation, species were continuously dispersed along the low-elevation valleys (e.g. the Yarlung Zangbo Valley and Hehuang Valley), huge mountain ranges (e.g. the Himalaya) or even the flat plateau platform. Under future climate change, the species on the TP could continue to move along these pathways, and thereby these areas should be paid more attention. Besides, the inferences in this study are just based on 20 species, we believe that more diverse range dynamics and dispersal corridors would be found if many species are involved in future studies.
Supporting information S1   Table 1). Pie chart corresponds to the sample locations of each species and the diverse colors from green to red are attributed to different chlorotypes or mitotypes.