Introduction {#s1}
============

The evolutionary history of plant species could be dramatically influenced by habitat fragmentation [@pone.0037196-Young1], [@pone.0037196-Fahrig1]. The process of fragmentation would strongly reduce the dispersal ability of individuals among resource patches, leading to the breakup of a large, genetically variable population into isolated small subpopulations [@pone.0037196-Templeton1]. An increasing number of studies indicated that the rapid decline of population size and restricted gene flow could result in genetic differentiation between populations and genetic diversity loss, giving rise to species extinction and biodiversity decrease [@pone.0037196-Young1], [@pone.0037196-Fahrig1], [@pone.0037196-Lowe1]. On the other hand, both evolutionary theory and empirical data indicate that global genetic variation could be maintained or even increased by a fragmented population, and allopatric speciation is prone to occur due to genetic isolation of fragmented populations [@pone.0037196-Carson1], [@pone.0037196-Xu1]. However, the genetic and ecological consequences of habitat fragmentation depend critically on the levels of gene flow between habitat fragments [@pone.0037196-Fahrig1], [@pone.0037196-Templeton1], [@pone.0037196-Lowe1]. So the analysis of genetic divergence and gene flow between fragmented populations would be greatly helpful in understanding the process of speciation or extinction under the effects of environmental changes.

The phased uplift of the Qinghai-Tibetan Plateau (QTP) continuing from the middle Tertiary has greatly changed the geology and topography of East Asia, leading to the unique geomorphological configuration and complex land conditions [@pone.0037196-Wang1]. The highest plateau and mountains ranges have been incised deeply by numerous valleys or rivers ([Fig. 1](#pone-0037196-g001){ref-type="fig"}), resulting in the rapid development of endemic, specialized montane species along mountain ranges and the concentration of the relatively old Tertiary flora in isolated low valleys [@pone.0037196-Chang1], [@pone.0037196-Qian1]. It would be of particular interest to unravel the roles of geographic and ecological changes in the species diversification by tracing the genetic footprint of environmental shifts in the QTP plants. Phylogeography, the analysis of the spatio-temporal pattern of population genetic variation, has been proved to be very efficient in retrieving the evolutionary history of species or close relatives on a relatively short evolutionary timescale [@pone.0037196-Hewitt1], [@pone.0037196-Hickerson1].

![A map showing the species\' ranges and a scenario of possible hybridization between species.\
Letters A--D indicate the sampling sites for *M. simplicifolia* (star) and *M. betonicifolia* (filled circle), A: Cuona, B: Linzhi, C: Milin, and D: Laojunshan Mountain. Triangle, solid circle, and pentagram represent some sampled populations of *M. integrifolia*, *M. betonicifolia*, and *M. simplicifolia*, respectively. Putative hybrid populations are indicated by the overlapping of solid circle and triangle.](pone.0037196.g001){#pone-0037196-g001}

To date, the geographic patterns of population genetic variation were investigated on more than ten plant species from different lineages in the QTP [@pone.0037196-Qiu1]. The results clearly indicated the effect of the Quaternary climatic oscillations on the development of geographic genetic structure [@pone.0037196-Qiu1], [@pone.0037196-Yang1], [@pone.0037196-Cun1], and the recolonization patterns from multiple refugia on the QTP were also outlined. At the same time, the discovery of some putative refugia or micro-refugia on the plateau implied that some of the populations might have survived the Last Glacial Maximum or even earlier glaciations in situ [@pone.0037196-Yang1], [@pone.0037196-Cun1]. Also, the time of lineage divergence within some species could be dated to the late Tertiary [@pone.0037196-Wang2], [@pone.0037196-Wang3], [@pone.0037196-Zhang1], when the dramatic uplift of the QTP occurred [@pone.0037196-Wang1]. However, compared with the well-discussed population evolutionary dynamics driven by the Quaternary climatic oscillations [@pone.0037196-Qiu1], much less attention was paid to the effects of tectonic events on plant divergence and speciation in the QTP [@pone.0037196-Xu1], [@pone.0037196-Cun1], [@pone.0037196-Wang4].

The genus *Meconopsis* comprises about 50 species that are mainly distributed in the QTP and the neighboring mountains, with the only European species *M. cambrica* inhabiting in the humid and shady deciduous forest from Ireland, south-west England to Northern Spain [@pone.0037196-Wu1]. *Meconopsis integrifolia*, one of the most widely distributed species of the genus, is a flagship species of the alpine scree in the QTP. It is scattered in the scree along mountain ridges, with an island-like distribution at altitudes of 4000--5200 m above sea level. In the eastern edge of its range, the species often grows in patchy habitat under azaleas at lower altitudes of 2700--4000 m above sea level. The bright yellow flower of this species makes it easy to be distinguished from the blue (or sky blue) flower relatives *M. betonicifolia* and *M. simplicifolia*, two quite narrowly distributed species in the southeast QTP ([Fig. 1](#pone-0037196-g001){ref-type="fig"}). A range-wide phylogeographical analysis of *M. integriflolia* and its close relatives is very helpful to reveal the underlying mechanism of the habitat fragmentation and to explore the genetic and ecological consequences of environmental changes for the QTP plants.

In the present study, the cpDNA *trn*L-*trn*F and *trn*fM-*trn*S regions were used to detect genetic variation and gene flow through seed between the highly fragmented populations of *M. integrifolia* and to retrieve the phylogeographical history of the species. Since the phylogenetic analysis on the two cpDNA markers indicated that one clade of *M. integrifolia* was nested within *M. betonicifolia*, it is interesting to know whether hybridization or ancestral polymorphism is responsible for it. Different from the uniparentally inherited cpDNA markers, the nuclear genes are biparentally inherited, with a rich source of phylogenetic information, so we tried to isolate some single/low copy nuclear genes to investigate phylogenetic relationships between the species and its relatives. However, complex polyploidization is quite prevalent in *Meconopsis* [@pone.0037196-Ying1] and it is difficult to determine orthology/paralogy of the genes. The internal transcribed spacer (ITS) of nuclear ribosomal DNA (nrDNA) is tandemly repeated with hundreds to thousands of copies, which could be homogenized by concerted evolution that eliminates sequence variation among the different copies [@pone.0037196-Baldwin1]. Therefore, the nrDNA ITS region was also used to retrieve the phylogenetic relationship of *M. integrifolia* and its close relatives.

Materials and Methods {#s2}
=====================

Ethics statement {#s2a}
----------------

No special permits were required for this study because all samplings were collected by researchers with introduction letters of IBCAS (Institute of Botany, Chinese Academy of Sciences) in Beijing.

Population sampling {#s2b}
-------------------

Population sampling was conducted throughout the range of *Meconopsis integrifolia* during the summers of 2005--2009. Fresh leaves were collected from 35 populations and, with few exceptions, 15--29 individuals were sampled from each population ([Table S1](#pone.0037196.s001){ref-type="supplementary-material"}). To avoid sampling closely related individuals, the samples were at least 30 m apart from each other. In total, 757 individuals of *M. integrifolia* were sampled and leaves were dried with silica gel. In addition, fifteen individuals of *M. betonicifolia* were sampled from Cuona and Milin Counties of Tibet, and the Laojunshan Mountain in Yunnan Province. These three populations can represent the easternmost, middle and westernmost distribution of the species, respectively ([Fig. 1](#pone-0037196-g001){ref-type="fig"}). Ten individuals of *M. simplicifolia* were sampled from Linzhi and Cuona Counties in Tibet, representing its easternmost and Himalayan populations, respectively. Also, four individuals from two populations of *M. chelidonifolia* were sampled as outgroups. The latitude, longitude and altitude of each collection location were measured using an eTrex Global Positioning System (Garmin).

DNA extraction, PCR amplification, cloning and sequencing {#s2c}
---------------------------------------------------------

Total genomic DNA was extracted from silica gel-dried leaves using the modified CTAB method [@pone.0037196-Rogers1] and used as template in the polymerase chain reaction (PCR). Two cpDNA intergenic spacers, *trn*fM-*trn*S [@pone.0037196-Shaw1] and *trn*L-*trn*F [@pone.0037196-Taberlet1], were amplified and directly sequenced for all samples. Amplification was conducted in a Tpersonal thermocycler or Biometra T1 themocycler (Biometra, Goettingen, Germany). The PCR products were purified using Gel Band Purification Kit (Tiangen), and then sequenced using the ABI Prism Bigdye Terminator Cycle Sequencing Ready Reaction Kit on an ABI PRISM 3730xl analyzer. To investigate the evolutionary history of the *M. integrifolia* clade that was nested within *M. betonicifolia* in the cpDNA phylogeny, 26 individuals of *M. simplicifolia*, *M. betonicifolia* and *M. cheilanthifolia* and 151 samples from 20 populations that can represent the seven geographic groups ([Table S1](#pone.0037196.s001){ref-type="supplementary-material"}) of *M. integrifloia* were further used in the nrDNA ITS analysis. The nrDNA ITS region was amplified and sequenced with the primers ITS1 and ITS4 of White *et al*. [@pone.0037196-White1]. However, one or more double-peaks were detected in the chromatograms of 17 individuals. All these individuals, and four additional ones without double-peaks in the chromatograms, were cloned with the pGEM-T Easy Vector System II (Promega), and 8--12 clones per individual were sequenced using the primer T7 [@pone.0037196-Wallace1]. The results showed that each double-peak represents a polymorphic site, and one to three different sequence types occur in each individual.

Data analyses {#s2d}
-------------

The DNA sequences were aligned with Clustal X [@pone.0037196-Thompson1] and manually adjusted in BioEdit version 7.0.9 [@pone.0037196-Hall1]. Segregating sites and haplotypes were extracted using DnaSP v5 [@pone.0037196-Librado1]. All haplotypes have been deposited in GenBank under accession numbers (JQ798299-JQ798389). To estimate the genetic differentiation within species, analysis of molecular variance (AMOVA) was performed based on pairwise differences and haplotype frequencies with Arlequin version 3.1 [@pone.0037196-Excoffier1]. Measures of DNA divergence between groups (*F* ~ST~) [@pone.0037196-Excoffier1] were calculated with Arlequin version 3.1, and the significance was tested using 10 000 permutations. To estimate differentiation between populations, we calculated the values of *G* ~ST~ and *N* ~ST~ using program HAPLONST [@pone.0037196-Pons1], with *G* ~ST~ taking into account haplotype frequency and *N* ~ST~ taking into account both frequency and similarity of the haplotypes. Besides, the *U*-statistic was performed to compare the *G* ~ST~ and *N* ~ST~, with a higher *N* ~ST~ than *G* ~ST~ indicating the presence of phylogeographic structure [@pone.0037196-Pons1]. Phylogenetic relationships among haplotypes were reconstructed by maximum parsimony and Bayesian inference analyses. For the maximum parsimony analysis that was conducted with PAUP version 4.0b10 [@pone.0037196-Swofford1], all characters were treated as unordered and equally weighted, with gaps treated as missing. A heuristic search was performed with 1 000 random addition sequence replicates, tree-bisection-reconnection branch swapping and the MULTREES option. To examine the robustness of clades in the most parsimonious trees, the bootstrap analysis was carried out with 1 000 replicates using the heuristic search options as described above. The Bayesian analysis was conducted with MrBayes version 3.1.2 (nst = 6; rates = invgamma; ngen = 1 000 000) [@pone.0037196-Ronquist1], using the best evolutionary model determined by the Akaike Information Criteria implemented in MrModeltest 2.3 [@pone.0037196-Nylander1]. To estimate the significance of difference between the topologies of chloroplast and nuclear gene phylogenies, a Shimodaira--Hasegawa test [@pone.0037196-Shimodaira1] was performed with PAUP version 4.0b10, using 10,000 bootstrap replicates with the resampling estimated log-likelihood method. The maximum likelihood tree generated from the cpDNA data was compared with the ITS data, and the maximum likelihood tree of the ITS data was assessed relative to the cpDNA data.

To examine the constancy of molecular evolution rate among lineages of the phylogenies, a likelihood ratio test [@pone.0037196-Felsenstein1] was performed with PAUP version 4.0b10, in which likelihood scores of the trees with and without an enforced molecular clock were compared. Significance was calculated by comparing two times the difference in log likelihoods to a χ^2^ distribution with *n* -- 2 degrees of freedom, where *n* is the number of taxa (haplotypes of *M. integrifolia* plus outgroups). When *M. chelidonifolia* served as outgroup, the results of the likelihood ratio tests rejected a clock-like evolution for both cpDNA (P\<0.01) and ITS haplotypes (P\<0.01). If *M. betonicifolia* and *M. simplicifolia* were chosen as outgroups, the results of the likelihood ratio tests rejected a clock-like evolution for cpDNA (P\<0.01) but not for ITS haplotypes (P = 0.096). As no fossil of *Meconopsis* and its close relatives was available, the crown group age (a mean of 9.2 Ma; 95% highest posterior density (HPD): 15.2--3.8 Ma) [@pone.0037196-Valtuena1] of *M. integrifolia* and its close relatives was used to date lineage ages of *M. integrifolia* in the cpDNA and nrDNA ITS phylogenies. Given the wide confidence interval of this constraint age, the average nrITS substitution rate for herbaceous plants (1.72--8.34×10^−9^ substitution/site/year, mean = 4.13×10^−9^ subs/site/yr) [@pone.0037196-Kay1] was also used to estimate the crown group age of *M. integrifolia* and its close relatives. Molecular clock analysis of the age and confidence interval was performed using BEAST v1.6.2 [@pone.0037196-Drummond1], with a Hasegawa-Kishino-Yano plus Gamma site heterogeneity model and a constant population size for the data set. To cover the nrITS substitution rate range (1.72--8.34×10^−9^ subs/site/yr) with a 95% confidence interval for the mean rate, the substitution rate was set to a normally distributed prior with a mean of 4.13×10^−9^ subs/site/yr and a standard deviation of 1.47×10^−9^ subs/site/yr. The Bayesian Markov chain Monte Carlo simulation was run for 20,000,000 generations with a sample frequency of 1000, and the first 2,000,000 generations were discarded as 'burnin'. Trace file and statistics of estimated parameters were created using Tracer version 1.5 [@pone.0037196-Rambaut1]. The final calibrated chronogram and node estimates were edited using FigTree version 1.3.1 [@pone.0037196-Rambaut2]. The most recent common ancestor of *M. integrifolia* and its close relatives was dated to 6.20 Ma (95% HPD: 15.23--2.78 Ma) (not shown). Then, a relaxed clock model [@pone.0037196-Drummond2] in BEAST was used to estimate nodal ages and credibility intervals for both cpDNA and nrDNA data sets using *M. chelidonifolia* as outgroup. An uncorrelated lognormal distribution model with a Yule speciation process and a Hasegawa-Kishino-Yano plus Gamma site heterogeneity model was applied in the analyses. To cover the time range (9.2--6.2 Ma) with a 95% confidence interval for the mean age, the crown group age of *M. integrifolia* and its relatives was set to a normally distributed prior with an average of 7.7 Ma and a standard deviation of 0.9 Ma. The remaining settings and procedures were the same as in the nrITS data set described earlier.

Results {#s3}
=======

CpDNA sequence analysis {#s3a}
-----------------------

Forty-five single nucleotide polymorphisms were detected among the 757 individuals of *M. integrifolia*, which defined 26 cpDNA haplotypes (A--Z) ([Table S1](#pone.0037196.s001){ref-type="supplementary-material"}). For phylogenetic relationships of the haplotypes, the topologies of trees generated by Bayesian ([Fig. 2](#pone-0037196-g002){ref-type="fig"}), BEAST ([Fig. 3a](#pone-0037196-g003){ref-type="fig"}), and maximum parsimony analyses (not shown) were congruent in the major clades. That is, the 26 haplotypes were clustered into six well supported clades ([Fig. 3a](#pone-0037196-g003){ref-type="fig"}), and five of which (Clades II--VI) further formed a strongly supported major clade, whereas clade I is sister to a clade of *M. betonicifolia*.

![Distribution of the cpDNA haplotypes.\
Phylogenetic relationships of the haplotypes based on the Bayesian analysis are indicated on the left of the map. Pie charts show the proportions of haplotypes in each population. Doted area indicates the region controlled by the Indian summer monsoon [@pone.0037196-Wang5].](pone.0037196.g002){#pone-0037196-g002}

![Phylogenetic chronograms of the cpDNA haplotypes (a) and the ITS haplotypes (b) generated from BEAST.\
Numbers on the branches indicate the Bayesian posterior probabilities (Left) and bootstrap values of 1000 replicates (Right) for the maximum parsimony analysis. Ages of the main clades are shown below the branches and horizontal bars represent estimates for node ages at the 95% highest posterior density.](pone.0037196.g003){#pone-0037196-g003}

The 26 cpDNA haplotypes are distributed in seven non-overlapped geographic regions ([Fig. 2](#pone-0037196-g002){ref-type="fig"}). The four haplotypes (A, B, C and D) of clade I are exclusively distributed in the population group WYN, and haplotypes E, F and G (clade II) only occur in the two populations (group LZ) located in the north Yarlung Zangbo Grand Canyon. Clade III comprises four haplotypes (H, I, J and K) that only occur in the eight populations (group SLL) along the Shalulishan and Bayankalashan Mountains. Clade IV includes five haplotypes (O, P, L, M and N) that are present exclusively in five populations (group ET) in the central and eastern Tibetan. Haplotypes of clade VI are mainly distributed in the easternmost regions of the Hengduan Mountains (Group EQa & EQb), but the distribution is separated into two parts by the Daxue mountain range and Qionglai mountain range, where the population group DX (Clade V) is located. Although the geographic distribution of haplotypes generally exhibits a lineage-specific pattern ([Fig. 2](#pone-0037196-g002){ref-type="fig"}), there are a couple of exceptions. That is, haplotypes Q (clade V) and U (VIa) are only present in population 18 (group SLL, corresponding to clade III), with a low frequency of 7%, respectively ([Table S1](#pone.0037196.s001){ref-type="supplementary-material"}). Also, haplotype T (clade V) occurs in population 19 (group SLL, corresponding to clade III) with a frequency of 10%, while it is fixed at 100% in populations 27 and 28 of group DX.

Of the studied 35 populations of *M. integrifolia*, twenty-two (62.86%) harbored pure haplotype, ten (28.57%) fixed two haplotypes, and only three (8.57%) had more than two haplotypes. Significant population differentiation was observed, with a *G* ~ST~ = 0.843 and a *N* ~ST~ = 0.887. The *U*-statistic showed that *G* ~ST~ and *N* ~ST~ were significantly different from each other (*N* ~ST~\>*G* ~ST~, *U* = 2.73, *P*\<0.005). The analysis of molecular variance (AMOVA) revealed that 80.18% of the total variation occurred among groups, 16.82% among populations, and 3.00% within populations ([Table 1](#pone-0037196-t001){ref-type="table"}). Pairwise *F* ~ST~ among geographic regions varied from 0.7397 to 0.9054, suggesting a high genetic differentiation between the geographic regions ([Table 2](#pone-0037196-t002){ref-type="table"}).

10.1371/journal.pone.0037196.t001

###### Hierarchical analysis of molecular variance (AMOVA) for populations of *Meconopsis integrifolia* based on chloroplast DNA haplotypes.

![](pone.0037196.t001){#pone-0037196-t001-1}

  Source of variation    d\. f.   Sum of squares   Variance components   Percentage of variation (%)
  --------------------- -------- ---------------- --------------------- -----------------------------
  Among groups             6         1798.661             2.742                     80.18
  Among populations        28        350.282              0.575                     16.82
  Within populations      722         74.098              0.102                     3.00
  Total                   756        2223.041             3.419         

10.1371/journal.pone.0037196.t002

###### Pairwise comparisons of *F* ~ST~ (chloroplast DNA haplotypes) between seven geographic regions of *Meconopsis integrifolia*.

![](pone.0037196.t002){#pone-0037196-t002-2}

                           ET                                        LZ                                        WYN                                       EQa                                       EQb                                       SLL                     DX
  ----- ----------------------------------------- ----------------------------------------- ----------------------------------------- ----------------------------------------- ----------------------------------------- ----------------------------------------- ----
  ET                                                                                                                                                                                                                                                                
  LZ     0.7397[\*](#nt101){ref-type="table-fn"}                                                                                                                                                                                                                    
  WYN    0.8463[\*](#nt101){ref-type="table-fn"}   0.7485[\*](#nt101){ref-type="table-fn"}                                                                                                                                                                          
  EQa    0.8632[\*](#nt101){ref-type="table-fn"}   0.7554[\*](#nt101){ref-type="table-fn"}   0.8655[\*](#nt101){ref-type="table-fn"}                                                                                                                                
  EQb    0.8787[\*](#nt101){ref-type="table-fn"}   0.8455[\*](#nt101){ref-type="table-fn"}   0.9054[\*](#nt101){ref-type="table-fn"}   0.7255[\*](#nt101){ref-type="table-fn"}                                                                                      
  SLL    0.7487[\*](#nt101){ref-type="table-fn"}   0.7809[\*](#nt101){ref-type="table-fn"}   0.8799[\*](#nt101){ref-type="table-fn"}   0.8390[\*](#nt101){ref-type="table-fn"}   0.8673[\*](#nt101){ref-type="table-fn"}                                            
  DX     0.8095[\*](#nt101){ref-type="table-fn"}   0.7932[\*](#nt101){ref-type="table-fn"}   0.8788[\*](#nt101){ref-type="table-fn"}   0.8889[\*](#nt101){ref-type="table-fn"}   0.8984[\*](#nt101){ref-type="table-fn"}   0.8116[\*](#nt101){ref-type="table-fn"}  

*P*\<0.001.

NrDNA ITS sequence analysis {#s3b}
---------------------------

Twenty-five nrDNA ITS haplotypes (alleles) were detected from 151 individuals representing the seven population groups ([Table S1](#pone.0037196.s001){ref-type="supplementary-material"}). Phylogenetic relationships of the haplotypes are consistent between the Bayesian tree (not shown), maximum parsimony tree (not shown), and the tree generated by the BEAST analysis ([Fig. 3b](#pone-0037196-g003){ref-type="fig"}). Monophyly of the 25 haplotypes is supported by the high bootstrap value and posterior probability ([Fig. 3b](#pone-0037196-g003){ref-type="fig"}), and four major clades (I~nu~, II~nu~, IV~nu~ and VII~nu~) with high bootstrap and posterior probability values were recognized. In the ITS tree ([Fig. 3b](#pone-0037196-g003){ref-type="fig"}), the group WYN (Clade I~nu~) is a natural group of *M. integrifolia* rather than being nested into *M. betonicifolia* as in the cpDNA tree ([Fig. 3a](#pone-0037196-g003){ref-type="fig"}). A significant topological incongruence between the two trees was supported by the Shimodaira-Hasegawa test (P\<0.001). Similar to the distribution of the cpDNA haplotypes ([Fig. 3a](#pone-0037196-g003){ref-type="fig"}), the clade I~nu~ and II~nu~ of the ITS haplotypes are exclusively restricted to population groups WNY and LZ (corresponding to clades I and II in the cpDNA phylogeny), respectively ([Table S1](#pone.0037196.s001){ref-type="supplementary-material"}). Clade IV~nu~ includes seven haplotypes (H8--H14) that only occur in population group ET (corresponding to clade IV in the cpDNA phylogeny), but two haplotypes (H18, H19) of clade VII~nu~ occur in population 2 of this group. The monophyletic relationships of groups SLL, DX and EQ were not supported by the nrDNA ITS analysis.

Molecular dating {#s3c}
----------------

When the crown group age of *M. integrifolia* and its close relatives was set to 7.7(±1.5) Ma as a constraint, the most recent common ancestor of *M. integrifolia* was dated to about 5.77 Ma (95% highest posterior density (HPD): 7.86--3.38 Ma) based on the cpDNA haplotypes ([Fig. 3a](#pone-0037196-g003){ref-type="fig"}) and around 5.33 Ma (95% HPD: 7.35--3.45 Ma) based on the ITS haplotypes ([Fig. 3b](#pone-0037196-g003){ref-type="fig"}). The crown age of the oldest clade II in the cpDNA phylogeny was dated to 3.92 Ma (95% HPD: 6.19--1.74 Ma), while the divergence times of the other main clades (I, III--VI) were estimated to range from 3.59 to 2.32 Ma. The divergence between clade I~nu~ and clade IV~nu~ was dated to 2.92 Ma (95% HPD: 4.45--1.73 Ma), and that between clade VII~nu~ and its sister clade was estimated to be 3.71 Ma (95% HPD: 5.37--2.26 Ma).

Discussion {#s4}
==========

Marked phylogeographic structure and genetic differentiation among highly fragmented populations of *M. integrifolia* and its implications for plant speciation in the QTP {#s4a}
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------

It is well known that complex topography could have great effects on the patterns of population genetic variation in space and time [@pone.0037196-Hewitt1]. To date, more than ten plant species in the QTP have been phylogeographically studied [@pone.0037196-Qiu1]. The common scenario is that the genetic variation was structured as a result of the interglacial re-colonization of the QTP platform from multiple glacial refugia, as retrieved by the mtDNA and/or cpDNA analyses of different kinds of plants. This conclusion was further supported by the perfect congruence of the mtDNA, cpDNA and nrDNA analyses on the conifer *Tsuga dumosa* [@pone.0037196-Cun1]. However, few studies revealed a pattern of population differentiation associated with the subtle landscape feature of the QTP [@pone.0037196-Qiu1], given that the majority of surveyed species are either conifers with strong dispersal ability or dominant species with a wide distribution on the plateau.

The distribution pattern of the 26 cpDNA haplotypes detected in *M. integrifolia* indicates a subtle geographic genetic structure and remarkable genetic differentiation among lineage-specific population groups of the species ([Fig. 2](#pone-0037196-g002){ref-type="fig"}), which is also supported by the results of analysis of molecular variance (AMOVA) ([Table 1](#pone-0037196-t001){ref-type="table"}), pairwise *F* ~ST~ ([Table 2](#pone-0037196-t002){ref-type="table"}), and the *U*-statistic of *G* ~ST~ and *N* ~ST~. As shown in [Fig. 2](#pone-0037196-g002){ref-type="fig"}, seven lineage-specific groups are well isolated in the eastern QTP. Three of them (ET, SLL & DX) are distributed exclusively in three mountain ranges with average altitudes above 4000 m above sea level, and the other four groups (WYN, LZ, EQa & EQb) inhabit in the edge of the QTP at relatively lower altitudes. That is, group ET includes five populations along the Nyainqentanglha Mountains, group SLL along the Shalulishan & Bayankalashan Mountains, and group DX in the Daxue Mountains & Qionglai Mountains. The exceptions, i. e., haplotypes U and Q in population 18 (group SLL) and haplotype T in population 19 (group SLL), might be resulted from gene flow from neighboring groups DX and EQa by seed dispersal in a low frequency. Analysis of the ITS data set also indicates a lineage-specific pattern of haplotype distribution, such as the three clades I~nu~, II~nu~, and VI~nu~ that correspond to groups WYN, LZ and ET, respectively ([Fig. 3b](#pone-0037196-g003){ref-type="fig"}). However, three lineage-specific groups (SLL, DX, and EQ) revealed in the cpDNA analysis are not supported by the ITS analysis. This could be attributed to the difference in gene flows through seed and pollen [@pone.0037196-Petit1], and may suggest that the mountain ranges could have acted as dispersal corridors for the native groups, although these corridors were isolated by unfavorable habitats, especially the heterogeneous environment along the deep-incised valleys ([Fig. 1](#pone-0037196-g001){ref-type="fig"}, [2](#pone-0037196-g002){ref-type="fig"}). On the other hand, these three mountains could also have blocked the spread of the four peripheral population groups (WYN, LZ, EQa & EQb) at relatively lower altitudes. The subtle phylogeographic structure among lineage-specific groups clearly indicates how the complex landscape, particularly the mountain ranges, has greatly affected the evolutionary history of the QTP plants ([Figs. 1](#pone-0037196-g001){ref-type="fig"}, [2](#pone-0037196-g002){ref-type="fig"}).

Moreover, individuals in relatively low altitudes (3000--4000 m) often inhabit in forest gaps or under shrub, a habitat quite different from the alpine scree at higher altitudes of the QTP. The great population genetic differentiation (*F* ~ST~) ([Table 2](#pone-0037196-t002){ref-type="table"}) and extensive morphological variation [@pone.0037196-Wu1] in *M. integrifolia* suggest that the population groups have been isolated for quite a long time, and the low genetic variation within population (3.00% from the AMOVA analysis) indicates that populations of the species have experienced genetic drift. An increasing number of studies indicated that habitat fragmentation could play important roles in the divergence and speciation for plants [@pone.0037196-Young1], [@pone.0037196-Fahrig1], [@pone.0037196-Templeton1]. Allopatric speciation is prone to occur in fragmented populations [@pone.0037196-Young1], whereas genetic drift resulting from the reduction in population size during habitat fragmentation could erase genetic diversity in isolated populations, leading to the decrease of population viability and extinction of species [@pone.0037196-Young1], [@pone.0037196-Templeton1]. Our molecular clock estimate based on cpDNA suggests that the divergence of main clades of *M. integrifolia* occurred during late Miocene (5.77 Ma) to early Quaternary (2.32 Ma) ([Fig. 3a](#pone-0037196-g003){ref-type="fig"}), a time span that is generally congruent with that estimated based on the ITS haplotypes (5.33--2.92 Ma) ([Fig. 3b](#pone-0037196-g003){ref-type="fig"}). The estimated divergence times of the main lineages are much earlier than the oldest Quaternary glaciation (about 0.71 Ma) [@pone.0037196-Zhou1] in the QTP. Therefore, the Quaternary climatic oscillations could not have contributed greatly to the divergence of the main clades of *M. integrifolia*, although it can not be completely ruled out. It is believed that the QTP reached an elevation like present at about 8.0 Ma [@pone.0037196-Harrison1], but was lowered by the following extensive faulting, and its recent rapid uplift occurred at around 3.6 Ma [@pone.0037196-Li1]. So, the divergence of the main clades of *M. integrifolia* during early to middle Quaternary is strongly correlated to the recent rapid uplift of the plateau [@pone.0037196-Li1]. That is, the rapid uplift of the QTP led to the strong fragmentation of the original continuous distribution of *M. integrifolia*, and the heterogeneous habitats along an elevation gradient hindered the dispersal of the species among the "sky islands", considering its limited dispersal ability (with insect-pollinated flower and gravity-dispersed seed). Given the uncertainty with the constraint age used in the present study, this scenario should be accepted with caution. However, our inference is consistent with previous reports that lineage divergence and speciation of plants in the QTP could be attributed to the topography changes resulting from the rapid uplift of the plateau [@pone.0037196-Wang3], [@pone.0037196-Zhang1], [@pone.0037196-Wang4]. As discussed above, the strong phylogeographic structure and great genetic differentiation among population groups of *M. integrifolia* may imply that habitat fragmentation arising from the uplift of the QTP has played a key role in speciation and extinction of the QTP plants.

Another interesting finding is the rough correlation between the haplogroups\' distribution and the regional climate in the QTP ([Fig. 2](#pone-0037196-g002){ref-type="fig"}). A typical pattern revealed from many previous studies [@pone.0037196-Qiu1] is that primitive lineages or ancestral haplotypes were distributed in the east edge of the QTP (refugia), from where northwestward re-colonization occurred. This pattern could be mostly attributed to the Quaternary climatic oscillations in the QTP [@pone.0037196-Qiu1]. However, evidence from the cpDNA and nuclear DNA (ITS) data set indicate that the older group (LZ) of *M. integrifolia* only occurs in southeast Tibet, with the distributions of younger groups (WYN, ET, SLL, DX, and EQ) extending northeastward in the eastern QTP ([Fig. 2](#pone-0037196-g002){ref-type="fig"}, [3](#pone-0037196-g003){ref-type="fig"}). There is palaeo-climatic evidence that the climatic shifts of the QTP and adjacent regions have been controlled by the Indian summer monsoon (southwest Asia monsoon) arising from the rapid uplift of the QTP [@pone.0037196-Wang5], [@pone.0037196-An1]. The Indian summer monsoon brings heavy summer rainfall for the south edge of the Himalayas, but leads to sharp decrease of precipitation in other area of the QTP due to the high topographic barrier ([Fig. 2](#pone-0037196-g002){ref-type="fig"}) [@pone.0037196-Wang5]. It is believed that the climatic divergence between the south edge and the other region of the QTP has developed since late Miocene when the Asia monsoon system initiated [@pone.0037196-An1]. The initial divergence of *M. integrifolia* dated to around 5.77 (95% HPD: 7.86--3.38) Ma based on cpDNA or 5.33 (95% HPD: 7.35--3.45) Ma based on ITS seems to have occurred synchronously with the climatic divergence in the QTP [@pone.0037196-An1]. The phylogeographic structure associated with the regional climate suggests that climatic divergence driven by the rapid uplift of the QTP likely led to the initial divergence of *M. integrifolia*.

Incongruence between cpDNA and nrDNA ITS haplotype phylogenies {#s4b}
--------------------------------------------------------------

Incongruence between phylogenies of nuclear and plastid DNA is generally caused by convergent evolution, lineage sorting, or hybridization/introgression [@pone.0037196-Rieseberg1]. Conflict between gene trees is a common phenomenon in plants, but few studies of the QTP plants revealed it and the underlying mechanism at the population level [@pone.0037196-Wang4]. In the present study, we found that the phylogenetic position of group WYN of *M. integrifolia* is very different between the cpDNA (clade I) and ITS (clade I~nu~) phylogenies constructed based on population sampling ([Fig. 3](#pone-0037196-g003){ref-type="fig"}), and the significance of this topological incongruence is also supported by the Shimodaira--Hasegawa test. Given the non-coding regions we used, convergent evolution should not be responsible for the reticulate evolution in *M. integrifolia*. Although incomplete lineage sorting could result in the stochastic fixation of ancestral polymorphisms in descendant populations, it is not supported by the strong geographic partition and early divergence of the cpDNA haplotypes, and the relatively antiquated ITS haplotypes of the population group WYN. However, hybridization frequently occurs in plants [@pone.0037196-Rieseberg1], and putative natural hybrids were reported in the genus *Meconopsis* [@pone.0037196-Taylor1]. Also, the alternation of topography resulting from the uplift of the QTP could have led to the contact of species originally distributed far from each other, providing opportunities for interspecific hybridization. So it is quite likely that the reticulate evolution was caused by hybridization between *M. integrifolia* and *M. betonicifolia*, and the chloroplast DNA of the latter was captured by the population group WYN of the former.

Although natural hybrids between *M. integrifolia* and its close relatives have never been reported so far, the cultivated hybrids ×*M. sarsonsii* (*M. integrifolia*×*M. betonicifolia*) and ×*M. harleyana* (*M. integrifolia*×*M. simplicifolia*) were reported about eighty years ago [@pone.0037196-Taylor1]. The hybrids could occur spontaneously or by deliberated crossing in the Royal Botanic Garden, Edinburgh. It is interesting that the flowers of the two hybrids were yellow as that of its parent *M. integrifolia* [@pone.0037196-Taylor1]. As a result, ×*M. harleyana* was identified as a yellow-flowered variety of *M. simplicifolia* rather than a hybrid, since it was often difficult to distinguish hybrids of *Meconopsis* only based on morphological characters [@pone.0037196-Taylor1]. In fact, the taxonomic status of the putative hybrid (group WYN) has confused botanists for about one hundred years owning to the intermediate morphological characters [@pone.0037196-Wu1], [@pone.0037196-Taylor1], [@pone.0037196-Prain1], [@pone.0037196-Zhang2], and it was even treated as a new species *M. pseudointegrifolia*, a variety of *M. integrifolia* or *M. pseudointegrifolia* [@pone.0037196-Taylor1], [@pone.0037196-Zhang2]. According to the present study, the most possible parents of the hybrid are *M. integrifolia* and *M. betonicifolia* ([Fig. 1](#pone-0037196-g001){ref-type="fig"}, [3](#pone-0037196-g003){ref-type="fig"}), but further investigation is necessary to corroborate this hypothesis considering the insufficient investigation of the genus in the geographically challenging region (QTP).

Field investigation found that the distributions of *M. integrifolia* and *M. betonicifolia* seldom overlap, and the former often inhabits at higher altitudes than the latter on the same mountain. Although the clade I includes four haplotypes widely distributed in northwest Yunnan ([Fig. 2](#pone-0037196-g002){ref-type="fig"}), the overlap distribution of the two species was only found in the Laojunshan Mountain. So it is unlikely that the hybridization has occurred recently, and the molecular dating also supports an ancient hybridization between the two species ([Fig. 3](#pone-0037196-g003){ref-type="fig"}).

Hybridization plays a great role in plant speciation, especially for groups with a high number of ployploids [@pone.0037196-Grant1]. The environmental and habitat heterogeneity are though to promote the rapid speciation by hybridization. Our phylogeographic analysis of *M. integrifolia* suggests that geographic isolation and hybridization could be two important mechanisms responsible for the divergence and speciation of *Meconopsis*, a species-rich genus with complex polyploids [@pone.0037196-Ying1], and even for speciation in other plant groups of the QTP.

Supporting Information {#s5}
======================

###### 

**Geographic origins, sample sizes, haplotypes and their frequencies of the 35** ***Meconopsis integrifolia*** **populations studied.**

(DOC)

###### 

Click here for additional data file.

We thank Dr. Dmitry A. Filatov, Dr. Maxim Kapralov and another anonymous reviewer for their insightful comments and suggestions on the manuscript. We also thank Prof. Jun Wen (Smithsonian Institution, National Museum of Natural History, USA) and her scientific expedition team (Drs. Ze-Long Nie, Yang Niu, Lei Xie, Bang Feng, Guo-Dong Li etc.) to Tibet for field work assistance; Dr. Ling-Feng Mao (Institute of Botany, CAS) for great help in the sample collection; and Ms. Wan-Qing Jin for assistance with lab work.

**Competing Interests:**The authors have declared that no competing interests exist.

**Funding:**This study was supported by the National Natural Science Foundation of China (Grant No. 30730010) and the National Basic Research Program of China (Grant No. 2007CB411602). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

[^1]: Conceived and designed the experiments: FSY XQW. Performed the experiments: FSY ALQ. Analyzed the data: FSY XQW. Contributed reagents/materials/analysis tools: FSY YFL. Wrote the paper: FSY XQW.
