Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Chloroplast DNA Phylogeography Reveals Repeated Range Expansion in a Widespread Aquatic Herb Hippuris vulgaris in the Qinghai-Tibetan Plateau and Adjacent Areas

  • Jin-Ming Chen ,

    Contributed equally to this work with: Jin-Ming Chen, Zhi-Yuan Du

    Affiliation Key Laboratory of Aquatic Botany and Watershed Ecology, Wuhan Botanical Garden, The Chinese Academy of Sciences, Wuhan, Hubei, China

  • Zhi-Yuan Du ,

    Contributed equally to this work with: Jin-Ming Chen, Zhi-Yuan Du

    Affiliation Key Laboratory of Aquatic Botany and Watershed Ecology, Wuhan Botanical Garden, The Chinese Academy of Sciences, Wuhan, Hubei, China

  • Shan-Shan Sun,

    Affiliations Key Laboratory of Aquatic Botany and Watershed Ecology, Wuhan Botanical Garden, The Chinese Academy of Sciences, Wuhan, Hubei, China, Graduate University of the Chinese Academy of Sciences, Beijing, China

  • Robert Wahiti Gituru,

    Affiliation Botany Department, Jomo Kenyatta University of Agriculture and Technology, Nairobi, Kenya

  • Qing-Feng Wang

    qfwang@wbgcas.cn

    Affiliation Key Laboratory of Aquatic Botany and Watershed Ecology, Wuhan Botanical Garden, The Chinese Academy of Sciences, Wuhan, Hubei, China

Abstract

Background

The Qinghai-Tibetan Plateau (QTP) is one of the most extensive habitats for alpine plants in the world. Climatic oscillations during the Quaternary ice age had a dramatic effect on species ranges on the QTP and the adjacent areas. However, how the distribution ranges of aquatic plant species shifted on the QTP in response to Quaternary climatic changes remains almost unknown.

Methodology and Principal Findings

We studied the phylogeography and demographic history of the widespread aquatic herb Hippuris vulgaris from the QTP and adjacent areas. Our sampling included 385 individuals from 47 natural populations of H. vulgaris. Using sequences from four chloroplast DNA (cpDNA) non-coding regions, we distinguished eight different cpDNA haplotypes. From the cpDNA variation in H. vulgaris, we found a very high level of population differentiation (GST = 0.819) but the phylogeographical structure remained obscure (NST = 0.853>GST = 0.819, P>0.05). Phylogenetic analyses revealed two main cpDNA haplotype lineages. The split between these two haplotype groups can be dated back to the mid-to-late Pleistocene (ca. 0.480 Myr). Mismatch distribution analyses showed that each of these had experienced a recent range expansion. These two expansions (ca. 0.12 and 0.17 Myr) might have begun from the different refugees before the Last Glacial Maximum (LGM).

Conclusions/Significance

This study initiates a research on the phylogeography of aquatic herbs in the QTP and for the first time sheds light on the response of an alpine aquatic seed plant species in the QTP to Quaternary climate oscillations.

Introduction

Profound climatic oscillations during the Quaternary ice age of the past ∼2 million years have left imprints on modern-day biota by effecting species' distributions, evolution and extinction [1][5]. Such range and demographic changes undoubtedly affect the geographic genetic structure of species [4]. Molecular-based phylogeographical methods have been widely used to unravel the evolutionary history of species by tracing the genetic footprints of Quaternary climatic shifts in different regions, particularly in Europe, North America, and Japan [5][7]. Numerous phylogeographical surveys in Europe, North America, and Japan have identified the locations of glacial refugia and routes of colonization/range expansion of organisms after glacial periods in these regions [5]. In contrast, knowledge of phylogeographical histories of organisms occurring in mainland China as well as the association of this history with climatic fluctuations, is relatively limited, particularly for plants [5], [8]. Such information would clarify why China is such an important center of global biodiversity [9], [10].

The Qinghai-Tibetan Plateau (QTP) located in western China is the highest and largest plateau on earth, with an average altitude of approximately 4,500 m and an area of 2.5×106 km2 [11]. The phased uplift of the QTP, which probably started with the collision of Indian and Eurasian plates about 50 million years ago [12][15], has greatly changed the geology and topography of East Asia, resulting in the significant climatic and environmental changes witnessed since the early Cenozoic Era [16][18]. The QTP has experienced several glacial and interglacial cycles since the Pliocene [19], [20]. The recent uplifts of much of the QTP during the mid-to-late Quaternary may have helped to drive the onset of mountain glaciers, ice caps, and valley glaciers [20][23]. However, it is unlikely that a unified ice sheet developed on the QTP in the late Pleistocene and Quaternary periods [23], [24]. Although no unified massive ice sheet could have developed in the QTP during glacial periods [23], [24], the alternating glacial and interglacial periods have strongly and periodically influenced the environments and vegetation of the QTP [24], [25].

It is estimated that the present-day alpine flora of the QTP (above 4,000 m) numbers ca.1,816 seed plant species in more than 339 genera, with ca.33 genera being endemic [26]. The Hengduan Mountains (HDM, reaching 7,556 m), located at the southeastern edge of the QTP [27], represent one of the world's biodiversity hotspots [28]. The vegetation in the QTP and adjacent highlands is considered to be highly sensitive and vulnerable to global climate change because plant growth and distribution in the region depend greatly on survivable temperatures [11]. The intensity of the uplift of QTP in the late Pliocene combined with repeated climate fluctuations during the Quaternary have dramatically influenced the distribution and evolution of many plant species in these regions [22], [26]. However, how the distribution ranges of plant species shifted on the QTP in response to Quaternary climatic changes is still in question. Recent phylogeographical studies have indicated remarkable differences in the demographic history of the plant species from the QTP. Some studies provide evidence favoring the hypothesis that the extant plant populations on the plateau are colonists from the low altitude refugia at the plateau edge (e.g., the HDM) during the Last Glacial Maximum (LGM) and /or previous glacial periods [29][32]. In contrast, several other studies support the hypothesis that species probably survived on high-altitude parts of the central Plateau throughout the Quaternary [33][39]. Alternatively, additional studies have indicated that the tectonic events rather than the Quaternary climatic oscillation may have played the most important role in plant divergence and speciation in the QTP [10], [40][42]. Species with different ecological requirements and distributions may have led to dissimilarities in refugia and migration routes [43]. The current plant phylogeographical studies on the QTP have left large gaps in taxonomic and habitat sampling; for example, as yet there are only few published data examining phylogeographical patterns on herbaceous plants [31][33], [38], [39], [42], [44], with no research on aquatic or wetland plant phylogeography.

In the QTP and adjacent areas, the large number of high altitude lakes, ponds, streams, rivers and wetlands [45] support ca.133 aquatic herbaceous seed plants species in 63 genera and 29 families [46]. Eleven of these species are endemic to China but a few (four species, e.g., Puccinellia degeensis, P. multiflora, Carex sagaensis, C. muliensis) are endemic to the QTP and adjacent areas [46]. Species richness of the aquatic plants exhibits a tendency whereby there is gradual decrease in total species number from the southeast edge of the QTP (e.g., the HDM) to the northwest region [46]. In the northern region, the flora of aquatic plants comprises only a few widespread and tolerant species [46]. Aquatic plants typically have much broader geographical distributions and ranges compared to their terrestrial counterparts [47], [48]. As opposed to trees, the range shifts of herbaceous aquatic plants could have responded more quickly to climate or vicariance [48], [49]. However, for aquatic plants in the QTP and adjacent regions, patterns of range shifts following climate change or vicariance are poorly understood.

In this study, we assessed the phylogeographical pattern of an aquatic perennial herb, Hippuris vulgaris (Hippuridaceae), a species with a worldwide distribution, in the QTP and adjacent areas. In China, the species is mainly distributed in the Southwest, Northeast and Northwestern regions [50]. The plant grows in ponds, ditches, shallow riverbeds, and other freshwater wetlands. Hippuris vulgaris is a major and widespread component of aquatic vegetation in the QTP and adjacent areas. It can be found at altitudes of up to 5,000 m and is a highly cold tolerant species [46]. We investigated the variations in chloroplast DNA (cpDNA) non-coding regions (ycf6-psbM, trnT-trnL, the rps16 intron and atpI-atpH) between 385 individuals from 47 populations from the QTP and adjacent areas. Our objectives were: (1) to elucidate the intraspecific divergence of H. vulgaris in relation to the Quaternary climatic fluctuations, and, (2) to infer historical population range expansion of this aquatic plant on the QTP.

Materials and Methods

Ethics Statement

This study was conducted in accordance with all People's Republic of China laws. No specific permits were required for the described field studies. No specific permissions were required for access to the locations described in this study. The location is not privately owned and neither is it protected in any way. This study species is also not protected by any law.

Plant sampling

Leaf material from each individual of H. vulgaris was collected from populations across the QTP and neighboring areas, from the southeast (the HDM) to the northeast (Qinghai province) to the north (Xinjiang province) (Fig. 1). The sampled individuals within populations were at least 3 m apart by distance. Collected fresh leaves were dried immediately with silica gel and stored at room temperature. In total, 385 individuals from 47 populations were included in the study. Details on the collection sites including the latitude, longitude, altitude and the sample size are provided in Table 1. Voucher specimens from all populations were deposited in the Wuhan Botanical Garden (WBG) herbarium.

thumbnail
Figure 1. Collection localities (population codes as in Table 1) and the geographical distributions of eight chloroplast haplotypes (A–H) found in 47 populations of H. vulgaris in the Qinghai-Tibetan Plateau and neighboring areas.

https://doi.org/10.1371/journal.pone.0060948.g001

thumbnail
Table 1. Details of sample locations, samples size (N), haplotypes and haplotype diversity (h) of 47 populations of Hippuris vulgaris surveyed for DNA sequence variation at four combined chloroplast regions.

https://doi.org/10.1371/journal.pone.0060948.t001

DNA extraction, PCR amplification and sequencing

Total genomic DNA was isolated from 0.5 g of silica-dried leaf tissue following the CTAB method [51]. After preliminary screening of 20 non-coding regions of cpDNA among six individuals from six geographically distant populations of H. vulgaris, polymorphism was observed in four regions, ycf6-psbM [52], trnT-trnL [53], rps16 intron [52] and atpI-atpH [54]. These four cpDNA fragments were amplified and sequenced in all individuals sampled. PCR reactions were carried out in a volume of 50 µl containing 0.25 mM of each dNTP, 5 µl of 10×Taq buffer [10 mM Tris-HCl (pH 8.3), 1.5 mM MgCl2, and 50 mM KCl], 1 mM of each primer, 2 U Taq Polymerase (TaKaRa, Dalian, China), and 60 ng of DNA template. Amplification of genomic DNA was carried out on a PTC-100™ thermocycler (Bio-Rad, Hercules, CA, USA). It was commenced with 5 min at 94°C and annealing at 60°C for 5 min, followed by 30 cycles of 1 min at 94°C, 2 min annealing at 50–60°C (depending on the kinds of primers), and 2 min extension at 72 °C, and a final extension cycle of 10 min at 72°C.

The size of PCR products was determined electrophoretically on 1.5% (W/V) agarose gels run at 100 V in 0.5×TBE (Tris-boric acid-EDTA) and visualized by staining with ethidium bromide. All PCR products were purified from agarose gel using the TIANquick Midi Purification Kit following the protocols provided by the manufacturer (Tiangen, Beijing, China). The purified PCR products were sequenced using the ABI Prism BigDye terminator cycle sequencing ready reaction kit (Applied Biosystems, Foster City, CA, USA) and performed on an ABI 3730 automated sequencer (Applied Biosystems, Foster City, CA, USA).

Phylogenetic analysis

CpDNA sequences were aligned by the program Clustal W [55] and then adjusted manually. Insertions/deletions (indels) were firstly treated as a single character resulting from one mutation. CpDNA haplotypes were determined from nucleotide substitutions and indels of the aligned sequences of ycf6-psbM, trnT-trnL, the rps16 intron, and atpI-atpH. For a comparison, the cpDNA haplotypes were also determined from the substitutions only. Phylogenetic analyses of cpDNA haplotype sequences were performed by the neighbor-joining (NJ) method. The NJ tree was constructed according to Kimura's two-parameter model [56] using MEGA 3 [57]. The mid-point rooting method which assuming approximately equal evolutionary rates over all the branches was used to root the NJ tree. Confidence values at the nodes were tested by performing 1,000 bootstrap replicates [58]. Phylogenetic analyses of cpDNA haplotype sequences were also reconstructed using a Bayesian approach. The best-fit model of DNA substitutions (F81) was selected with the Akaike information criterion (AIC) in MrModeltest v2.3 [59]. The Bayesian Markov Chain Monte Carlo (BMCMC) estimate of phylogeny was inferred using MrBayes v. 3.1.2 [60] with default priors. Bayesian analysis was performed for 10,000,000 generations. Two simultaneous chains were run with sampling every 1,000 generations. Tracer v1.5 [61] was used to check the likelihood scores of sampled trees. The first 25% trees were discarded and > 50% posterior probability consensus trees from the remaining trees were calculated. Pairwise differences between DNA haplotypes were calculated using ARLEQUIN ver. 3.1 [62]. A statistical parsimony haplotype network based on the matrix of pairwise differences between DNA haplotypes was obtained with the aid of the TCS 1.06 [63] using the 95% connection probability limit and treating gaps as single evolutionary events.

Population and phylogeographic analyses

To measure levels of genetic variation, haplotype diversity (h) and nucleotide diversity (pi) for each population and for the entire species were calculated using ARLEQUIN ver. 3.1 [62]. The average gene diversity within populations (HS), total gene diversity (HT), proportion of total diversity caused by differences between populations (GST), and the number of substitution types (NST) were estimated from cpDNA haplotypes following the methods described by Pons and Petit [64] using the program PERMUT (available at http://www.pierroton.inra.fr/genetics/labo/Software/Permut/; 1000 permutations test). A comparison was made between GST and NST using the U-statistic, which is approximated by a Gaussian variable taking into account the covariance between GST and NST. GST depends only on the haplotype frequency, whereas NST takes into account differences between haplotypes. If NST is significantly higher than GST, closely related haplotypes occur more often in the same populations than less closely related haplotypes, indicating the presence of phylogeographical structure [64].

The spatial genetic structure of haplotypes was analyzed by spatial analysis of molecular variance (SAMOVA) using SAMOVA v. 1.0 [65]. SAMOVA iteratively seeks the composition of a user-defined number (K) of groups of geographically adjacent populations that maximizes FCT, i.e., differences among groups of populations. The most likely number of groups (K) was determined by repeatedly running SAMOVA with 2–10 groups and choosing those partitions with a maximum FCT value, as suggested by Dupanloup et al. [65]. For each K, the configuration with the largest FCT value after the 100 repetitions was retained as the best grouping of populations. Hierarchical analysis of molecular variance (AMOVA; [66]) was also conducted to quantify the proportion of total genetic variance explained by differences between regional population groups (as identified by SAMOVA) and between populations within groups. The AMOVA was performed with the program ARLEQUIN ver. 3.1 [62] and significance of variance components were tested with 10, 000 permutations.

The F statistic (FST) between populations was calculated using the program ARLEQUIN ver. 3.1 [62], and 1,000 permutations were used for significance testing in all cases. To test the significance of isolation by distance between populations, the Mantel test with 1,000 random permutations on matrices of pairwise population FST values and the geographical distances was performed using ARLEQUIN ver. 3.1 [62].

To infer the historic demographic expansion events within each inferred haplotype clade, we carried out mismatch distribution analyses in ARLEQUIN ver. 3.1 [62]. The goodness-of-fit under a sudden-expansion model was tested with the sum of squared deviations (SSD) and Harpending's raggedness index (Rag) [67]. A parametric bootstrap approach [68] with 1,000 replicates was used to test the observed mismatch distribution's fit to the sudden expansion model. If the sudden-expansion model was not rejected, the expansion parameter (τ) was converted to an estimate of time (T, in number of generations) since expansion began using T  = τ/2u [69], [70], where τ is the mode of the mismatch distribution and u is the neutral mutation rate for the entire sequence per generation. The value for u was calculated as u = µkg, where µ is the substitution rate in substitutions per site per year (s/s/y), k is the average sequence length of the DNA region under study (each indel was treated as one mutation) and g is the generation time in years (i.e. the average age at which reproduction first occurs). Rate-constancy of cpDNA haplotype evolution in H. vulgaris was evaluated by relative rate tests in MEGA 3 [57]. After discovery of rate-constancy, the average cpDNA mutation rate of 1.52×10−9 substitutions per site per year (s/s/y) [71] was used to estimate the expansion time. This substitution rate value used here is approximate estimate but fall within the range commonly reported for non-coding chloroplast regions of seed plants [72], [73]. The generation time of H. vulgaris was assumed to be one year (J.M. Chen, pers. obs.).

Divergence time between lineages

A Bayesian approach was used to estimate the divergent times for different H. vulgaris lineages using BEAST v 1.5.4 [74]. CpDNA haplotype sequences were found to follow the GTR+Γ+G model as inferred from MrModeltest v. 2.3 [59] with six Gamma categories; a starting tree was randomly generated and a Yule process was performed prior. We assumed an evolutionary rate of 1.52×10−9 s/s/y [71] for the four combined chloroplast non-coding regions. This rate was then used for estimating lineage divergence time in BEAST under an uncorrelated lognormal relaxed molecular clock assumption. Two separate MCMC analyses were run for 10,000,000 generations with sampling at every 1,000 generations to ensure all the effective sample size (ESS) values greater than 200. Tracer v1.5 [61] was used to check the parameters and the first 25% generations were discarded as burn-in.

Results

Chloroplast DNA diversity

The aligned sequences of ycf6-psbM, trnT-trnL, the rps16 intron and atpI-atpH were 515, 427, 513, and 651 bp in length, respectively. The total length of combined alignments was 2,106 bp. A total of eight polymorphic sites (four substitutions and four indels) of the combined four cpDNA non-coding regions resulted in the resolution of eight haplotypes (A–H) across the 385 individuals (47 populations) of H. vulgaris sampled (Table 1 and Table 2). GenBank accession number (Nos. JX415319 to JX425330) for each of the cpDNA non-coding region of each chloroplast haplotype is listed in Table S1. At the species level, the estimated haplotype diversity was h = 0.604. Haplotype diversity varied across populations, ranging from 0.000 to 0.733, the QHH and KS populations have the highest h values (Table 1). Nucleotide diversity was estimated within the species as a whole (pi = 0.0020±0.0010) and within populations, ranging from 0.0000 to 0.2167, the QHH population has the highest pi value (data not shown). Within-population gene diversity (HS) was 0.111 and total gene diversity (HT) was 0.615. Parameters of VS and VT were 0.071 and 0.484, separately. From the 47 populations sampled, 35 were fixed for a single chlorotype and only 12 populations were polymorphic with two or three haplotypes (Fig. 1, Table 1). Among the eight cpDNA haplotypes revealed in H. vulgaris in this study, haplotypes A and D were found to be the most widespread. Thirty of the 47 populations investigated harbored haplotype A and 18 populations harbored haplotype D. The haplotype E was found in five populations (Fig. 1, Table 1). When the indels were excluded, there were only five haplotypes remained (A, C, D, F and H) (Table 2). The distribution pattern of these haplotypes was similar to that of haplotypes based on substitutes and indels, e.g., haplotypes A and D were found to be widely distributed in the studied populations; 10 populations harbored two or three haplotypes, and haplotypes F and H were unique to a particular population, separately (Fig. 2). Thus, the following analyses of population structure and phylogeny were based only on the haplotypes realized from substitutes and indels.

thumbnail
Figure 2. The geographical distributions of the five chloroplast haplotypes (A, C, D, F and H) based on substitutions in only 47 populations of H. vulgaris in the Qinghai-Tibetan Plateau and neighboring areas.

https://doi.org/10.1371/journal.pone.0060948.g002

thumbnail
Table 2. Chloroplast DNA sequence polymorphism detected in four non-coding regions of Hippuris vulgaris identifying eight haplotypes (A–H).

https://doi.org/10.1371/journal.pone.0060948.t002

Population structure

The test for phylogeographical structure in the variation of cpDNA haplotype showed that

NST was not significantly higher than GST (0.853 and 0.819, respectively) indicating no phylogeographic structure. In the SAMOVA analyses, FCT values increased progressively as the value of K was increased from 2 to 10, which did not allow us to unambiguously identify the number K of groups of populations. The Mantel test revealed no significant correlation between genetic and geographical distances of cpDNA haplotypes (r = −0.01582, P = 0.575) over all populations.

Phylogenetic relationships and geographical distribution of haplotypes

The neighbor-joining tree of the eight cpDNA haplotypes detected from the 385 individuals of H. vulgari (Fig. 3A) was identical to the tree reconstructed under a Bayesian approach. The unrooted TCS network of all H. vulgaris haplotypes (Fig. 3B) is consistent with the trees created from the NJ and Bayesian methods. Among these haplotypes, haplotype A was widely distributed in all the sampled regions (30 populations). Haplotype F was unique to population QHH from the east edge of the Plateau and G was confined to population KS from the interior of the Plateau. Haplotype B was found in one population (HBH) from the Xinjiang province and one population (RT) from the interior of the Plateau. Haplotype C was also found in three populations from the interior region (Fig. 1); HPG II contained three haplotypes (D, E and H) (Fig. 3). Among these haplotypes, haplotype D was widely distributed on the east edge of the QTP, e.g., the HDM (10 populations) and all the sampled populations from Qinghai province (3 populations). In this region, haplotype H was unique to population GZ5 and haplotype E was found in two populations (WXH and GZ4). In addition, haplotype D was found in five populations from the interior of the QTP, e.g., Xizang province. Haplotype E was also found in four populations from the interior region (Fig. 1).

thumbnail
Figure 3. Phylogeny of the eight chloroplast haplotypes (A–H) detected in H. vulgaris.

(A) Neighbor-joining clustering of the eight chloroplast haplotypes. Numbers above branches indicate the bootstrap values (>50% are shown) for NJ (Left, 1,000 replicates) and Bayesian analyses (Right). HPG: Haplogroup; (B) 95% plausible network of the eight chloroplast haplotypes. Each solid line represents one mutational step that interconnects two haplotypes for which parsimony is supported at the 95% level. The distribution of a certain haplotype is marked in the circles. The small open circle indicates an inferred intermediate haplotype not deteched in this investigation. The size of each circle is proportional to the haplotype frequency.

https://doi.org/10.1371/journal.pone.0060948.g003

Demographic analyses

Mismatch analyses indicated that the distribution of pairwise differences for HPG I and HPG II were unimodal. Both the observed variance (SSD) and the HRag for each of the two groups was not significantly different from that expected under the population expansion model (P>0.05), which further supported that HPG I and HPG II had experienced rapid expansions. The expansions for HPG I and HPG II were estimated to have been at 0.12 and 0.17 Myr (Table 3).

thumbnail
Table 3. Mismatch distribution analyses and estimation of expansion time.

https://doi.org/10.1371/journal.pone.0060948.t003

Lineages divergence time

The divergence times between the haplotypes estimated by Bayesian method ranged from 0.480 Myr to 0.026 Myr (Fig. 4). The divergence time of the two main groups of haplotypes (HPG I and II) were estimated at 0.480 Myr, which suggests that the divergence among haplotypes of H. vulgaris falls into mid-to-late Pleistocene.

thumbnail
Figure 4. Bayesian analyses of the divergent time (Myr) between chloroplast DNA haplotypes.

https://doi.org/10.1371/journal.pone.0060948.g004

Discussion

For the aquatic herb H. vulgaris from the QTP and adjacent areas, inter-population differentiation was high, with GST = 0.819 occurred among populations, indicating high genetic differentiation across the sampled region. Unlike most previous studies, which show high genetic differentiation coupled with distinct phylogeographical structure [29], [75], the phylogeographical structure among the sampled populations of H. vulgaris was not obvious. Although two haplotype groups were recognized from the NJ tree, Bayesian methods and the haplotype network, the distributions of haplotype groups did not correspond to distinct geographical areas. SAMOVA analyses also failed to detect any meaningful geographical groups of the sampled populations based on the cpDNA haplotypes. A lack of significant geographical structure of populations was also confirmed by the results of haplotype-identification permutation tests, which indicated that NST was not significantly larger than GST (NST = 0.853>GST = 0.819, P>0.05). Lack of significant phylogeographical structure in H. vulgaris from the QTP and adjacent areas may result mainly from the occurrence of widespread haplotypes, i.e., haplotypes A and D, across the sampled regions. As a widespread, clonal and wind-pollinated herbaceous plant species, H. vulgaris populations were expected to be less differentiated [76], [77]. However, contrary to the expectation, populations were greatly differentiated although they shared several common haplotypes. The high population differentiation was probably because of the different frequencies of these common haplotypes occurred in each region.

In this study, the divergence times for all the haplotypes were estimated to be less than 0.5 Myr. The genetic divergence between the two intraspecific lineages of the species occurred at 0.48 Myr, corresponding to the mid-to-late Pleistocene [20]. This divergence within the species falls within the most extensive glaciations during the Naynayxungla Glaciation (0.8–0.5 Myr), when there were many large ice caps, glacier complexes and great valley glaciers, covering a total area≥500,000 km2 [20]. This extensive glacial advance likely continued until 0.17 Myr after the penultimate glaciations (0.3–0.13 Myr) [20], [21][23]. During these glaciations, the ice sheet is thought to have covered an area five to seven times larger than it does today, and even during the interglacial warming periods, the central QTP region retained glaciers [20], [78][80]. The glaciers and/or extremely low temperature in the high mountains during these stages might have created barriers to gene flow between geographically isolated populations of H. vulgaris, which therefore promoted its intraspecific divergence.

For each of the intraspecific lineages, our mismatch distribution analyses showed a unimodal distribution, suggesting a recent range expansion had occurred. Both the observed variance (SSD) and the HRag for each of the two groups was also not significantly different from that expected under the population expansion model. In addition, the haplotype network analysis revealed two “star-like” clusters originating from the ancestral haplotypes A and D, respectively, which also supports range expansion. These two expansions were dated to be 0.12 and 0.17 Myr, respectively, which fell within the interglacial stage between the Penultimate Glaciation and the Last Glaciation [20]. Previous studies showed that glacial retreat has occurred since 0.17 Ma [20], [22]. Our results suggest that H. vulgaris populations expanded as a consequence of this retreat. The widespread haplotypes A and D in each lineage probably originated from the recent expansions during the last interglacial period. Hippuris vulgaris is a clonal aquatic herb; vegetative reproduction of this species is by rhizomes, which is potentially a mode of long distance colonization. The simple haplotype composition (many populations were composed of a single haplotype) might be attributed to founder effects. Most of the haplotypes at tip positions were unique to their particular population. Haplotypes A and D may represent the ancestral genotypes and probably produced other mutational derivatives following recent expansions [81], [82].

The two independent range expansions could have begun from different refugee populations. The phylogeographical pattern of cpDNA haplotypes demonstrates two “star-like” clusters that were the result of haplotypes being linked to a central haplotype, but this relatively simple pattern can be explained by populations preserved in refugia that have experienced population expansion after glaciations [83]. However, no specific refugia areas could be deduced because there was no significant geographical structure of the genetic variation (although several haplotypes were unique to a single population). The HDM range of the southeastern QTP is known as an important glacial refugium of many plants [27]. Based on the aquatic plant species richness of different regions of the QTP and adjacent areas, e.g., the total species number decreasing gradually from the HDM region to the interior and to the northwest regions of the plateau [46], we also supported the hypothesis that the HDM region might be a refugium for aquatic plants and that aquatic plant species on the QTP might have colonized their present-day distribution ranges from this refugium. However, we could not exclude the possibility of H. vulgaris survival in refugia in the interior QTP, because this region possesses a similar cpDNA haplotype diversity and haplotype uniqueness as the HDM region. Plant species survival throughout the Quaternary on the QTP has also been reported in several previous studies [33][39]. Six populations distributed in the central region (e.g., LZ population), northeast region (e.g., QHH population) and also the southeast edge (e.g., MGC, LZG, GZ1 and HY2 populations) of the QTP possessed the two widespread haplotypes A and D indicating the populations' admixture might be a result of two independent range expansions.

In three adjacent populations (HBH, QH1 and QH2) of the QTP, only two haplotypes were found. One is the widespread haplotype A and another is the narrow distributed and newly generated haplotype B, which is shared only by HBH and RT populations. Long distance colonization throughout the QTP and northward to the Northwest region of China might be a possible reason for the occurrence of these two haplotypes in northern region of the QTP, which is consistent with previous findings in other plants, e.g., Pedicularis longiflora [32] and Lepisorus clathratus [8]. However, these two populations (HBH and RT) currently are separated by the Tarim basin, which contains the vast sandy Taklimakan desert, the driest desert in Asia. This region had already begun aridification by the early Pleistocene; the desert and gobi terrain further expanded and the dry climate continued to develop during the middle Pleistocene [84]. Thus, it seems unlikely that H. vulgaris had colonized through this wide arid region in Northern China, which would have acted as a barrier to dispersal. It was possible that H. vulgaris had colonized the northwest region of China through other routes, e.g., from the western or eastern edge of the Tarim basin, but this hypothesis still needs further study, with more samples around the Tarim basin included.

In conclusion, from the phylogeographic analyses of H. vulgaris, we revealed two independent range expansions in response to the Quaternary climatic oscillations on the QTP and adjacent areas. These two expansions might have begun from different refuge populations from before the LGM, which suggested that the H. vulgaris had survived the LGM and/or previous glaciers since its origins. These findings also partly support the assumption that there was “no massive ice-sheet on the QTP” [24] at least during LGM of the Quaternary. Although our findings on the aquatic herb H. vulgaris were basically similar to the results obtained for several terrestrial herbs [31][33], [38], [39], our study reveals the range shifts during glacial and interglacial periods of aquatic plants occurring on the QTP for the first time, which improves our understanding of phylogeographical patterns of aquatic plants in the QTP and adjacent areas. Although our methods are able to show only the latest range expansions, it seens probable that expansions also occurred during earlier interglacials. Admixture of populations due to repeated range expansions, therefore, might be an explanation of why no significant phylogeographical structure was found and no specific refugia areas could be deduced. Whether other aquatic plant species also display similar evolutionary histories on the QTP as revealed by H. vulgaris remains a future research question.

Supporting Information

Table S1.

GenBank accession numbers for each of the four non-coding regions of each chloroplast haplotype (A-H) identified in Hippuris vulgaris.

https://doi.org/10.1371/journal.pone.0060948.s001

(DOC)

Acknowledgments

The authors thank Dun Wang, Kuo Liao, Ling-Yun Chen, Shu-Ping Gu, Zhi-Cheng Long, Ming Chen for their help in fieldwork, Yang-Yang Yuan, Hua Liu and Qi Mi for their assistance in the laboratory, Victoria L. Sork (University of California, Los Angeles, USA) for correcting the English language, and the two anonymous reviewers for their constructive comments.

Author Contributions

Conceived and designed the experiments: JMC QFW. Performed the experiments: JMC ZYD SSS. Analyzed the data: JMC ZYD. Contributed reagents/materials/analysis tools: JMC ZYD. Wrote the paper: JMC RWG.

References

  1. 1. Hewitt GM (1996) Some genetic consequences of ice ages and their role in divergence and speciation. Biol J Linn Soc 58: 247–276.
  2. 2. Hewitt GM (2000) The genetic legacy of the Quaternary ice ages. Nature 405: 907–913.
  3. 3. Hewitt GM (2004) Genetic consequences of climatic oscillations in the Quaternary. Phil Trans Roy Soc London Ser B Biol Sci 359: 183–195.
  4. 4. Comes HP, Kadereit JW (1998) The effects of Quaternary climatic changes on plant distribution and evolution. Trends Plant Sci 3: 432–438.
  5. 5. Qiu YX, Fu CX, Comes HP (2011) Plant molecular phylogeography in China and adjacent regions: Tracing the genetic imprints of Quaternary climate and environmental change in the world's most diverse temperate flora. Mol Phylogenet Evol 59: 225–244.
  6. 6. Soltis DE, Morris A, McLachlan JS, Manos PS, Soltis PS (2006) Comparative phylogeography of unglaciated eastern North America. Mol Ecol 15: 4261–4293.
  7. 7. Hickerson MJ, Carstens BC, Cavender-Bares J, Crandall KA, Graham CH, et al. (2010) Phylogeography's past, present and future: 10 years after Avise, 2000. Mol Phylogenet Evol 54: 291–301.
  8. 8. Wang L, Wu ZQ, Bystriakova N, Ansell SW, Xiang QP, et al. (2011) Phylogeography of the Sino-Himalayan fern Lepisorus clathratus on “The roof of the world”. PloS One 6: e25896.
  9. 9. Ying TS (2001) Species diversity and distribution pattern of seed plants in China. Biodiver Sci 9: 393–398.
  10. 10. Lopez-Pujol J, Zhang FM, Ge S (2006) Plant diversity in China: richly varied, endangered, and in need of conservation. Biodivers Conser 15: 3983–4026.
  11. 11. Zheng D (1996) The system of physico-geographical regions of the Qinghai–Tibet (Xizang) Plateau. Sci China (Ser D) 39: 410–417.
  12. 12. Molnar P, Tapponnier P (1975) Cenozoic tectonics of Asia: Effects of a continental collision. Nature 189: 419–426.
  13. 13. Harrison TM, Copeland P, Kidd WSF, Yin AN (1992) Raising Tibet. Science 255: 1663–1670.
  14. 14. Tapponnier P, Xu Z, Roger F, Meyer B, Arnaud N, et al. (2001) Oblique stepwise rise and growth of the Tibet Plateau. Science 294: 1671–1677.
  15. 15. Royden LH, Burchfiel BC, Van Der Hilst RD (2008) The geological evolution of the Tibetan Plateau. Science 321: 1054–1058.
  16. 16. Ruddiman WF, Kutzbach JE (1991) Plateau uplift and climate change. Sci Am 264: 66–74.
  17. 17. Raymo ME, Ruddiman WF (1992) Tectonic forcing of Late Cenozoic climate. Nature 359: 117–122.
  18. 18. An ZS, Kutzbach JE, Prell WL, Porter SC (2001) Evolution of Asian monsoons and phased uplift of the Himalaya-Tibetan Plateau since Late Miocene times. Nature 411: 62–66.
  19. 19. Li JJ, Shi YF, Li BY (1995) Uplift of the Qinghai-Xizang (Tibet) Plateau and Global change. Lanzhou: Lanzhou University Press.
  20. 20. Zheng BX, Xu QQ, Shen YP (2002) The relationship between climate change and Quaternary glacial cycles on the Qinghai–Tibetan Plateau: review and speculation. Quat Internat 97: 93–101.
  21. 21. Zheng Z, Yuan BY, Petit-Maire N (1998) Paleoenvironments in China during the last glacial maximum and the Holocene optimum. Episodes 21: 152–158.
  22. 22. Zhang DF, Li FQ, Ben JM (2000) Eco-environmental effects of the Qinghai–Tibet plateau uplift during the Quaternary in China. Environ Geol 39: 1352–1358.
  23. 23. Shi YF (2002) Characteristics of late Quaternary monsoonal glaciation on the Tibetan Plateau and in East Asia. Quatern Internat 97/98: 79–91.
  24. 24. Shi YF, Li JJ, Li BY (1998) Uplift and Environmental Changes of Qinghai-Tibetan Plateau in the Late Cenozoic. Guangzhou, China: Guangdong Science and Technology Press.
  25. 25. Tang LY, Shen CM (1996) Late Cenozoic vegetational history and climatic characteristics of Qinghai–Xizang Plateau. Acta Micropalaeon Sin 13: 321–337.
  26. 26. Wu SG, Yang YP, Fei Y (1995) On the flora of the alpine region in the Qinghai-Xizang (Tibet) plateau. Acta Bot Yunnan 17: 233–25.
  27. 27. Wu ZY (1988) Hengduan Mountain flora and her significance. J Japanese Bot 63: 297–311.
  28. 28. Mittermeier RZ, Gil PR, Hoffman M, Pilgrim J, Brooks T, et al.. (2005) Hotspots revisited: Earth's biologically richest and most endangered terrestrial ecoregions. Cemex, Conservation International and Agrupación Sierra Madre, Monterrey, Mexico.
  29. 29. Zhang Q, Chiang TY, George M, Liu JQ, Abbott RJ (2005) Phylogeography of the Qinghai-Tibetan Plateau endemic Juniperus przewalskii (Cupressaceae) inferred from chloroplast DNA sequence variation. Mol Ecol 14: 3513–3524.
  30. 30. Meng L, Yang R, Abbott RJ, Miehe G, Hu T, et al. (2007) Mitochondrial and chloroplast phylogeography of Picea crassifolia (Pinaceae) in the Qinghai–Tibetan Plateau and adjacent highlands. Mol Ecol 16: 4128–4137.
  31. 31. Chen SY, Wu GL, Zhang DJ, Gao QB, Duan YZ, et al. (2008) Potential refugium on the Qinghai–Tibet Plateau revealed by the chloroplast DNA phylogeography of the alpine species Metagentiana striata (Gentianaceae). Bot J Linn Soc 157: 125–140.
  32. 32. Yang FS, Li YF, Ding X, Wang XQ (2008) Extensive population expansion of Pedicularis longiflora (Orobanchaceae) on the Qinghai–Tibetan Plateau and its correlation with Quaternary climate change. Mol Ecol 17: 5135–5145.
  33. 33. Wang LY, Abbott RJ, Zheng W, Chen P, Wang YJ, et al. (2009) History and evolution of alpine plants endemic to the Qinghai–Tibetan Plateau: Aconitum gymnandrum (Ranunculaceae). Mol Ecol 18: 709–721.
  34. 34. Wang LY, Ikeda H, Liu TL, Wang YJ, Liu JQ (2009) Repeated range expansion and glacial endurance of Potentilla glabra (Rosaceae) in the Qinghai-Tibetan Plateau. J Integr Plant Biol 51: 698–706.
  35. 35. Li C, Shimono A, Shen HH, Tang YH (2010) Phylogeography of Potentilla fruticosa, an alpine shrub on the Qinghai–Tibetan Plateau. J Plant Ecol 3: 9–15.
  36. 36. Sun YS, Ikeda H, Wang YJ, Liu JQ (2010) Phylogeography of Potentilla fruticosa (Rosaceae) in the Qinghai-Tibetan Plateau revisited: a reappraisal and new insights. Plant Ecol Diver 3: 249–257.
  37. 37. Shimono A, Ueno S, Gu S, Zhao X, Tsumura Y, et al. (2010) Range shifts of Potentilla fruticosa on the Qinghai–Tibetan Plateau during glacial and interglacial periods revealed by chloroplast DNA sequence variation. Heredity 104: 534–542.
  38. 38. Yang LC, Zhou GY, Li CL, Song WZ, Chen GC (2011) Potential refugia in Qinghai-Tibetan Plateau revealed by the phylogeographical study of Swertia tetraptera (Gentianaceae). Pol J Ecol 59: 753–764.
  39. 39. Gao QB, Zhang DJ, Duan YZ, Zhang FQ, Li YH, et al. (2012) Intraspecific divergences of Rhodiola alsia (Crassulaceae) based on plastid DNA and internal transcribed spacer of fragments. Bot J Linn Soc 168: 204–215.
  40. 40. Wang H, Qiong LQ, Sun K, Lu F, Wang YG, et al. (2010) Phylogeographic structure of Hippophae tibetana (Elaeagnaceae) highlights the highest microrefugia and the rapid uplift of the Qinghai-Tibetan Plateau. Mol Ecol 19: 2964–2979.
  41. 41. Zhang TC, Comes HP, Sun H (2011) Chloroplast phylogeography of Terminalia franchetii (Combretaceae) from the eastern Sino–Himalayan region and its correlation with historical river capture events. Mol Phylogenet Evol 60: 1–12.
  42. 42. Yang FS, Qin AL, Li YF, Wang XQ (2012) Great genetic differentiation among populations of Meconopsis integrifolia and its implication for plant speciation in the Qinghai-Tibetan Plateau. PloS One 7: e37196.
  43. 43. Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF (1998) Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol 7: 453–464.
  44. 44. Zhang YH, Volis S, Sun H (2010) Chloroplast phylogeny and phylogeography of Stellera chamaejasme on the Qinghai-Tibet Plateau and in adjacent regions. Mol Phylogenet Evol 57: 1162–1172.
  45. 45. Wang SM, Dou HS (1998) Records of Lakes in China. Beijing: Science Press.
  46. 46. Wang D (2003) The Geography of Aquatic Vascular Plants of Qinghai Xizang (Tibet) Plateau. Doctoral Dissertation of Wuhan University, Wuhan.
  47. 47. Barrett SCH, Eckert DG, Husband BC (1993) Evolutionary processes in aquatic plant populations. Aquat Bot 44: 105–145.
  48. 48. Dorken ME, Barrett SCH (2004) Chloroplast haplotype variation among monoecious and dioecious populations of Sagittaria latifolia (Alismataceae) in eastern North American. Mol Ecol 13: 2699–2707.
  49. 49. Liu X, Gituru WR, Wang QF (2004) Distribution of basic diploid and polyploid species of Isoetes in East Asia. J Biogeogr 31: 1239–1250.
  50. 50. Chen JR (2000) Flora Reipublicae Popularis Sinicae. Vol. 53 . Beijing: Science Press.
  51. 51. Doyle JJ, Doyle JL (1987) A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochemistry 19: 11–15.
  52. 52. Shaw J, Lickey EB, Beck JT, Farmer SB, Liu W, et al. (2005) The tortoise and the hare II: relative utility of 21 noncoding chloroplast DNA sequences for phylogenetic analysis. Amer J Bot 92: 142–166.
  53. 53. Taberlet P, Gielly L, Pautou G, Bouvet J (1991) Universal primers for amplification of three non-coding regions of chloroplast DNA. Plant Mol Biol 17: 1105–1109.
  54. 54. Shaw J, Lickey EB, Schilling EE, Small RL (2007) Comparison of whole chloroplast genome sequences to choose noncoding regions for phylogenetic studies in angiosperms. The tortoise and the hare III. Amer J Bot 94: 275–288.
  55. 55. Thompson JD, Higgins DG, Gibson TJ (1994) CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res 22: 4673–4680.
  56. 56. Kimura M (1980) A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol 16: 111–120.
  57. 57. Kumar S, Tamura K, Nei M (2004) MEGA3: integrated software for molecular evolutionary genetics analysis and sequence alignment. Brief Bioinform 5: 150–163.
  58. 58. Felsenstein J (1985) Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39: 783–791.
  59. 59. Nylander JAA (2004) MrModeltest v2. Program Distributed by the Author.Evolutionary Biology Centre, Uppsala University.
  60. 60. Ronquist F, Huelsenbeck JP (2003) MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574.
  61. 61. Rambaut A, Drummond AJ (2007) Tracer v1.4., Available fromhttp://beast.bio.ed.ac.uk/Tracer.
  62. 62. Excoffier L, Laval G, Schneider S (2006) ARLEQUIN version 3.1: A Software for Population Genetic Data Analysis. Computational and Molecular Population Genetics Laboratory, Institute of Zoology, University of Berne, Switzerland.
  63. 63. Clement M, Posada D, Crandall K (2000) TCS: a computer program to estimate gene genealogies. Mol Ecol 9: 1657–1660.
  64. 64. Pons O, Petit RJ (1996) Measuring and testing genetic differentiation with ordered versus unordered allels. Genetics 144: 1237–1245.
  65. 65. Dupanloup I, Schneider S, Excoffier L (2002) A simulated annealing approach to define the genetic structure of populations. Mol Ecol 11: 2571–2581.
  66. 66. Excoffier L, Smouse PE, Quattro JM (1992) Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131: 479–491.
  67. 67. Harpending HC (1994) Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol 66: 591–600.
  68. 68. Schneider S, Excoffier L (1999) Estimation of demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA. Genetics 152: 1079–1089.
  69. 69. Rogers AR (1995) Genetic evidence for a Pleistocene population explosion. Evolution 49: 608–15.
  70. 70. Rogers AR, Harpending H (1992) Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol 9: 552–569.
  71. 71. Yamane K, Yano K, Kawahara T (2006) Pattern and rate of indel evolution inferred from whole chloroplast intergenic regions in sugarcane, maize and rice. DNA Research 13: 197–204.
  72. 72. Wolfe KH, Li WH, Sharp PM (1987) Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast, and nuclear DNAs. Proc Natl Acad Sci USA 84: 9054–9058.
  73. 73. Richardson JE, Pennington RT, Pennington TD, Hollingsworth PM (2001) Rapid diversification of a species-rich genus of neotropical rain forest trees. Science 293: 2242–2245.
  74. 74. Drummond AJ, Rambaut A (2007) BEAST: bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7: 214–221.
  75. 75. Avise JC (2004) Molecular Markers, Natural History and Evolution, 2nd edn. Sunderland, MA: Sinauer Associates.
  76. 76. Hamrick JL, Godt MW (1996) Effects of life history traits on genetic diversity in plant species. Phil Trans Roy Soc London Ser B Biol Sci 351: 1291–1298.
  77. 77. Nybom H (2004) Comparison of different nuclear DNA markers for estimating intraspecific genetic diversity in plants. Mol Ecol 13: 1143–1155.
  78. 78. Shi YF, Ren JW (1990) Glacier recession and lake shrinkage indicating a climatic warming and drying trend in Central Asia. Ann Glaciol 14: 261–265.
  79. 79. Wu YQ, Cui ZJ, Liu GN, Ge DK, Yin JR, et al. (2001) Quaternary geomorphological evolution of Kunlun Pass area and uplift of the Qinghai-Xizang (Tibet) Plateau. Geomorphology 36: 203–216.
  80. 80. Owen LA, Finkel RC, Barnard PL, Ma HZ, Asahi K, et al. (2005) Climatic and topographic controls on the style and timing of Late Quaternary glaciation throughout Tibet and the Himalaya defined by 10Be cosmogenic radionuclide surface exposure dating. Quatern Sci Rev 24: 1391–1411.
  81. 81. Donnelly P, Tavare S (1986) The ages of alleles and a coalescent. Adv Appl Prob 18: 1–19.
  82. 82. Crandall KA, Templeton AR (1993) Empirical tests of some predictions from coalescent theory with applications to intraspecific phylogeny reconstruction. Genetics 134: 959–969.
  83. 83. Dynesius M, Jansson R (2000) Evolutionary consequences of changes in species' geographical distributions driven by Milankovitch climate oscillations. Proc Natl Acad Sci USA 97: 9115–9120.
  84. 84. Williams MAJ, Dunkerley DL, De Dekker P, Kershaw AP, Stokes T (1993) Quaternary Environments. London: Edward Arnold.