Climatic Factors Drive Population Divergence and Demography: Insights Based on the Phylogeography of a Riparian Plant Species Endemic to the Hengduan Mountains and Adjacent Regions

Quaternary climatic factors have played a significant role in population divergence and demography. Here we investigated the phylogeography of Osteomeles schwerinae, a dominant riparian plant species of the hot/warm-dry river valleys of the Hengduan Mountains (HDM), Qinling Mountains (QLM) and Yunnan-Guizhou Plateau (YGP). Three chloroplast DNA (cpDNA) regions (trnD-trnT, psbD-trnT, petL-psbE), one single copy nuclear gene (glyceraldehyde 3-phosphate dehydrogenase; G3pdh), and climatic data during the Last Interglacial (LIG; c. 120–140 ka), Last Glacial Maximum (LGM; c. 21 ka), and Current (c. 1950–2000) periods were used in this study. Six cpDNA haplotypes and 15 nuclear DNA (nDNA) haplotypes were identified in the 40 populations of O. schwerinae. Spatial Analysis of Molecular Variance, median-joining networks, and Bayesian phylogenetic trees based on the cpDNA and nDNA datasets, all suggested population divergence between the QLM and HDM-YGP regions. Our climatic analysis identified significant heterogeneity of the climatic factors in the QLM and HDM-YGP regions during the aforementioned three periods. The divergence times based on cpDNA and nDNA haplotypes were estimated to be 466.4–159.4 ka and 315.8–160.3 ka, respectively, which coincide with the time of the weakening of the Asian monsoons in these regions. In addition, unimodal pairwise mismatch distribution curves, expansion times, and Ecological Niche Modeling suggested a history of population expansion (rather than contraction) during the last glaciation. Interestingly, the expansion times were found being well consistent with the intensification of the Asian monsoons during this period. We inferred that the divergence between the two main lineages is probably caused by disruption of more continuous distribution because of weakening of monsoons/less precipitation, whilst subsequent intensification of the Asian monsoons during the last glaciation facilitated the expansion of O. schwerinae populations.

Hot/warm-dry river valleys are one of the most prominent natural landscape features in the HDM and adjacent regions [15]. In comparison to studies of (sub)alpine plants, however, few riparian plants of hot/warm-dry river valleys in these regions have been studied. The phylogeographic studies of riparian plants conducted so far have focused on two species: Terminalia franchetii (Combretaceae) and Buddleja crispa (Buddlejaceae), both of which indicated that the discontinuous distribution patterns of these plants are strongly linked to historical drainage re- organization events [16,17]. More phylogeographic studies of riparian plants in the HDM and adjacent regions (particularly the QLM) are required to obtain a better understanding of their recent population history; and in particular, to detect whether some riparian species in these regions have been affected by climatic factors (i.e. heterogeneity and changes) during the Quaternary. River valleys in different regions probably differ in climate because of different monsoon systems. The QLM, which is located on the north-eastern margin of the HDM (Fig 1), is dominated by the East Asian monsoon, originating from the Pacific Ocean [18,19]. The east-west orientated mountains form a natural obstacle, preventing the wet summer East Asian monsoon from penetrating northern China, but provide good corridors for it to penetrate into the western QLM or even HDM. In contrast, The HDM is affected by both the southwest monsoon from the Indian Ocean and the East Asian monsoon from the Pacific Ocean, and the YGP is mainly affected by the East Asian monsoon [20]. Similarly, the south-north orientated river valleys in the HDM could also provide convenient corridors for the warm, moist air from the Indian and Pacific Oceans to penetrate into the northern and/or north-eastern HDM [21,22]. During the Quaternary, the Asian monsoons appear to have fluctuated with the alternating interglacial-glacial periods, with more intense summer monsoons during interglacial periods and weaker monsoons during glaciations [23][24][25]. However, some strong summer monsoons seem to have also occurred during glacial periods (e.g. the last glaciation) in addition to their anticipated occurrence during interglacial periods. Conversely, some significantly weak summer monsoons have also been found to have occurred during interglacial periods (e.g. the last interglacial period) in addition to their occurrence during glacial periods [26][27][28]. Given the river valleys' role as corridors for monsoons, it seems plausible that such climatic heterogeneity and changes due to Asian monsoons would have an impact on the population divergence and demography of riparian plants.
Osteomeles schwerinae Schneid. (Rosaceae) is an iconic riparian plant species endemic to the dry river valleys of the HDM, YGP, and QLM [29][30][31] (Fig 1), and probably ideal to study such impacts. It is a small tree species, characterized by imparipinnate leaves with 5 to 10 pairs of elliptic or elliptic-oblong leaflet blades, a campanulate hypanthium possessing sepals, and stony pomes [32]. According to the Flora of China [32] and our field observations, it commonly grows on the slopes of dry river valleys. It includes two varieties (O. schwerinae var. schwerinae and O. schwerinae var. microphylla Rehder & E. H. Wilson). The former is mainly confined to the river valleys of the HDM-YGP and has larger leaflet blades (5-10 mm), while the latter is mostly present in the river valleys of the QLM and has smaller leaflet blades (3-5 mm) [32]. Although there is now much molecular evidence, as mentioned above, that the population divergence and demographic history of taxa (e.g. alpine plants) in the HDM and adjacent regions were closely linked with historical climatic events [9,11,20], very little climatic data have been explicitly used to examine these inferences. In order to understand better how climatic factors contributed to the present-day distribution and population divergence of plants in the HDM and adjacent regions, we constructed the phylogeography of O. schwerinae based on molecular approaches and paleo-and current climatic data. We had the following objectives: (1) to quantify whether there is genetic divergence between the populations of O. schwerinae; (2) to explore whether the genetic divergence coincide with the regional morphological differentiation of O. schwerinae in the QLM and HDM-YGP regions; (3) to date divergences and determine whether these dates relate to geological events (e.g. historical drainage re-organization) or climatic changes (e.g. the monsoon fluctuation and/or the glacial-interglacial alternation); and (4) to unravel the demographic history of O. schwerinae.

Ethics Statement
Osteomeles schwerinae is not protected nor endangered in the sampled areas, and all samples were collected by researchers following current Chinese regulations. None of the sampled locations are privately owned or protected by law.

Sample Collection, DNA Extraction, and Sequencing
We sampled 304 individuals from 40 populations of O. schwerinae covering its geographical range; out of these, eight populations originated from the QLM region, and 32 populations from the HDM-YGP region (Table 1). Individuals were randomly sampled and were at least 100 m apart within each population. Fresh green leaves were collected from each individual and dried in silica gel for DNA extraction.

Molecular Data Analysis
Haplotypes of the combined three cpDNA sequences were produced in DnaSP v5.0 [36]. Method of Clark (1990) was used to identify the haplotypes for G3pdh heterozygotes. In this method, sequences of heterozygotes were compared with those of homozygotes until we could account for the observed combination of double-banded sites. ArcMap v9.3 (ESRI, Inc.) was used to plot the distribution of haplotypes on a relief map coving the range. Original map in the figures was downloaded from http://www.diva-gis.org/Data. The elevation dataset in the figures was derived from NASA's SRTM data (http://www2.jpl.nasa.gov/srtm). Spatial Analysis of Molecular Variance (SAMOVA) was conducted to analyze the spatial genetic structure of O. schwerinae using SAMOVA v1.0 [37]. This program uses a simulated annealing approach to define K groups of populations that are geographically homogenous and maximally differentiated from each other. In this analysis, we assumed that K varied from 2 to 13, and computed a corresponding F CT index, which measures the proportion of genetic differentiation among the K groups [38]. The configuration with the F CT value that achieves a plateau first is usually considered the optimal grouping of populations. An analysis of molecular variance (AMOVA) was conducted within ARLEQUIN v3.5 [39] to partition variation within and among defined groups and populations.
Genealogical relationships among haplotypes were inferred from an unrooted statistical parsimony network within Network v4.6.1.2 [40] (available at http://www.fluxus-engineering. com) using a median-joining method. Indels (gaps) were treated separately as binary states, namely single events, following Caicedo and Schaal (2004) [41]. Mononucleotide repeats and duplicated indels of sequences in this study were removed as they introduced homoplasious reticulations into the network, due to their high levels of bidirectional mutation [42]. Observed and inferred haplotypes were nested following the rules in Posada and Crandall (2001) [43].
Phylogenetic relationships among haplotypes were reconstructed by Bayesian inference. This phylogenetic analysis was conducted in MrBayes v3.1.2 [44] using substitution models (GTR+G) explored by Modeltest v3.7 [45]. Four independent Markov chain Monte Carlo analyses with runs of 10 million generations were carried out, sampling at every 1000 generations. The first 20% of generations was discarded as burn-in, while the remaining data was used to construct a 50% majority rule consensus tree. The posterior probabilities found in this consensus tree were given to evaluate the robustness of the trees. The divergence time of haplotypes was estimated via a molecular dating method, and estimated under a strict molecular clock in BEAST v1.7.5 [46], using the substitution rate range for angiosperm species (cpDNA: 1.01-3.0×10 −9 substitutions per site per year (s/s/y) [47,48]; nDNA: 1.58-3.15×10 −8 s/s/y [48]). Results from BEAST were compiled and visualized in the program Tracer v1.5 [49].
Pairwise mismatch distribution was carried out using ARLEQUIN to infer the demographic history of O. schwerinae. In this analysis, an expected distribution was generated from 1000 parametric bootstrap replicates to test 'demographic expansion' and 'spatial expansion' model, respectively [2]. Generally, a unimodal mismatch distribution curve indicates a recent population expansion, whereas a multimodal distribution usually suggests that populations are at demographic equilibrium. The populations for which the hypothesis of expansion was not rejected, the value for the mode of the mismatch distribution (Tau) was assessed by 1000 parametric bootstrap replicates in ARLEQUIN, and was converted into estimates of time since expansion (t, in years BP) using t = Tau/2μkg [2]. In this formula, μ is the substitution rate, and was assumed to be based on the substitution rate range for cpDNA mentioned above; k is the average sequence length used for analysis (2039 bp in this study); and g is the generation time in years, and was assumed to be 6 years here based on personal observations of on the age of first reproduction of O. schwerinae in cultivation at Kunming Botanical Garden.

Climatic Data Analysis
Climatic data [50] that summarized temperature and precipitation information (S1 Table), were downloaded from the WorldClim database (http://www.worldclim.org/), and were used to identify climatic factors potentially associated with the genetic structure and demographic history identified in this study. DIVA-GIS v7.5 [51] was used to extract the 19 climatic factors during the LIG, LGM, and Current periods (S2 Table). The mean values and standard errors of these factors for populations of the QLM and HDM-YGP regions were calculated with SPSS v18.0.0 [52], and were further analyzed using independent-Samples t-tests to validate the impact of the climatic factors on the genetic differentiation between the QLM and HDM-YGP regions.
Ecological Niche Modeling (ENM) can simulate past species distributions, and further elucidate demographic history through time [53]. Therefore, ENM was also performed in this study to validate potential present and past climate envelopes for this species. Assuming the species has not changed with respect to its climatic preference, we re-constructed the range of O. schwerinae during the LIG, LGM, and Current periods. Species occurrence data for O.  Table) and sampling records collected in this study (Table 1). ENM based on the 19 BioClim factors and species occurrence data was performed using MAXENT v3.3.3k [54] with the default parameters for number of iterations (500) and convergence threshold (10 −5 ). We set up 10 replicate runs in each analysis in order to ensure more reliable results. Meanwhile, in order to test the performance of each model, 25% of the data in each run was randomly chosen by MAXENT and compared with the model output created with the remaining data. To reduce the effects of spatial autocorrelation, duplicate records from the same locality were also removed. The area under the receiver operating characteristic curves (AUC) was used to compare model performance. An AUC value of 0.5 indicates that the performance of the model is no better than random, while values closer to 1.0 indicate better model performance [54].

Chloroplast and Nuclear DNA Sequencing Data
Variations in the length of two sequenced regions of cpDNA were detected (953-958 bp, trnT-trnD; 658-697 bp, psbD-trnT), while the petL-psbE was uniformly 449 bp long in our samples. Removal of mononucleotide repeats and duplicated indels for the alignment of three cpDNA regions from the 304 samples resulted in an overall alignment length of 2039 bp, and produced six cpDNA haplotypes (designated H1-6; S4 Table). The aligned sequences of the G3pdh gene were 550 bp long, and produced 15 nDNA haplotypes (designated Hap1-15) based on 18 nucleotide substitutions (S5 Table). These sequences were deposited in the GenBank database (accession numbers KT724331-KT724335, KT750890-KT750920).

Genetic Divergence and Demography
In the SAMOVA analysis based on cpDNA sequences, two groups corresponding to the QLM and HDM-YGP regions were well defined when K = 2. Although the F CT value was not the highest, the F CT values reached a plateau, and fluctuated little with an increasing number of groups (S1 Fig). In addition, when divided into two groups, the AMOVA revealed high genetic variation partitioned between the groups (68.57%), and low variations among populations within groups (18.60%) as well as within populations. As shown in Fig 2A, 8 populations located in the QLM region were contained in one group, whilst the remaining 32 populations from the HDM-YGP were included in the other group. The SAMOVA analysis based on nDNA sequences showed relatively low F CT values, and fluctuated little with the increase of K (S2 Fig). Nonetheless, when divided into two groups, the AMOVA revealed that still 30.55% of the total nDNA variation was due to differences among the two groups identified, comparing with 22.61% among populations and 46.84% within populations. In this scenarios of grouping, 8 populations from the QLM (pops. 1-8) and 2 populations from the HDM-YGP (pops. 9 and 11 of the Dadu River) formed one group, whereas the remaining 30 populations from the HDM-YGP were included in the other group (Fig 2B). Genealogical analysis of cpDNA haplotypes identified a common haplotype (H4) shared by populations along all rivers surveyed in QLM, HDM, and YGP regions (Fig 3A). Apart from this common haplotype, the remaining haplotypes were distinctly regional, and could be divided into two clusters from the QLM and HDM-YGP regions (Fig 3A). In addition, it is worth noting that the intermediate zone between these two regions (namely the drainage basin of the Dadu River) possess high cpDNA haplotype diversity, containing the common haplotype (H4), the main regional haplotypes from the QLM (H1) and HDM-YGP (H2), as well as its own unique haplotype (H5) ( Table 1 ; Fig 2A). Genealogical analysis of nDNA haplotypes showed similar structure, which also identified a common haplotype (Hap3) and two clusters of regional haplotypes from the QLM and HDM-YGP, respectively (Fig 3B). Simultaneously, the intermediate zone between QLM and HDM-YGP (the Dadu River basin) also showed high nDNA haplotype diversity, possessing the common haplotype (Hap3), regional haplotypes both from the QLM (Hap2, 5, 11) and HDM-YGP (Hap7), as well as its own unique haplotype (Hap12; Table 1; Fig 2B).
The Bayesian phylogenetic tree clustered five of the cpDNA haplotypes (all except H4) into two lineages, respectively from the QLM and HDM-YGP regions (Fig 4A). In the phylogenetic tree based on nDNA haplotypes, one haplotype (Hap15) and the dominated haplotype (Hap10) in HDM-YGP region fell in one lineage, while three haplotypes (Hap 1, 4, 14) from QLM in another lineage (Fig 4B). In addition, this tree showed several polytomies, however, it is not surprising since the nuclear gene is two parents-inherited. Based on molecular dating, the divergence times of the cpDNA and nDNA haplotypes were estimated to be 466.4-159.4 ka and 315.8-160.3 ka, respectively.
The mismatch distribution curves for the total, QLM and HDM-YGP populations were clearly unimodal (Fig 5). Therefore, the hypothesis of population expansion could not be rejected. Based on the Tau value, the demographic and spatial expansion times for the total, QLM and HDM-YGP populations fell within the last glacial period (Table 2).

Climatic Data
Out of the 19 bioclimatic factors, all but five (bio2, 5, 8, 10, 15) during the LGM and the current periods showed highly significant differences in mean values between populations of the QLM and HDM-YGP regions (Table 3). When comparing with the LGM and Current periods, two more bioclimatic factors (bio8 and 15) during the LIG exhibited highly significant differences in mean values between the two regions ( Table 3). As shown in Table 4, there was reduction in temperature during the LGM. For instance, the annual average temperature (bio1) in the QLM and HDM-YGP regions decreased 2.3 and 2.6°C during the LGM than the LIG, respectively. In contrast, the majority of the bioclimatic factors relating to the precipitation (6/8 = 75%) for both the QLM and HDM-YGP during the LGM period were higher than during the LIG, with most variables (5/6 = 83.33%) being significantly higher (Table 4). It suggested a more moist environment in the QLM and HDM-YGP regions during the LGM, whereas a more drought environment in the LIG.
The AUCs of the ENM for the total, QLM and HDM-YGP populations under all climate scenarios were 0.987, indicating quite better model performance of ENM. The current  schwerinae. Paleo-distribution modeling for the total, QLM and HDM-YGP populations, all indicated a much more restricted distribution during the LIG than the Current, with subsequent expansion during the LGM. This predicted subsequent expansion covered a slightly greater range than that predicted based on the current climatic factors (Fig 6).

Genetic and Morphological Divergence Relative to Climatic Heterogeneity and Changes
In this study, the SAMOVA, median-joining networks, and Bayesian phylogenetic trees, all suggested population genetic divergence between the QLM and HDM-YGP regions (Figs 2, 3 and 4), which coincide with the morphological divergence of the two varieties of O. schwerinae in the QLM and HDM-YGP regions. Climatically, the QLM has been primarily dominated by the East Asian monsoon, and is situated near the northern limit of this monsoon system [19] (Fig 1). However, the HDM-YGP has been affected both by the East Asian monsoon from the Pacific and the southwest monsoon from the Indian Ocean [20] (Fig 1). In addition, the HDM-YGP has experienced more intense monsoons than the QLM owing to its more southerly latitude. Therefore, the climate in the HDM-YGP is generally more warm and moist than that in the QLM. Our climate data further confirmed this, demonstrating that the temperature and precipitation values in the QLM were significantly lower than that in the HDM-YGP during the LIG, LGM, and Current periods (Table 3). Considering the two varieties of O. schwerinae, the leaflet blades of the one growing in the QLM (3-5 mm) are smaller than the one (5-10 mm) growing in the HDM-YGP. Since the sizes and shapes of leaves are quite sensitive to moisture, with leaves in dry climates generally smaller than in wet climates [55], it is quite likely that the difference in morphology of O. schwerinae between the two regions is the result of long-term adaptation to the significant climatic heterogeneity between the two regions.
It has been speculated that drainage re-organization was a main cause of genetic divergence for riparian plants in the HDM and adjacent regions [16,17], but we did not detect its effects on O. schwerinae. When using the same substitution rate range for cpDNA (1.01-3.0×10 −9 s/s/ y) as previous studies on riparian plants [16,17] in these regions, the genetic divergence time in this study was estimated to be 466.4-159.4 ka (Fig 4A). In addition, the divergence time based on nDNA was calculated to be 315.8-160.3 ka (Fig 4B). They are much more recent than the timings found in previous riparian plants' studies (4.24-0.4 Ma, Terminalia franchetii [17]; 3.683-1.228 Ma, Buddleja crispa [16]) and the geological drainage re-organization events (no later than the later Pliocene; 5.332-2.588 Ma) [56][57][58]. However, our divergence times largely coincide with the time of the weakening of the monsoons in these regions. According to paleoclimatic analysis, the southwest monsoon was sharply weakened from 342-118 ka [59]. Coincidentally, the East Asian monsoon was also found to decrease during the period 200-130 ka [23,    , and, as expected, displayed two more significantly different bioclimatic factors between the QLM and HDM-YGP than the LGM and the Current periods (Table 3). Finally, these two reasons together led to the genetic divergence between the QLM and HDM-YGP regions.

Demographic History of Osteomeles schwerinae
Usually, the demographic scenarios reported for species from the HDM and adjacent regions suggest distribution contraction during glacial periods, and expansions during interglacial periods [8][9][10][11]. However, our results revealed an unusual demographic history of O. schwerinae populations: the opposite pattern to most demographic scenarios in the HDM and adjacent regions. Before the period "342-118/200-130 ka" (including the LIG; c. 140-120 ka), the O. schwerinae populations probably had a more continuous or wider distribution linking the two geographic regions because of stronger monsoons/more precipitation in both regions. However, during the LIG, the O. schwerinae populations probably had a less continuous or more narrow distribution because of weaker monsoons/less precipitation. The ENM analysis further confirmed this, showing more narrow distribution during the LIG, but larger distribution ranges during the LGM period (Fig 6), which implied the "contraction during interglacial periods, and expansions during glacial periods" demographic pattern of O. schwerinae populations. Unimodal pairwise mismatch distribution curves of O. schwerinae (Fig 5) were identified in this study, which fit the sudden expansion model particularly well. Moreover, the expansion times for the total, QLM and HDM-YGP populations of O. schwerinae fell into the last glacial period (Table 2), which further supported this pattern. This pattern is similar to the "phalanx" model [61], which has only been sporadically found in a few cold-adapted species at high-elevation in (sub)tropical areas [62][63][64][65]. These plants were found to expand to lower elevations during glacial periods, allowing isolated patches to join together, because of the retreat of the (sub)tropical forests occupying lower altitudes. In contrast, these plants would be forced back to 'refugia' at higher elevations during interglacial periods as a result of the re-expansion of the (sub)tropical forests at low elevations [66]. In the HDM and adjacent regions, the majority of plants have been found to retreat to refugia because of cooling during the glacial periods, which probably left a large space for O. schwerinae populations to expand during the last glacial period. Meanwhile, because of the Foehn effect [67,68], the climate in river valleys is warmer than in neighboring areas [15]. Therefore, the cooling during the last glacial period might have had a limited impact on the populations of O. schwerinae. For instance, the annual global-mean continental cooling during the LGM has been found to have been 4.3-9.8°C [69,70]. In contrast, the climatic data for our sampling sites reveals a lower annual average temperature reduction during the LGM (only 2.3 and 2.6°C lower than the LIG for the QLM and HDM-YGP, respectively; Table 4), which probably reflects the warm climate in river valleys. In addition, owing to the intensification of the Asian monsoons affecting precipitation during the last glacial period in these regions [23,27,28,59], the precipitation at our sampling sites was significantly increased during the LGM compared to the LIG (Table 4). This may have relieved or even eliminated the impact of the drought since 342-118/ 200-130 ka [23,27,59,60].
In this context, the river valleys could have provided warm and moist habits for O. schwerinae to expand during the last glacial period, resulting in a second contact between the populations from the QLM and HDM-YGP. This second contact probably just explains the high cpDNA and nDNA haplotype diversity within their intersecting distributions in the Dadu River basin, where there is the common haplotype (cpDNA: H4; nDNA: Hap3), regional haplotypes both from the QLM (cpDNA: H1; nDNA: Hap2, 5, 11) and HDM-YGP (cpDNA: H2; nDNA: Hap7), as well as its own exclusive haplotypes (cpDNA: H5; nDNA: Hap12) ( Table 1 ;  Fig 2A and 2B).

Conclusions
Our analyses of O. schwerinae, based on molecular phylogeography and climatic data, reveal a major population divergence between the QLM and HDM-YGP regions, and its correlation to the climatic heterogeneity and changes within both regions since the late Pleistocene. We inferred that weakening of Asian monsoons since 342-118/200-130 ka increased the climatic differences between the two regions, caused the loss/reduction of gene flow, and triggered the genetic divergence between them. Subsequent intensification of the monsoons during the last glaciation facilitated the expansion of O. schwerinae populations, contributed to the second contact between the populations from the QLM and HDM-YGP, and finally led to the development of the current distribution pattern of O. schwerinae populations. This study provides a new example in which riparian plant species in the HDM and adjacent regions are shown to be affected by climatic heterogeneity and changes since the late Pleistocene. However, more studies of riparian plants in these regions are needed to explore this hypothesis.