Introduction

Global climate fluctuated greatly during the Pleistocene, causing the distributions of most species to change markedly, including range expansion and contraction and by founder events. Such range and demographic changes undoubtedly affect the geographical patterns of genetic variation within and among populations (Hewitt, 2004). Extensive molecular-based research has helped to identify areas where species survived and also to elucidate their subsequent spread, which resulted in their present-day distributions (Avise, 2000; Petit et al., 2003; Soltis et al., 2006).

Molecular studies have clearly shown, for example, the demographic histories of arctic and alpine plants in high-latitude regions (Hewitt, 2000; Abbott and Brochmann, 2003; Schonswetter et al., 2005). These plant species survived in unglaciated peripheral or nunatak refugia during glacial periods, when massive ice sheets covered these regions, and alpine species expanded from these glacial refugia after the last glacial period.

Contrasting demographic patterns to those in high-latitude regions have been suggested for alpine plants in both temperate and tropical regions, where the cooler periods of the Pleistocene had a crucial role in their expansion (Ikeda et al., 2006; Assefa et al., 2007). During warm periods, alpine species probably retreated to a number of refugia, and they probably experienced several cycles of range contraction and expansion in response to the climatic fluctuations of the Pleistocene (Koch et al., 2006). However, genetic patterns of alpine plants in temperate and tropical regions are complex, and it is difficult to detect clear range shifts because alpine habitats are highly fragmented, having never been connected during glacial or interglacial periods. High-elevation populations may have mainly migrated upslope and downslope, thus remaining isolated throughout the palaeoclimatic cycles, with the result that there was little opportunity for inter-population gene flow (DeChaine and Martin, 2005). Moreover, stochastic demographic events, including extinction, population isolation and long-distance dispersal, can have more noticeable effects in fragmented habitats than in connected habitats, further obscuring geographic genetic patterns.

The Qinghai-Tibetan Plateau occupies an area of 2.5 million km2, with an average altitude of 4500 m. Although it is the plateau with highest altitude in the world, it was not covered by massive ice sheets during the last glacial period (Shi, 2002). Limited ice coverage during glaciations resulted in less habitat fragmentation, enabling the Plateau to serve as a refuge for many alpine plants. More than 12 000 species of plants grow on the Plateau, and it is especially rich in endemic species and genera (Lopez-Pujol et al., 2006). Therefore, it is one of the most biologically important and outstanding examples of alpine habitat. The patterns of genetic variation in present-day plant populations on the Plateau are expected to exhibit deep divergence and outstanding diversification, which should allow us to reconstruct their demographic history in more detail.

A rough demographic history of alpine plants can be deduced from palaeoenvironmental changes. Fossil pollen records suggest that the typical vegetation of the Plateau during the last glacial period was adapted to permafrost-steppe, tundra or desert habitats (Yu et al., 2000), although the temperature was 6–9 °C lower than at present (Shi, 2002). These plant communities also developed in the interior of the Plateau during the postglacial period, when forests became widespread on the eastern Plateau (Herzschuh et al., 2006). After the forest retreat, populations may have then expanded from the interior of the Plateau towards its eastern margin. Therefore, results of palaeovegetation studies indicate that the interior Plateau, with its relatively high altitudes, may have provided habitats for alpine plants throughout the glacial and interglacial periods. However, fossil pollen data from the Plateau before the last glacial period are extremely rare and insufficient for the construction of palaeovegetation maps covering the entire region, and the extent of the glacial advance on the Plateau is still under debate (Zheng et al., 2002; Lehmkuhl and Owen, 2005).

The geographical pattern of genetic variation within species of the Qinghai-Tibetan Plateau should provide insight into the above-mentioned glacial and interglacial demographic changes. Recent phylogeographic studies have examined the demographic history of alpine plants across the Plateau (Chen et al., 2008; Yang et al., 2008; Wang et al., 2009). Wang et al. (2009) suggested that Aconitum gymnandrum has survived on high-altitude parts of the central Plateau throughout the Quaternary. In contrast, Metagentiana striata and Pedicularis longiflora probably expanded from refugia on the eastern edge to the central Plateau during the interglacial and postglacial periods (Chen et al., 2008; Yang et al., 2008). Therefore, these three investigated alpine species do not share a common phylogeographical history. To understand climatic changes and past distributions of alpine plant species on the vast Plateau, further study of the extensive high-elevation area is required, for which sampling is difficult even today.

For phylogeographic study, we selected the alpine shrub Potentilla fruticosa L. (Rosaceae) as a representative of the alpine plants that currently grow throughout the Qinghai-Tibetan Plateau. P. fruticosa is a small, long-lived, deciduous shrub species that is widely distributed in alpine meadows, grasslands and forest margins. Its altitude ranges from 400 to 5000 m, according to specimen records in China (Li et al., 2003). Its flowers are self-incompatible and are pollinated not only by Hymenoptera, mainly wild bees, but also by Diptera, Lepidoptera and occasionally by Hemiptera (Guillen et al., 1995). The mature achenes are covered with long hairs, which presumably aid in their dispersal by wind (Elkington and Woodell, 1963). We examined the spatial pattern of chloroplast DNA variation in P. fruticosa populations in habitats across the northeastern and interior Plateau. Our objective was to test the following two hypotheses: (1) that the interior of the Plateau served as a refugium for P. fruticosa during interglacial periods and (2) that the range of this species expanded eastwards from the interior of the Plateau.

Materials and methods

Plant material and study sites

Leaves from P. fruticosa were collected in 2006 and 2007 from 23 sites distributed across the Qinghai-Tibetan Plateau, from the interior to the northeastern region (Figure 1). At each site, leaves were collected from spatially separated individuals; 508 samples were collected in total (Table 1). The leaf material was dried in silica gel and stored at room temperature.

Figure 1
figure 1

Map of the locations of the populations sampled in this study and the geographical distribution of cpDNA haplotypes detected in Potentilla fruticosa. cpDNA, chloroplast DNA.

Full size image
Table 1 Population number, location, coordinates, altitude, number of samples (N), number of haplotypes (nh), nucleotide diversity (Ï€) and haplotype diversity (h) with s.d. of the 23 populations of Potentilla fruticosa
Full size table

DNA extraction and sequencing

Total DNA was extracted from dried leaf material using a modified version of the hexadecyltrimethyl ammonium bromide (CTAB) method (Murray and Thompson, 1980).

Part of the matK gene of the chloroplast DNA was amplified by PCR using the AF and 8R primers (Ooi et al., 1995). PCR was performed in a total volume of 15 μl, containing 1 × PCR buffer, 2.5 mM MgCl2, 0.2 μM of each primer, 0.2 mM of each dNTP, 0.25 U Ampli Taq Gold (Applied Biosystems, Foster City, CA, USA) and ∼20 ng of template DNA. PCR products were cleaned using the QIAquick PCR purification Kit (Qiagen, Tokyo, Japan).

Forward and reverse strands were cycle-sequenced with the same PCR primers, using a BigDye Terminator version 3.1 Cycle Sequencing Kit (Applied Biosystems). For the cycle-sequencing reaction, the thermal profile was 25 cycles at 96 °C for 10 s, 50 °C for 5 s and at 60 °C for 4 min. Products were analysed on an ABI 3730 automated sequencer (Applied Biosystems). Forward and reverse sequences were assembled using Vector NTI software (Invitrogen, Carlsbad, CA, USA). Sequences were compared visually with the original chromatograms to avoid reading errors. All sequences are accessible at GenBank (accession nos. AB458562–AB458594).

Sequence analysis

The DNA sequences were aligned using the CLUSTAL_X 1.81 software (Thompson et al., 1997) and refined manually. Indels were treated as a single character resulting from one mutation event.

Relationships among haplotypes were resolved by constructing a statistical parsimony network and a maximum parsimony tree. The statistical parsimony network was produced using the TCS version 1.13 software (Clement et al., 2000) with a 95% statistical confidence limit, as outlined by Templeton et al. (1992). Maximum parsimony trees were constructed using the PAUP 4.0b10 software (Swofford, 2003), applying a tree bisection and reconnection branch-swapping algorithm. The consensus tree was estimated using a bootstrap approach with 1000 replicates. To determine the ancestral haplotype, Fragaria vesca was selected as an outgroup by referring to the phylogeny of the Rosaceae (Eriksson et al., 2003). Eriksson et al. (2003) reported that the clade including P. fruticosa comprises a number of Potentilla species along with several other groups (namely Fragaria, Alchemilla, Chamaerhodos and Sibbaldia) and that F. vesca is the sister group of P. fruticosa. The out-group sequence for matK (Potter et al., 2002) was downloaded from GenBank (accession no. AF288102). We also tested another related species (Potentilla anserine, accession no. AF288113) as an outgroup, to confirm whether we would obtain the same topology.

To measure levels of genetic variation, nucleotide and haplotype diversities (Nei, 1987) were calculated for each population using ARLEQUIN version 2.0 software (Schneider et al., 2000).

An AMOVA (analysis of molecular variance) (Excoffier et al., 1992) was also conducted using ARLEQUIN version 2.0 (Schneider et al., 2000). The F-statistic (Fst) was calculated, and significance was tested by using 10 000 permutations.

Spatial genetic analyses

We visualized detailed patterns of spatial genetic structure among populations by genetic landscape shape interpolation using the AIS (Alleles in Space) software (Miller, 2005). This procedure is initiated by constructing a Delauney triangulation-based connectivity network of sampling areas, and it creates three-dimensional genetic landscapes by interpolation of measured genetic distances between sampled individuals. As AIS interpolates using residuals from an isolation-by-distance regression across individuals, unusually high genetic distances between individuals within an area result in high elevations on the structural map. We selected a grid size of 50 × 46 to best represent the rectangular shape of the study area and a distance-weighting parameter of 1. We also tested other grid sizes and a range of distance-weighting parameters (0.5–2) to confirm that qualitatively similar results would be obtained.

Demographic analyses

To test the hypothesis that the range of P. fruticosa expanded eastwards, we divided the populations into two regions, namely the northeast (sites 1–17) and the interior (sites 18–23), and analysed the demographic change according to a sudden expansion model by examining the mismatch distribution (the distribution of pairwise differences between haplotypes). This model assumes a sudden population growth from N0 to N1 that occurred t generations ago, after which the population reaches demographic equilibrium (Rogers and Harpending, 1992). The mismatch distribution is usually multimodal in samples collected from populations at demographic equilibrium that have been of relatively stable size over time, but it is unimodal in populations that have experienced a recent demographic expansion (Rogers and Harpending, 1992). We used ARLEQUIN version 2.0 software for this implementation (Schneider et al., 2000). This software estimates three demographic parameters by a generalized nonlinear least-squares approach (Schneider and Excoffier, 1999), θ0=2uN0, θ1=2uN1 and τ=2ut, where u is the mutation rate. The validity of this stepwise expansion model is tested by a parametric bootstrap approach simulating sum of squared deviations, which is also used to estimate confidence interval values of the estimated parameters.

Results

Geographical distribution and phylogenetic analysis of haplotypes

We analysed a total of 508 individuals from 23 populations of P. fruticosa on the Qinghai-Tibetan Plateau and detected 33 haplotypes by examining ∼1120 bp of chloroplast DNA from the matK region. Variable sites among the 33 haplotypes were observed as 56 nucleotide substitutions and 2 indels, of which 18 were parsimony informative (Table 2). Haplotype A was dominant and widely distributed over the northeastern Plateau (Figure 1). Haplotypes A1 and B1 were also detected frequently in the northeastern and interior Plateau, respectively, whereas the other haplotypes tended to be restricted to local areas. Haplotype B especially occurred with low frequency in two disjunct areas, that is, in site 4 in the northeast and in sites 19 and 21 in the interior region.

Table 2 Variable sites of the aligned sequences of the chloroplast DNA fragment matK region in the haplotypes of Potentilla fruticosa
Full size table

Parsimony analysis produced 12 most-parsimonious trees, which required 80 steps (consistency index=0.950; retention index=0.944). The topology of the strict consensus tree of the 12 most-parsimonious trees is shown in Figure 2. The same topology was obtained using another outgroup, P. anserine. The haplotype network drawn by the statistical parsimony method also clearly supports this topology. Three haplogroups are apparent in both the parsimony tree and the haplotype network. The relationships among the haplotypes within each haplogroup are unresolved polytomies; that is, each haplotype network exhibits a star-shaped phylogeny. Geographically, the basal haplotype I and haplogroup 1, which consists of haplotypes C–H, were confined to the interior region (sites 18–22). Haplogroup 2, consisting of haplotype B and haplotypes derived from haplotype B, was widespread throughout the Plateau, but occurred with low frequency in the northeastern region. Haplogroup 3, consisting of haplotype A and the haplotypes derived from haplotype A, was dominant and widely distributed over the northeastern region.

Figure 2
figure 2

Phylogenetic relationships inferred among the cpDNA haplotypes of Potentilla fruticosa based on the sequences of the cpDNA matK region. Left: strict consensus tree of the 12 most-parsimonious trees. The numbers on the branches are bootstrap support values for 1000 replicates. Branches corresponding to partitions reproduced in <50% of the bootstrap replicates are collapsed. Right: the haplotype networks. The sizes of circles correspond to the frequency of each haplotype. The small black circles on the branches represent the inferred intermediate haplotypes not observed in the samples. Grey line and dots indicate the network that was not supported by the phylogenetic tree. cpDNA, chloroplast DNA.

Full size image

Genetic diversity and spatial genetic structure

Most interior populations (sites 18–23), which had widely divergent haplotypes (haplotypes C–I), exhibited high levels of nucleotide diversity (Table 1). In contrast, northeastern populations (sites 1–17) exhibited relatively low levels of nucleotide diversity, a tendency that was particularly noteworthy in populations occurring at altitudes below 4000 m (Table 1).

The AMOVA showed that ∼48.7% of the genetic variation occurred among populations, whereas ∼51.3% of the variation occurred within populations in the study area. The Fst value was 0.487 (P<0.0001).

Genetic landscape shape interpolation analysis produced a clear surface plot, showing that genetic distances increased towards the interior region (Figure 3a). Four high ridges were observed between site 19 and the adjacent sampling sites (Figure 3b).

Figure 3
figure 3

(a) Genetic landscape shape interpolation analysis using a 46 × 50 grid. X and Y axes correspond to geographical locations within the overall physical landscape examined in this study (Figure 1). Surface plot height (Z axis) reflects residual genetic distance. (b) Delaunay triangulation (fine lines) and the primary barriers (thick line). Dots indicate the sampling points. The greatest genetic barrier was detected between site 19 and its neighbours.

Full size image

Demographic analyses

The distribution of pairwise differences between plants from the northeastern populations (sites 1–17) was unimodal, indicating that the sudden expansion model is applicable (Figure 4a). The main expansion event was estimated to have occurred at τ=0.8 (95% confidence interval, 0–2.1). The newly derived haplotypes (haplogroup 3) were mainly found in these populations.

Figure 4
figure 4

Mismatch distribution analysis for the (a) northeastern and (b) interior populations. Histograms correspond to observed frequencies of pairwise nucleotide differences, and lines represent the expected frequencies under the sudden expansion model. Goodness-of-fit of the observed versus the theoretical mismatch distributions under a sudden expansion model was tested using the sum of squared deviations (SDD). Upper and lower 95% confidence limits around estimates of Ï„ are given in parentheses.

Full size image

The pairwise differences in the interior populations (sites 18–23), in contrast, exhibited a bimodal distribution, although this result did not lead us to reject the sudden expansion model (Figure 4b). The main expansion event was estimated to have occurred at τ=13.2 (95% confidence interval, 6.3–23.7). The estimated expansion time in the interior region was thus 17 times that in the northeastern region. These estimates suggest that interior populations existed for a long period of time before they expanded towards the northeast.

Discussion

The interior as a refugium on the Qinghai-Tibetan Plateau

Refugia, where species are able to survive for a long time, generally represent long-term reservoirs of the genetic variation of a species. In these locations, evolution often produces unique genotypes and high levels of diversity (Comes and Kadereit, 1998; Hewitt, 2000). We identified the interior as one such location on the Plateau for P. fruticosa. The following three pieces of evidence led to this conclusion: (1) the interior populations (sites 18–22) had widely divergent haplotypes (Figure 2); (2) these populations also exhibited high levels of genetic diversity (Table 1) and (3) within these populations, unique haplotypes were generally dominant and restricted to small areas (Figure 1). The high levels of diversity and unique genotypes in the interior region may also explain the clear tendency of genetic distances to increase towards the interior in the genetic landscape shape interpolation analysis (Figure 3). In particular, population 19 seems to be the longest surviving population in the sampling area of this study, because (1) it shows the highest genetic diversity among all populations (Table 1), (2) it harbours the ancestral haplotype I (Figures 1 and 2), and (3) clear genetic barriers exist between it and adjacent populations (Figure 3).

Palaeoenvironmental evidence also supports our view that the interior region may have served as a refugium for alpine species. The fossil pollen record shows that, during the mid-Holocene, when temperatures exceeded present-day temperatures by at least 1–2 °C (Yu et al., 2000), forests were located at markedly higher elevations than they are today (Herzschuh et al., 2006; Ren, 2007). This phenomenon had probably occurred during each interglacial period. Therefore, at these times, alpine species in the eastern region must have withdrawn to unforested areas of relatively high altitude, such as in the interior region, where altitudes exceed 4500 m.

The interior region probably provided a suitable habitat for P. fruticosa during glacial periods as well, because the Plateau was apparently not covered by massive ice sheets even during the last glacial maximum (Zheng et al., 2002). During the last glacial maximum, glaciers are estimated to have descended to ∼3500 m in the northeast and to 5000 m in the interior of the Plateau, reflecting the decreasing precipitation from east to west (Shi, 2002). These altitudes are more than 1000 m higher than those in the other regions at the same latitude (for example, Japan), partly as a result of the low precipitation in winter on the Plateau.

Demographic history of P. fruticosa

Past demographic and distributional changes of the alpine species are expected to be reflected in the genetic variations within and between populations. Our data strongly suggest that populations of P. fruticosa experienced drastic range expansions in response to past climatic fluctuations. The current northeastern populations apparently underwent a recent expansion from a restricted source area. This conclusion is supported by the unimodal mismatch distribution of the northeastern populations (Figure 4a), suggesting historical population expansion, and by the star-shaped phylogeny of haplotypes derived from haplotype A (Figure 2), the reduced levels of haplotypic and nucleotide diversity (Table 1), and the prevalence of haplotype A in the northeastern populations (Figure 1). The homogeneous genetic pattern in the northeastern populations is consistent with the AMOVA result, which showed relatively low genetic variation among populations (48.7%) in comparison with other alpine plant species on the Plateau (Chen et al., 2008; Yang et al., 2008; Wang et al., 2009). These patterns are consistent with the results of a simulation study in which sudden expansion from a few populations produced reduced allele diversity and areas of genetic homogeneity (Avise, 2000).

This recent expansion most likely occurred after the beginning of the Holocene, as suggested by the fossil pollen record for the region. For example, a site at 3400 m on the eastern margin was covered by a cool mixed forest at 6 ka B.P. (Yu et al., 2000). After 6 ka B.P., the forest cover decreased and alpine meadow became widely distributed in this region (Ren, 2007). At present, coniferous forest is located between 2000 and 2900 m.

Phylogeographic studies suggest a similar history of coniferous trees (Zhang et al., 2005; Meng et al., 2007), the dominant species in the forested zone and which occasionally grow within the alpine zone on the northeastern Plateau. These scattered conifer populations probably reflect colonization from a refugial population that existed in the past along the Plateau's edge. In addition, results of a phylogeographic study of the alpine herb A. gymnandrum, which is widely distributed in alpine meadows across the Plateau, suggest the existence of refugia at high altitudes on the central Plateau (Wang et al., 2009). Conversely, refugia of another two alpine species were inferred to have been located at the Plateau's edge (Chen et al., 2008; Yang et al., 2008).

Further comparative studies with other regions and other species are required to confirm the history of alpine species on the Qinghai-Tibetan Plateau. Extensive studies conducted in high-latitude regions, especially in Europe, allow us to extract shared phenomena among species, although each species has its own unique history during glacial and interglacial periods (Schonswetter et al., 2005).

We also inferred that P. fruticosa populations probably experienced a demographic expansion before the recent expansion described above, and haplotype B prevailed across the Plateau. This scenario is supported by the star-shaped phylogeny exhibited by haplotypes derived from haplotype B (Figure 2) and by the occurrence of haplotype B in two disjunct areas, that is, in site 4 in the northeast and in sites 19 and 21 in the interior region (Figure 1). After the expansion, the P. fruticosa population would have contracted to the interior of the Plateau in response to climatic fluctuations, with the result that lineages derived from haplotype B would have been distributed mainly in interior regions, with only a few examples left in the northeast region (Figure 1). However, further studies are required to confirm this early range shift of the P. fruticosa population.

Conclusions

We reconstructed the demographic history of P. fruticosa on the Plateau by using genetic variation data obtained from spatially distributed populations. Our results suggest that periods of population expansion were associated with climatic cooling, and that during climatic warming, populations contracted to the interior of the Qinghai-Tibetan Plateau. Thus, the interior region acted as a refugium during interglacial periods, and greatly contributed to the diversification of P. fruticosa.