Skip to main content

Phylogeny and biogeography of Primula sect. Armerina: implications for plant evolution under climate change and the uplift of the Qinghai-Tibet Plateau



The historical orogenesis and associated climatic changes of mountain areas have been suggested to partly account for the occurrence of high levels of biodiversity and endemism. However, their effects on dispersal, differentiation and evolution of many groups of plants are still unknown. In this study, we examined the detailed diversification history of Primula sect. Armerina, and used biogeographic analysis and macro-evolutionary modeling to investigate a series of different questions concerning the evolution of the geographical and ecological distribution of the species in this section.


We sequenced five chloroplast and one nuclear genes for species of Primula sect. Armerina. Neither chloroplast nor nuclear trees support the monophyly of the section. The major incongruences between the two trees occur among closely related species and may be explained by hybridization. Our dating analyses based on the chloroplast dataset suggest that this section began to diverge from its relatives around 3.55 million years ago, largely coinciding with the last major uplift of the Qinghai-Tibet Plateau (QTP). Biogeographic analysis supports the origin of the section in the Himalayan Mountains and dispersal from the Himalayas to Northeastern QTP, Western QTP and Hengduan Mountains. Furthermore, evolutionary models of ecological niches show that the two P. fasciculata clades have significantly different climatic niche optima and rates of niche evolution, indicating niche evolution under climatic changes and further providing evidence for explaining their biogeographic patterns.


Our results support the hypothesis that geologic and climatic events play important roles in driving biological diversification of organisms in the QTP area. The Pliocene uplift of the QTP and following climatic changes most likely promoted both the inter- and intraspecific divergence of Primula sect. Armerina. This study also illustrates how niche evolution under climatic changes influences biogeographic patterns.


Understanding the processes that shape geographical and ecological distribution of biodiversity is one of the most challenging questions in evolutionary biology and ecology. This is particularly true for regions that have experienced rapid habitat changes and harbor high species diversity. These characteristics are present in many mountainous areas and historical orogenesis has been proposed to play an important role in shaping their current biodiversity [13]. The alteration of topography and climatic changes associated with mountain uplifts can cause fragmentation of species distributions, thus limiting gene flow between isolated populations and initiating allopatric divergence and speciation [47]. However, extreme environmental changes and fragmented distributions can also lead to the extinction of lineages and species (e.g., [8, 9]). The processes occurring during mountain uplifts are therefore complex and we need to better understand the mechanisms that are at play during these events.

The fragmentation of species distributions can be due to the presence of limits on dispersal due, for example, to geographical barriers. Such limitations can induce a reduction in the movement of individuals into new locations and will result in distinct biogeographic patterns in the extant species [10]. However, fragmentation can also occur because of a lower success of establishment of individuals in some areas, which will limit the range of species [11]. This process is primarily set by ecological factors, potentially including both abiotic and biotic variables [1012]. The dynamics of species range evolution will be constrained by phylogenetic niche conservatism, which is defined as the tendency of species to retain their ancestral ecological niche, thus shaping the geographic ranges of species over time (e.g., [13, 14]). However, evidence for rapid shifts in climatic preferences among species also exists [15, 16] and macro-evolutionary modeling should be used to characterize the processes driving the evolution of ecological niches [17]. A complete assessment of these processes, coupled with detailed analyses of biogeographic patterns of species distribution, should then be used to help understand the distribution of species diversity [10].

One region that experienced drastic habitat changes and harbors extremely rich species diversity and endemism is the Qinghai-Tibet Plateau (QTP; [18]). While the start of its uplift dates from approximately 50 million years ago (Ma; [19]), the extensive uplifts of the QTP occurred in at least four periods since the early Miocene, specifically between 25–17 Ma, 15–13 Ma, 8–7 Ma, and 3.4-1.6 Ma [9, 2023]. At present, the QTP, with an average altitude of more than 4000 m (a.s.l.), is the highest and one of the most extensive plateaus on Earth [20]. About 9,000–12,000 species of vascular plants in ca. 1,500 genera are present in this plateau, and at least 20 % of these species and ca. 50 genera are endemic [3, 18]. The historical sequence of uplifts of the QTP has been suggested to partly account for the occurrence of high levels of biodiversity and endemism in the region [24]. However, the potential effects of climatic changes during the Quaternary on the diversification and distribution of many groups of plant species in the QTP are not very well known (see Review [2, 3, 25]).

Primula L. (Primulaceae) is one of the genera that exhibit high levels of species diversity in the QTP. The group, with a predominantly northern hemisphere distribution, contains ca. 500 species. About 60 % of the species are present in the QTP and its adjacent regions [26, 27]. Although this genus represents an important floristic element of alpine meadows in the region, it remains unclear whether the uplift of the QTP and the following climatic changes affected its diversification and distribution. In this context, a better understanding of the historical biogeography of key floristic elements of the region is an important way to illuminate the evolutionary history of these organisms in space and time. Available studies mainly utilize genus- or family-level phylogenies to elucidate the biogeographic connections between the QTP and neighboring regions [2832]. However, the presence of a single sample per species hardly provides insights into the biogeographic patterns of species distributions within the QTP. Therefore, sampling multiple individuals per species and focusing on endemic species may help to better understand the mechanisms that were responsible for biogeographic patterns within the QTP.

In this study, we include several samples per species to investigate the historical biogeography of Primula sect. Armerina Lindley (Primulaceae), which exhibits a typical Sino-Himalayas distribution. According to the most recent global monographic treatment of the genus, Primula sect. Armerina comprises 14 species [26]. Eight species (P. fasciculata, P. tibetica, P. conspersa, P. gemmifera, P. zambalensis, P. pumilio. P. pamirica and P. involucrata; Fig. 1) are endemic to the QTP, with different geographic distributions [26, 27]. Among them, there has been some confusion between P. tibetica and P. fasciculata because of their morphological similarities at high altitude ([26, 27]; field observation). The two species can be easily distinguished when bracts are present. Primula tibetica has oblong and pouched bracts, while the bracts of P. fasciculata are linear and non-pouched (Fig. 1a, d). However, at high altitude, bracts are usually missing in P. fasciculata (Fig. 1b, c), while in P. tibetica, they can also be absent in small individuals with single flower (Fig. 1e, f). Both species have wide altitude distributions, ranging from 2900 m to 5000 m [26, 27] and the use of molecular data combined with macro-evolutionary modeling may provide useful insights into the dynamics of their range evolution. The four remaining species of this section (P. iljinskyii, P. chrysostoma, P. knorringiana and P. valentinae) have very restricted areas in regions adjacent to the QTP. Primula nutans has the most widespread distribution in the section, including N Europe, W & E Siberia, NW America to N Mongolia, NW China and NW QTP. All species from sect. Armerina are considered to be diploid (2n = 18, 20 or 22) [26, 27], except P. egaliksensis, which is the only tetraploid species (2n = 36, 40) and occurs mainly in North America. It was assigned to sect. Armerina based on morphological features [33, 34], and might be of hybrid origin between P. mistassinica (sect. Aleuritia) and P. nutans [3537].

Fig. 1

The five species of sect. Armerina which showed mainly incongruence between the two trees. (a) P. fasciculata with linear and non-pouched bracts, (b) P. fasciculata without bracts, (c) one photo of P. fasciculata collected from populations of clade F2 (see Results), (d) P. tibetica with oblong and pouched bracts at low altitude, (e) and (f) P. tibetica with and without bracts at high altitude, respectively, (g) P. nutans, (h) P. gemmifera, (i) P. conspersa. Bracts for P. fasciculata and P. tibetica are indicated by red arrows. All photos were taken by the first author in the field

Most species of the Armerina section are thus prominent floristic elements of alpine meadows at high altitudes in the QTP and most are endemic to the QTP and its adjacent regions [26, 27]. This section of Primula hence represents a good candidate to assess the biogeographic history of the QTP and to understand the effects of its uplift and associated climatic changes on the geographical distribution of biodiversity. We analyzed both nuclear and chloroplast DNA sequences of multiple samples per species in the Armerina section to reconstruct a comprehensive phylogenetic tree of this group. The aims of our study are to: i) test the inter-specific relationships of sect. Armerina to obtain a detailed and resolved phylogenetic tree for the section; ii) assess whether the diversification of this section was influenced by the uplifts of the QTP; iii) combine biogeographic analyses with macro-evolutionary modeling of ecological niches to better understand the role of dispersal and ecological constraints during the diversification of the three main species in the section (P. fasciculata, P. nutans and P. tibetica).


Sequence characteristics

Five chloroplast (matK, rpl16, rps16, trnLF and trnH-psbA) and one nuclear (translin family protein, tfp) markers were sequenced in this study for phylogenetic analyses. The matK dataset comprised 892 characters, 815 of which were constant, 22 variable but parsimony-uninformative, 55 variable and parsimony-informative. The rpl16 dataset comprised 1063 characters, 903 of which were constant, 90 variable but parsimony-uninformative, 70 variable and parsimony-informative. The rps16 dataset comprised 877 characters, 789 of which were constant, 24 variable but parsimony-uninformative, 64 variable and parsimony-informative. The trnLF dataset comprised 968 characters, 840 of which were constant, 54 variable but parsimony-uninformative, 74 variable and parsimony-informative. The trnH-psbA dataset comprised 629 characters, 512 of which were constant, 46 variable but parsimony-uninformative, 71 variable and parsimony-informative. We combined the five plastid regions for all subsequent analyses, modeling them as five partitions. It was not possible to obtain these sequences for P. watsonii and four chloroplast sequences (matK-DQ378314, rpl16-DQ378443, rps16-FJ786584 and trnLF-FJ794215) were downloaded from GenBank for this species.

The aligned nuclear dataset comprised 648 characters, 445 of which were constant, 91 variable but parsimony-uninformative, and 112 variable and parsimony-informative. Despite repeated attempts, the tfp sequences for three samples of P. tibetica, as well as the sample of P. pamirica, P. pumilio and two outgroup species (P. watsonii and P. pinnatifida) failed to amplify. Two copies were identified in the samples of P. fasciculata, P. conspersa and P. egaliksensis and these clones were added to the sequences obtained directly from PCR in subsequent phylogenetic analyses.

Phylogenetic analyses and molecular dating

The maximum likelihood (ML) and Bayesian analyses done on each data set resulted in congruent topologies, but discrepancies were obtained between the two types of markers. The only tetraploid species, P. egaliksensis, was included in a well-supported clade with P. mistassinica and P. farinosa in the chloroplast tree. This result is in agreement with previous studies [35, 37, 38]. The node subtending the rest of the samples of Primula sect. Armerina received very low support (posterior probability, PP 0.18, ML 6 %) in the choloroplast phylogenetic tree and the relationships between species remained partly unresolved (Fig. 2). Three main clades were inferred in the chloroplast tree. The clade involucrata (including P. involucrata, P. pamirica, P. fasciculata, P. nutans and P. tibetica) and the clade conspersa (including P. conspersa, P. gemmifera and P. zambalensis) were strongly supported in both ML and Bayesian analyses, while the clade pumilio (P. pumilio) was not well-supported by ML (74 %), but received very high posterior probabilities in the Bayesian analyses (PP 1.0). Overall, well-supported clades (PP > 0.95) in the chloroplast tree grouped sequences from the same species, except for P. fasciculata, which was separated into two groups (Fig. 2).

Fig. 2

The maximum clade credibility (MCC) tree derived from BEAST analyses of five chloroplast genes. Maximum likelihood (ML) bootstrap values and Bayesian posterior probabilities (PP) are indicated at major nodes. Bootstrap values ≥ 80 and PP ≥ 0.95 are indicated with thicker branches. Outgroup species are shown in bold

In contrast to the plastid dataset, Primula sect. Armerina and two nested outgroup species received very high node support (PP 1.0, ML 100 %) in the nrDNA phylogenetic tree, but the relationships between species were less well supported (Fig. 3). Three main clades within the section identified in the chloroplast tree were also inferred in the nuclear tree (Fig. 3). The clade involucrata was well-supported (PP 1.0, ML 86 %), while the clades conspersa (except for P. farinosa, P. mistassinica and P. egaliksensis) and pumilio received very weak nodal support in both types of analyses. The relationships within each clade were further incongruent between the trees obtained by the two datasets. Primula fasciculata was divided into three clades in the nrDNA tree (Fig. 3). One clade included samples from P. fasciculata that cluster with a moderately supported clade representing P. involucrata. A second clade included all samples of P. tibetica and P. fasciculata and one copy of P. fasciculata. Finally, the third clade included all samples of P. nutans and P. pamirica, one copy of P. egaliksensis and the remaining samples of P. fasciculata (Fig. 3). Similarly, P. gemmifera separated into two groups, either with P. zambalensis or in a clade including all samples of P. conspersa (Fig. 3). Two copies of P. egaliksensis were clustered with either P. nutans or P. mistassinica, corroborating the hypothesis of the allopolyploid origin of this species [3537].

Fig. 3

The maximum clade credibility (MCC) tree derived from MrBayes analyses of the nuclear dataset. Maximum likelihood (ML) bootstrap values and Bayesian posterior probabilities (PP) are indicated at major nodes. Bootstrap values ≥ 80 and PP ≥ 0.95 are indicated with thicker branches. Outgroup species are shown in bold. Two nuclear gene copies for some samples are indicated with “-1” or “-2”

Previous dating analyses at the level of the family used low intra-sectional sampling and suggested that sect. Armerina diverged from its relatives about 5 Ma [38]. This date is generally congruent with the results of our dating analysis, which indicated that the section (except P. egaliksensis) diverged from its two relatives, P. watsonii and P. pinnatifida, 3.55 Ma (1.76–5.93 Ma, 95 % highest probability density, HPD; Fig. 3). Most cladogenetic events in this section occurred during the past 3.4 million years (Fig. 4). The crown age of the three closely related species, P. nutans, P. fasciculata and P. tibetica, was about 1.19 Ma (95 % HPD: 0.51–2.13 Ma; Fig. 4).

Fig. 4

Dispersal–vicariance scenarios for sect. Armerina and the outgroup speices based on the chloroplast dataset reconstructed by Statistical Dispersal–Vicariance Analysis (S-DIVA) optimization with the maximum number of area units set to two. Triangle: dispersal event; diamond: vicariance event. Letters denoting area units are indicated on the map. Pie charts at internal nodes represent the marginal probabilities for each alternative ancestral area. Alternative ancestral areas (letters on nodes) are indicated for the major nodes. The grey bars on the nodes represent the 95 % highest posterior density intervals of the dates obtained from BEAST analyses. Time scale is shown at the bottom. Three groups (F1, F2 and NT) are used for the evolutionary niche models: groups F1 and F2 are two clades of P. fasciculata in the chloroplast tree; group NT includes all samples of P. tibetica and samples of P. nutans that were only collected from the QTP

Biogeographic inference

Biogeographic analysis based on the chloroplast dataset was reconstructed by Statistical Dispersal–Vicariance Analysis (S-DIVA). Fourteen dispersal and 15 vicariance events for the section were identified in this analysis (Fig. 4). The origin of this section was inferred with high confidence in the Himalayan Mountains (B, 91 %). We found that one clade (P. zambalensis, P. gemmifera and P. conspersa) colonized the Northeast QTP (C) and subsequently diversified and dispersed to the Hengduan Mountains (A), while P. pamirica colonized the Mountains of Central Asia (D). The common ancestral area of P. fasciculata, P. tibetica and P. nutans was inferred to be in the Himalayan Mountains (B, 86 %).

Evolution of ecological preferences

We fitted a series of macro-evolutionary models based on 19 bioclimatic variables (i.e., climatic niches) to better understand the biogeographic patterns of three closely related species, P. fasciculata, P. tibetica and P. nutans. We extracted the 19 bioclimatic variables from the sampled localities of the three species (Additional file 1). For P. nutans, we used only the samples that were collected in the QTP. The first two axes of the principal component (PC) analysis based on this dataset explained 53.2 % and 25.3 % of variance, respectively. The first axis (PC1) was strongly and positively correlated with temperature seasonality (BIO4, WorldClim variables) and negatively correlated with temperature in coldest and driest Quarter (BIO6 and BIO9). The second axis (PC2) was correlated strongly and positively with precipitation in coldest and driest Quarter (BIO14, BIO17 and BIO19), and strongly and negatively with precipitation seasonality and mean diurnal range (BIO2 and BIO15).

We used the values obtained for PC1 and PC2 (Additional file 2) to test for the evolution of the ecological niche in P. fasciculata, P. tibetica and P. nutans. The Brownian motion model was rejected for both PC1 and PC2 in all species sets tested (Additional file 3). For PC1, the best-performing models were OU1 for SET1, SET2 and SET3, and OUMV for SET4. Average AICc weights were 0.46, 0.36, 0.66 and 0.48, respectively (see Additional file 3 for all AICc weights). The OUM was the second-best model for SET1 (Average AICc weights = 0.25). The OUMV, OUMA and OUM models that allow different niche optima for SET2 also received non-negligible AICc weights (0.29, 0.18, 0.12). For PC2, all four sets were best modeled under OUMV (AICc weights 0.97, 0.93, 0.78 and 0.64 respectively; Additional file 3).

The parameters (niche optimum θ, rate of niche evolution σ2 and strength of selection α) estimated for the three species groups (F1, F2 and NT) from all supported models based on the four group sets were congruent (Additional file 4) and we showed the parameters estimated based on SET2 (Fig. 5). We used model averaging to estimate the parameter values for PC1 over the supported models OUMV, OUMA and OUM. The averaged niche optima (θ) across models for group F1, F2 and NT were −0.17, −2.0 and 0.55, respectively (Fig. 5). The averaged rate parameter (σ2) across models for group F2 was two times slower than that for the groups F1 and NT (59 vs. 131 and 112). Finally, the averaged strength of selection estimated across models for the three groups was similar (6.9, 6.3, 6.9). For PC2, model OUMV, which allows for different niche optima and rates of niche evolution among groups, was the only supported model. The optimum values estimated based on this model for the three groups were also different from each other (F1: 0.2, F2: −0.99, NT: −0.33). The group F2 still exhibited the slowest rate of niche evolution (F1: 228, F2: 94, NT: 1723; Fig. 5).

Fig. 5

Parameter estimates of models of niche evolution for the three groups (F1, F2 and NT). For PC1, averaged parameters are obtained based on three supported models (OUM, OUMV and OUMA). The averaged strength of selection (α) estimated across models for the three groups is similar and not shown. For PC2, parameter estimates are from the only supported OUMV model (different rates σ2 and niche optima θ among the three groups)


Non-monophyly of Primula sect. Armerina

The phylogenetic analyses of Primula sect. Armerina presented here contain samples of several individuals per species and cover most of the geographic distributions of the species. Neither the chloroplast tree nor the nuclear tree supports the monophyly of sect. Armerina. The section and the two outgroup species, P. watsonii and P. pinnatifida (sect. Muscarioides) form a well-supported clade in the chloroplast tree, despite the fact that these two outgroup species are distinguished from sect. Armerina by clear morphological traits (e.g., spicate inflorescence vs. umbel; [26, 27]). Similarly, the two outgroup species, P. farinosa and P. mistassinica (sect. Aleuritia), are grouped with the section and form a well-supported clade in the nrDNA tree. The non-monophyly of sect. Armerina is in agreement with previous family-wide analyses [38, 39]. Moreover, non-monophyly of sections in genus Primula seems pervasive in phylogenetic trees [38, 39].

Phylogenetic relationships within the section

The relationships among some of the basal nodes of the section in the nuclear tree are uncertain (Fig. 3), which may result from low sequence divergence within the section. The use of a single nuclear gene is thus clearly not sufficient to resolve the relationships within the group, which is a pattern often found also in other lineages (e.g., [40, 41]). Multiple nuclear genes or genomic data are therefore needed to resolve the precise relationships between the main clades in this group. However, both phylogenetic trees show three main clades within sect. Armerina, which is in agreement with previous phylogenetic studies [39] as well as morphological based taxonomy [26, 27].

Phylogenetic relationships inferred from the nuclear and chloroplast datasets were incongruent (Figs. 2, 3). The tree obtained from the latter is in agreement with morphology-based taxonomy, which contrasts with other studies that showed a better congruence of taxonomy with the trees inferred from nuclear datasets (e.g., [42]). Incongruence between different plant genomic markers is found in numerous studies and can be explained by incomplete lineage sorting, hybridization and introgression [40, 4245]. Introgression represents the transfer of genes between species mediated primarily by backcrossing [46], but it does not seem a likely explanation for the incongruence that we observed. Maternally inherited chloroplast loci with relatively low rates of intraspecific gene flow should be more frequently introgressed [46]. In contrast, biparentally inherited nuclear loci that experience high rates of intraspecific gene flow should enhance species delimitation [46]. We find the opposite pattern in our results (Figs. 2 and 3). The chloroplast tree has much clearer species delimitation than the nuclear tree and this pattern seems incompatible with the assumption that the incongruence results from introgression.

Although introgression cannot occur without hybridization, hybridization followed by no backcrossing and introgression could still occur and such phenomenon has been detected in numerous studies (e.g., [47, 48]). Natural hybridization in Primula is common and has been confirmed by several studies [37, 4951], although, it is currently unclear to what degree species within sect. Armerina hybridize with each other. The incongruent placement of P. egaliksensis between chloroplast and nuclear gene trees can be explained by hybridization [3537]. Moreover, our results provide further evidence in support of the hypothesis that P. egaliksensis originated from an intersectional allopolyploidization event, which places the two tfp copies of P. egaliksensis with either P. nutans or P. mistassinica (Fig. 3), confirming previous results by Guggisberg et al. [3537]. From this perspective, similar incongruence detected for sample number 14 of P. fasciculata (two tfp copies grouped with either P. nutans or P. tibetica) may also result from hybridization. Beside hybridization, incomplete lineage sorting is another important explanation for the incongruence between data sets, but the two processes are often difficult to distinguish from each other [5254]. Although incomplete lineage sorting could also be involved in the incongruences found in our results, the occurrence of such a process would imply that the origin of the haplotypes of P. pamirica preceded the speciation events of the whole clade [52]. Such extensive levels of incomplete lineage sorting may yield gene trees with random patterns of relationships among taxa [55]. The patterns of relationships are, however, non-random in the nuclear tree. The major incongruences result mainly from the division of P. fasciculata and P. gemmifera in different lineages. We thus consider hybridization as the most likely explanation for the major incongruences between chloroplast and nuclear trees.

However, it should be noted that using a single nuclear gene that provides low resolution of the phylogenetic relationships might not be sufficient to elucidate the reasons of the genealogical incongruences between different genomic markers. Although our results tend to suggest a more probable role of hybridization as the most likely explanation for the major incongruence between two trees, incomplete lineage sorting and introgression cannot be completely excluded. We therefore recognize that gene trees/species trees analyses involving multiple nuclear loci or population genomic approaches would be necessary to clearly discriminate among these possible scenarios.

Biogeographic history

The biogeographic reconstruction based on the chloroplast dataset showed that sect. Armerina originated in the Himalayas and subsequently dispersed to the Hengduan Mountains, Northeastern QTP and Western QTP (Fig. 4). The lineages involved in these dispersal events further diversified in the Hengduan Mountains, Northeastern QTP and Western QTP, respectively, and gave rise to several of the extant species. Our dating analysis estimates that this section diverged from its closest relatives in the Pliocene about 3.55 Ma (95 % HPD: 1.76–5.93 Ma, Fig. 4). The timeframe of this event coincides with the recent uplift of the QTP, which occurred between 1.6 and 3.4 Ma [21, 22, 56]. A similar time of divergence was also observed in other groups of plants distributed in the QTP [32, 5759]. It has been suggested that the uplifts of the QTP might have limited the spread of many species, but accelerated speciation via vicariance [60]. The timeframe of the uplift also coincides with a period of high climatic oscillations that could have reinforced the processes initiated by the uplifts [57, 58].

Vicariance and dispersal triggered by the uplift of the QTP and associated climatic changes are common mechanisms in the diversification of plants in the QTP (e.g., [59, 61]), and also in other mountain areas (e.g., [1, 62, 63]). Based on the S-DIVA analysis, five of the 15 vicariance and seven of the 14 dispersal events account for cladogenetic events, and both events occurred during and after the Pliocene uplift of the QTP (Fig. 4). Vicariance and dispersal triggered by the uplift of the QTP and Quaternary climatic oscillations may accelerate the early diversification of sect. Armerina, and further shape the biogeographic patterns [59]. Furthermore, ten “vicariance” (for ease of notation, here we still keep the word “vicariance” for the isolation of populations of the same species) and seven dispersal events are identified within species-specific clades, which might play a role in promoting intraspecific divergence. Extensive inter- and intra-specific divergence took place in the QTP within the Pliocene and Quaternary climatic changes in many groups of plants (e.g., [6468]). Our analyses together with previous studies thus highlight the importance of the Pliocene uplift of the QTP and Quaternary climatic changes in promoting the diversification of plants in this mountain area.

Niche evolution of P. fasciculata

The S-DIVA analysis shows different biogeographic patterns for the two P. fasciculata clades. One clade (F2; Fig. 4) occupies only Northern Tibet, while samples from the other clade (F1) can be found in the Hengduan Mountains, Eastern Tibet and Northeastern QTP. Wiens & Donoghue [69] argued that phylogenetic niche conservatism and niche evolution might be critical in the biogeographic history of many groups. In contrast to most previous studies that have suggested the importance of niche conservatism in setting range limits and creating biogeographic patterns (e.g., [70, 71]), niche evolution under climatic changes seems to be the major factor explaining the biogeographic patterns detected here. Although the OU1 model that allows a single niche optimum is the best model along the temperature gradient (PC1), models that allow different niche optima received together higher AICc weights (AICcOUM+OUMV+OUMA = 0.59). This result suggests that ecological differentiation (i.e., different niche optima) is occurring in this group.

The two P. fasciculata clades and their two closely related species are estimated to diverge from each other during the Quaternary after the uplift of the QTP (Fig. 4). Climatic oscillations during the Quaternary had a dramatic effect on species distribution ranges [72]. Many species have repeatedly retreated and expanded their distributions following these climatic oscillations (e.g., [57, 58, 72, 73]. In the context of a changing environment, dispersal plays a crucial role in tracking favorable environmental conditions through space [74]. It can also help adaptation of small populations through both demographic and genetic rescue effects [75, 76]. Two dispersal events may have provided the opportunities for populations of clade F1 to occupy wide ranges and also invade new habitats and climatic regimes (Fig. 4). These events are associated with relatively relaxed niches (i.e., niche optima are not strongly correlated with temperature and precipitation gradient; Fig. 5) and fast niche evolution (Fig. 5) and these characteristics might have allowed these populations to adapt to the changing environmental conditions [11, 77]. In contrast, populations of clade F2 occur at higher altitudes (average 4600 m) compared to those of clade F1 (average 4200 m). These populations of clade F2 might have been adapted to a colder climate characterized by lower temperature seasonality (i.e., cooler summers) and less precipitation in the coldest Quarter (see Results of PCA and Fig. 5). Clade F2 displays a lower rate of niche evolution than F1 populations and this lower rate could have limited its dispersal into lower and warmer places. A similar pattern was observed in tropical treefrogs [78], which were unable to extend their ranges further North into temperate regions. Furthermore, recent climatic changes are involved in a shift toward higher elevations in the climatic envelopes of two closely related monkey-flower species in the direction of higher elevations [79]. However, given the harsh environmental conditions, it is plausible that climatic warming in the future might adversely affect the populations of the clade F2 and cause their distributions to shrink [57]. Our results also indicate that contrasting evolutionary processes can occur within closely related lineages, reinforcing the idea that phylogenetic niche conservatism is unlikely to hold at lower spatial scales [80].

While we focus on climatic variables (i.e., temperature and precipitation) to explain the biogeographic patterns detected here, additional ecological factors such as edaphic variables, competition, seed bank and seed number could be involved in creating biogeographic patterns [10, 12]. As argued by Hoskin et al. [81], geographic isolation of populations within species and variation in ecological factors are major precursors to cryptic speciation. The ecological differences and biogeographic patterns found between the two P. fasciculata clades may have given rise to some degree of differential adaptation to their respective environmental conditions, as also suggested in Taxus wallochiana [58]. However, our data is not appropriate to gain a detailed knowledge of the processes at play here and further studies involving a finer sampling of populations associated with large scale genomic data should be employed to better understand the mechanisms involved in the separation of the P. fasciculata clades.


Our phylogenetic analyses, based on both chloroplast and nuclear datasets, show non-monophyly of Primula sect. Armerina, corroborating the results of previous family-level studies [38, 39]. The topologies inferred from nuclear gene and concatenated chloroplast datasets are incongruent, which may mainly result from hybridization. This section was suggested to originate in the Himalayas during the Pliocene uplift of the QTP. Subsequent dispersals to the Hengduan Mountains, Northeastern QTP and Western QTP were considered as the consequence of the Pliocene uplift of the QTP and following climatic changes. We further provide a practicable framework for the first time to test the relationship between biogeographic patterns and ecological factors in the QTP area. Our evolutionary models suggest that niche evolution, rather than niche conservatism, seems to explain the biogeographic patterns of the two P. fasciculata clades.


Sampling and extraction

We collected in total 57 samples representing 10 of the 14 species belonging to Primula sect. Armerina (Additional file 5). We could not obtain plant material for P. iljinskyii, P. chrysostoma, P. knorringiana and P. valentinae, which have small distributions in Central Asian Mountains and are difficult to obtain due to their geographical locations. Widespread species were collected from different localities across their geographical ranges. For example, P. nutans was represented by two samples from N America, two from N Europe, one from NW Mongolia and four from China. Seven outgroup species were sampled based on the large phylogenetic tree of Primulaceae (Additional file 5; [38]). All samples were dried and stored in silica gel after collection, except for P. pamirica, which was obtained from Harvard University herbaria. The leaf tissues were ground to dust using an electric tissue homogenizer. Total genomic DNA was then isolated using the DNeasy Plant Mini Kit (Qiagen AG, Hombrechtikon, Switzerland) following the manufacturer’s instructions.

Amplification and sequencing

Five chloroplast DNA regions and one nuclear gene were sequenced. Three cpDNA loci (rpl16 intron; trnL-F region, which comprises the trnL intron and the trnL-trnF intergenic spacer; rps16 intron) were amplified and sequenced using the published primers [34]. The matK gene and trnH-psbA intergenic spacers were amplified and sequenced following the protocol described in Li et al. [82]. For the nuclear gene, we designed three pairs of exon-primed-intron-crossing (EPIC) primers based on an Arabidopsis thaliana translin family protein locus (tfp, AT2G03780) and a Primula sieboldii seedling cDNA library (FS228429). Only one pair of primers: tfp_e1.F (5’-CGAGAAAGGGTGGTAAAAGC-3’) and tfp_e1.R (5’-CTGGGGAGTAAGCTCGTCTG-3’), was amplified successfully for sect. Armerina. Polymerase chain reactions (PCR) generated double bands and direct sequencing of tfp_e1 amplicons produced electropherograms with double peaks and non-complementarity between sequenced strands in the following accessions: P. fasciculata (populations 9, 14, 16), P. conspersa (population 3) and P. egaliksensis. These PCR products were applied on a 1.5 % agarose gel, then excised and purified using a QIAquick Gel Extraction Kit (Qiagen, cat. no. 28704). The purified products were subsequently cloned into a pTZ57R/T vector and sequenced. Eight clones were sequenced per band.

All PCR reactions were performed in 25 μL volumes containing 1 × buffer (including 1.5 mM MgCl2), 2 mM MgCl2, 300 μM dNTPs, 0.2 μM of each primer and one unit Taq polymerase (GoTag DNA Polymerase, Promega, Madison, WI, USA). Amplifications were carried out on a thermocycler (Biometra, Goettingen, Germany) using the following conditions: a first cycle at 94 °C for 3 min; 36 cycles at 94 °C for 40 s, 55 °C for 1 min and 72 °C for 1.2 min; a final cycle of 7 min at 72 °C. All sequencing reactions used the Big Dye 3.1 Terminator cycle sequencing kit (Applied Biosystems, Foster City, CA, USA), then sequenced on an ABI Prism 3100 genetic analyzer (Applied Biosystems). DNA sequences were aligned with Geneious 6.1.6 (Biomatters) using MAFFT [83] and revised manually. The nuclear gene data generated from direct sequencing were scanned carefully and edited when necessary to ensure that all double peaks were identified correctly with standard degeneracy codes (e.g., Y means C or T; R means G or A; W means A or T; K means G or T; M means C or A). When double peaks were detected at a site, the site was ascertained as ambiguous only if the weakest signal reached at least 25 % of the peak signal strength [84, 85]. For individuals that contained multiple clones for the tfp gene, we randomly chose a single representative sequence for the phylogenetic analysis if all the clones formed a well-supported clade in a preliminary analysis, while multiple sequences were retained otherwise. All sequences were submitted to GenBank (accessions KT259477-KT259852).

Phylogenetic reconstruction and molecular dating

The five chloroplastic genes were concatenated into a single dataset using SequenceMatrix 1.7.8 [86], while the chloroplast and nuclear datasets were analyzed separately. The GTR + G model of sequence evolution was selected on the basis of the Akaike information criterion (AIC) for all DNA regions as estimated by jModelTest 2.1.4 [87]. Maximum likelihood analyses were done with PhyML (ver. 3.0; [88]) using the BEST algorithm for branch swapping and 103 bootstrap replicates to assess node support. We estimated tree topology by Bayesian inference using MrBayes 3.2 [89] with the GTR + G model of evolution and default priors. We unlinked the parameters of the GTR + G model between the five different genes for the analysis of the chloroplast dataset. We repeated the MrBayes analyses three times for each analysis (i.e., chloroplast and nuclear dataset) and each analysis consisted of four chains of 107 generations, sampling every 103 steps with temperature parameter set to 0.1. We determined convergence by examining trace plots of the log-likelihood values for each parameter in Tracer 1.5.

We used the chloroplast dataset for dating analysis with a secondary calibration strategy, as described in de Vos et al. [38], However, age estimation obtained from this kind of calibration may be inherently subjected to bias and errors [90]. We addressed this concern by comparing our estimated age with previously published ones, but it should be noted that they are estimates that should be treated with caution. Divergence time analysis was performed in BEAST (ver. 1.7; [91]). The fossil record of Primulaceae is too sparse to provide multiple and reliable calibrations within the family [92, 93]. The only available fossil that can be used as minimum-age estimate for the split between Primula and Soldanella is represented by seeds from Primula riosiae from the Miocene that are dated at 15.97 Ma (the early-mid Miocene boundary; [94]). Therefore, we performed a completely separate divergence-time analysis from a taxonomically more inclusive sample of six plastid gene regions (matK, ndhF, rbcL, trnL-F, rps16 and rpl16) available in GenBank (Additional file 6). We included P. fasciculata, P. involucrata, P. sikkimensis and P. alpicola in the larger analyses to obtain a root age estimate for Primula sect. Armerina. The resulting data matrix comprised 7978 aligned sites and 13 species of Primulaceae, with 8.3 % missing data (Additional file 6). Sequence alignment and model specification proceeded as described above, unless otherwise stated. The GTR + G model of sequence evolution was selected by jModelTest 2.1.4 for rpl16, trnL-F, ndhF and matK, GTR + I model for rps16 and HKY + G + I model for rbcL. A normally distributed prior with a mean of 39.996 Ma and a standard deviation of 11.492 Ma [38] was used to constrain the root of the Soldanella/Androsace divergence to be within the interval 21.09–58.90 Ma with 95 % probability. The calibration point between Primula and Soldanella based on the fossil of Primula riosiae was set to a lognormal prior with an offset of 15.97, a mean of 2.1 and a standard deviation of 0.63. The analyses were run using a random starting tree for 108 generations sampling every 103 generations under the uncorrelated lognormal relaxed clock model, a birth-death tree prior and the selected models of substitution for different partitions. The analyses were repeated three times to verify convergence by examining the posterior distribution of parameters in Tracer 1.5. After the removal of the burn-in (10 million generations in each analysis, corresponding to 10 % of the samples), the inferred age distribution of the node separating the groups containing either P. fasciculata and P. involucrata or P. sikkimensis and P. alpicola was estimated in Tracer 1.5.

The age obtained for the Armerina section was then used as a calibration point for the root age of the Armerina analysis and modeled as a γ prior with a shape of 9.7, a scale of 0.61 and an offset of 1.4. We used similar settings as described above and the samples retained after removal of the burn-in from the three runs were summarized as a maximum clade credibility tree with mean divergence times using TreeAnnotator (part of the BEAST package).

Biogeographic reconstruction

We ran Statistical Dispersal Vicariance Analysis (S-DIVA) using RASP v.2.1 [95, 96] to infer the biogeographic history of this section based on the phylogenetic trees constructed only from our concatenated chloroplast dataset. We did not use the tfp nuclear dataset since two homologous copies were obtained from some samples, but multiple copies were not present in all species. We defined seven biogeographic regions for the individuals that were collected: A (East Tibet and Hengduan Mountains), B (Himalayas Mountains), C (Northeast QTP), D (Monutains of Central Asia), E (North Europe), F (North America) and G (Mongolian Plateau). Regions A-C were defined according to the biogeographic divisions of China [97], and had been applied in other studies (e.g., [59, 98]). Region D was defined based on the distribution area of P. pamirica. Regions E-G were defined based on the distribution of P. nutans and some outgroup samples used in this study. To account for uncertainties in phylogenetic reconstructions, we randomly chose 20,000 trees from the posterior distribution of trees obtained by BEAST. The number of maximum areas was set to 2 and we estimated the possible ancestral ranges at each node of the selected phylogenetic trees.

Evolution of ecological preferences

Climatic niche is one of the main factors for setting historically biogeographic patterns, especially during drastically climatic changes, such as Quaternary climate oscillations [1012]. In order to better understand the biogeographic patterns obtained above, we fitted a series of macro-evolutionary models based on 19 bioclimatic variables. We focused on the clades formed by the species P. fasciculata, P. tibetica and P. nutans because they represent the main lineages in the group, and tested whether the evolutionary trajectories of the climatic niches differed among the different clades (F1, F2; Fig. 4) obtained for P. fasciculata (see Results) and its two closely related species P. nutans and P. tibetica (NT; Fig. 4). For this test, we used only the samples of P. nutans that were collected in the QTP.

We extracted the 19 bioclimatic variables of WorldClim (; [99]) for all samples of the three groups (F1, F2, NT) using the package raster [100] in R. All the 19 bioclimatic variables were then summarized into principle components using the prcomp function in the stats package of R [101]. We used the R package OUwie [102] to compare the fit of a series of models (see Table 1 for detailed interpretation for each model) to explain the differences in niche evolution between species inhabiting similar or different biogeographic regions. We tested these models on different sets of groups: (1) F1/F2 vs. NT (SET1); (2) F1 vs. F2 vs. NT (SET2); (3) F1 vs. F2/NT (SET3); and (4) F2 vs. F1/NT (SET4).

Table 1 Models of niche evolution relevant to different group-sets with their parameters and interpretation, indicating for each model whether the optimal niche value, θ, the intensity of random fluctuations in the evolutionary trajectory, σ2, and the strength of selection toward the optimal value, α, are modeled with one global parameter or with two or three parameters that are group-specific

Stochastic mapping for all model tests were run 10 times for 100 trees randomly selected from the posterior distribution of trees from the BEAST analysis to account for possible uncertainty in the estimated values. Model fit was determined using AICc weights calculated from ΔAICc scores [103]. The highest value of AICc weight represents the best model. Finally, we calculated an average AICc weight and lower (2.5 %) and upper (97.5 %) quantiles of the distributions of AICc weights for each evolutionary niche model.



Qinghai-Tibet Plateau


Maximum Likelihood


Posterior Probability


Principal Component Analysis




Polymerase chain reactions


Akaike information criterion


Statistical Dispersal Vicariance Analysis


Maximum Clade Credibility


  1. 1.

    Hoorn C, Wesselingh FP, Ter Steege H, Bermudez MA, Sevink J, et al. Amazonia through time: Andean uplift, climate change, landscape evolution, and biodiversity. Science. 2010;330(6006):927–31.

    CAS  PubMed  Google Scholar 

  2. 2.

    Favre A, Päckert M, Pauls SU, Jähnig SC, Uhl D, Michalak I, et al. The role of the uplift of the Qinghai-Tibetan plateau for the evolution of Tibetan biotas. Biol Rev. 2014;90(1):236–53.

    PubMed  Google Scholar 

  3. 3.

    Wen J, Zhang JQ, Nie ZL, Zhong Y, Sun H. Evolutionary diversifications of plants on the Qinghai-Tibetan Plateau. Front Genet. 2014;5:4.

    PubMed Central  PubMed  Google Scholar 

  4. 4.

    Mayr E. Animal Species and Evolution. Cambridge, Mass: Harvard University Press; 1963.

    Google Scholar 

  5. 5.

    Coyne JA. Genetics and speciation. Nature. 1992;355:511–5.

    CAS  PubMed  Google Scholar 

  6. 6.

    Rice WR, Hostert EE. Laboratory experiments on speciation: what have we learned in 40 years. Evolution. 1993;47:1637–53.

    Google Scholar 

  7. 7.

    Richardson JE, Pennington RT, Pennington TD, Hollingsworth PM. Rapid diversification of a species-rich genus of neotropical rain forest trees. Science. 2001;293:2242–5.

    CAS  PubMed  Google Scholar 

  8. 8.

    Huang WP, Ji HX. Discovery of Hipparion fauna in Xizang. Chin Sci Bull. 1979;19:006.

    Google Scholar 

  9. 9.

    Spicer RA, Harris NB, Widdowson WM, Herman AB, Guo S, Valdes PJ, et al. Constant elevation of southern Tibet over the past 15 million years. Nature. 2003;421:622–4.

    CAS  PubMed  Google Scholar 

  10. 10.

    Wiens JJ. The niche, biogeography and species interactions. Philos Trans R Soc B Biol Sci. 2011;366:2336–50.

    Google Scholar 

  11. 11.

    Sexton JP, McIntyre PJ, Angert AL, Rice KJ. Evolution and ecology of species range limits. Ann Rev Ecol Evol Syst. 2009;40:415–36.

    Google Scholar 

  12. 12.

    Lomolino MV, Riddle BR, Brown JH. Biogeography. 3rd ed. Sunderland, MA: Sinauer Associates; 2006.

    Google Scholar 

  13. 13.

    Peterson AT, Soberón J, Sánchez-Cordero V. Conservatism of ecological niches in evolutionary time. Science. 1999;285(5431):1265–7.

    CAS  PubMed  Google Scholar 

  14. 14.

    Wiens JJ, Graham CM. Niche conservatism: integrating evolution, ecology, and conservation biology. Ann Rev Ecol Evol Syst. 2005;36:519–39.

    Google Scholar 

  15. 15.

    Pfenninger M, Nowak C, Magnin F. Intraspecific range dynamics and niche evolution in Candidula land snail species. Biol J Linn Soc. 2007;90(2):303–17.

    Google Scholar 

  16. 16.

    Evans MEK, Smith SA, Flynn RS, Donoghue MJ. Climate, Niche Evolution, and Diversification of the “Bird-Cage” Evening Primroses (Oenothera, Sections Anogra and Kleinia). Am Nat. 2009;173(2):225–40.

    PubMed  Google Scholar 

  17. 17.

    Kostikova A, Litsios G, Salamin N, Pearman PB. Linking life-history traits, ecology, and niche breadth evolution in North American eriogonoids (Polygonaceae). Am Nat. 2013;182(6):760–74.

    PubMed  Google Scholar 

  18. 18.

    Wu CY. Flora of Tibet. Beijing: Science Press; 1987.

    Google Scholar 

  19. 19.

    Royden LH, Burchfiel BC, van der Hilst RD. The geological evolution of the Tibetan plateau. Science. 2008;321:1054–8.

    CAS  PubMed  Google Scholar 

  20. 20.

    Harrison TM, Copland P, Kidd WSF, Yin A. Raising Tibet. Science. 1992;255:1663–70.

    CAS  PubMed  Google Scholar 

  21. 21.

    Li JJ, Shi YF, Li BY. Uplift of the Qinghai-Xizang (Tibet) Plateau and Global Change. Lanzhou: Lanzhou University Press; 1995.

    Google Scholar 

  22. 22.

    Shi YF, Tang MC, Ma YZ. The relation of second rising in Qinghai-Tibetan Plateau and Asian monsoon. Sci China Ser D. 1998;28:263–71.

    Google Scholar 

  23. 23.

    Guo ZT, Ruddiman WF, Hao QZ, Wu HB, Qiao YS, Zhu RX, et al. Onset of Asian desertification by 22 Myr ago inferred from loess deposits in China. Nature. 2002;416:159–63.

    CAS  PubMed  Google Scholar 

  24. 24.

    Wu ZY, Wu SG. A proposal for a new floristic kingdom (realm) – the E Asiatic kingdom, its delimitation and characteristics. In: Zhang A, Wu S, editors. Floristic Characteristics and Diversity of East Asian Plants. Beijing: China Higher Education Press; 1996. p. 3–42.

    Google Scholar 

  25. 25.

    Liu JQ, Duan YW, Hao G, Ge XJ, Sun H. Evolutionary history and underlying adaptation of alpine plants on the Qinghai-Tibet Plateau. J Syst Evol. 2014;52(3):241–9.

    Google Scholar 

  26. 26.

    Richards AJ. Primula. Portland, Oregon: Timber Press, USA; 2003.

    Google Scholar 

  27. 27.

    Hu CM, Kelso S. Primulaceae. In: Wu ZY, Raven PH, editors. Flora of China, Vol 15. Beijing: Science Press; St Louis: Missouri Botanical Garden Press; 1996. p. 99–185.

    Google Scholar 

  28. 28.

    Mao KS, Hao G, Liu JQ, Adams RP, Milne RI. (2010). Diversification and biogeography of Juniperus (Cupressaceae): variable diversification rates and multiple intercontinental dispersals. New Phytol. 2010;188(1):254–72.

    CAS  PubMed  Google Scholar 

  29. 29.

    Mao KS, Milne RI, Zhang LB, Peng YL, Liu JQ, Thomas P, et al. Distribution of living Cupressaceae reflects the breakup of Pangea. Proc Natl Acad Sci USA. 2012;109:7793–8.

    PubMed Central  CAS  PubMed  Google Scholar 

  30. 30.

    Sun YS, Wang AL, Wan DS, Wang Q, Liu JQ. Rapid radiation of Rheum (Polygonaceae) and parallel evolution of morphological traits. Mol Phylogenet Evol. 2012;63:150–8.

    PubMed  Google Scholar 

  31. 31.

    Zhang XL, Wang YJ, Ge XJ, Yuan YM, Yang HL, Liu JQ. Molecular phylogeny and biogeography of Gentiana sect. Cruciata (Gentianaceae) based on four chloroplast DNA datasets. Taxon. 2009;58(3):862–70.

    Google Scholar 

  32. 32.

    Liu JQ, Gao TG, Chen ZD, Lu AM. Molecular phylogeny and biogeography of the Qinghai–Tibetan Plateau endemic Nannoglottis (Asteraceae). Mol Phylogenet Evol. 2002;23:307–25.

    CAS  PubMed  Google Scholar 

  33. 33.

    Kelso S. Taxonomy of Primula sects. Aleuritia and Armerina in North America. Rhodora. 1991;93:67–99.

    Google Scholar 

  34. 34.

    Kelso S. The genus Primula as a model for evolution in the Alaskan flora. Arct Antarct Alp Res. 1992;24:82–7.

    Google Scholar 

  35. 35.

    Guggisberg A, Mansion G, Kelso S, Conti E. Evolution of biogeographic patterns, ploidy levels, and breeding systems in a diploid–polyploid species complex of Primula. New Phytol. 2006;171(3):617–32.

    CAS  PubMed  Google Scholar 

  36. 36.

    Guggisberg A, Baroux C, Grossniklaus U, Conti E. Genomic origin and organisation of the allopolyploid Primula egaliksensis investigated by in situ hybridization. Ann Bot. 2008;101:919–27.

    PubMed Central  CAS  PubMed  Google Scholar 

  37. 37.

    Guggisberg A, Mansion G, Conti E. Disentangling reticulate evolution in an arctic–alpine polyploid complex. Syst Biol. 2009;58:55–73.

    CAS  PubMed  Google Scholar 

  38. 38.

    de Vos JM, Hughes CE, Schneeweiss GM, Moore BR, Conti E. Heterostyly accelerates diversification via reduced extinction in primroses. P Roy Soc B-Biol Sci. 2014;281(1784):20140075.

    Google Scholar 

  39. 39.

    Mast AR, Kelso S, Conti E. Are any primroses (Primula) primitively monomorphic? New Phytol. 2006;171(3):605–16.

    CAS  PubMed  Google Scholar 

  40. 40.

    Zhang YX, Zeng CX, Li DZ. Complex evolution in Arundinarieae (Poaceae: Bambusoideae): Incongruence between plastid and nuclear GBSSI gene phylogenies. Mol Phylogenet Evol. 2012;63(3):777–97.

    PubMed  Google Scholar 

  41. 41.

    Alström P, Höhna S, Gelang M, Eriscon PG, Olsson U. Non-monophyly and intricate morphological evolution within the avian family Cettiidae revealed by multilocus analysis of a taxonomically densely sampled dataset. BMC Evol Biol. 2011;11(1):352.

    PubMed Central  PubMed  Google Scholar 

  42. 42.

    Yu WB, Huang PH, Li DZ, Wang H. Incongruence between Nuclear and Chloroplast DNA Phylogenies in Pedicularis Section Cyathophora (Orobanchaceae). PLoS ONE. 2013;8(9), e74828.

    PubMed Central  CAS  PubMed  Google Scholar 

  43. 43.

    Degnan JH, Rosenberg NA. Gene tree discordance, phylogenetic inference and the multispecies coalescent. Trends Ecol Evol. 2009;24:332–40.

    PubMed  Google Scholar 

  44. 44.

    Casazza G, Granato L, Minuto L, Conti E. Polyploid evolution and Pleistocene glacial cycles: a case study from the Alpine species Primula marginata. BMC Evol Biol. 2012;12(1):56.

    PubMed Central  PubMed  Google Scholar 

  45. 45.

    Schmidt-Lebuhn AN, de Vos JM, Keller B, Conti E. Phylogenetic analysis of Primula section Primula reveals rampant non-monophyly among morphologically distinct species. Mol Phylogenet Evol. 2012;65:23–34.

    PubMed  Google Scholar 

  46. 46.

    Petit RJ, Excoffier L. Gene flow and species delimitation. Trends Ecol Evol. 2009;24:386–93.

    PubMed  Google Scholar 

  47. 47.

    Hansson B, Tarka M, Dawson DA, Horsburgh GJ. Hybridization but no evidence for backcrossing and introgression in a sympatric population of Great Reed Warblers and Clamorous Reed Warblers. PLoS ONE. 2012;7(2), e31667.

    PubMed Central  CAS  PubMed  Google Scholar 

  48. 48.

    Rieseberg LH. Hybrid origins of plant species. Annu Rev Ecol Syst. 1997;28:359–89.

    Google Scholar 

  49. 49.

    Woodell SRJ. Natural hybridization in Britain between Primula vulgaris Huds (the primrose) and P. elatior (L.) Hill (the oxlip). Watsonia. 1969;7:115–27.

    Google Scholar 

  50. 50.

    Zhu XF, Li Y, Wu GL, Fang ZD, Li QJ, Liu JQ. Molecular and morphological evidence for natural hybridization between Primula secundiflora Franchet and P. poissonii Franchet (Primulaceae). Acta Biol Cracov Ser Bot. 2009;51:29–36.

    Google Scholar 

  51. 51.

    Ma YP, Tian XL, Zhang JL, Wu ZK, Sun WB. Evidence for natural hybridization between Primula beesiana and P. bulleyana, two heterostylous primroses in NW Yunnan, China. J Syst Evol. 2014;52(4):500–7.

    Google Scholar 

  52. 52.

    Buckley TR, Cordeiro M, Marshall DC, Simon C. Differentiating between hypotheses of lineage sorting and introgression in New Zealand alpine cicadas (Maoricicada Dugdale). Syst Biol. 2006;55(3):411–25.

    PubMed  Google Scholar 

  53. 53.

    Doyle JJ. Gene trees and species trees: molecular systematics as one-character taxonomy. Syst Bot. 1992;17(1):144–63.

    Google Scholar 

  54. 54.

    Joly S, McLenachan PA, Lockhart PJ. A statistical approach for distinguishing hybridization and incomplete lineage sorting. Am Nat. 2009;174(2):54–70.

    Google Scholar 

  55. 55.

    Liu L, Pearl DK. Species trees from gene trees: reconstructing Bayesian posterior distributions of a species phylogeny using estimated gene tree distributions. Syst Biol. 2007;56(3):504–14.

    CAS  PubMed  Google Scholar 

  56. 56.

    Shi YF, Li JJ, Li BY. Uplift of the Qinghai–Tibetan Plateau and East Asia environmental changes during the late Cenozoic. Acta Geogra Sin. 1999;54:10–21.

    Google Scholar 

  57. 57.

    Li L, Abbott RJ, Liu BB, Sun YS, Li LL, Zou JB, et al. Pliocene intraspecific divergence and Plio-Pleistocene range expansions within Picea likiangensis (Lijiang spruce), a dominant forest tree of the Qinghai-Tibet Plateau. Mol Ecol. 2013;22:5237–55.

    CAS  PubMed  Google Scholar 

  58. 58.

    Liu J, Moller M, Provan J, Gao LM, Poudel RM, Li DZ. Geological and ecological factors drive cryptic speciation of yews in a biodiversity hotspot. New Phytol. 2013;199:1093–108.

    PubMed  Google Scholar 

  59. 59.

    Zhou Z, Hong D, Niu Y, Li G, Nie Z, Wen J, et al. Phylogenetic and biogeographic analyses of the Sino-Himalayan endemic genus Cyananthus (Campanulaceae) and implications for the evolution of its sexual system. Mol Phylogenet Evol. 2013;68(3):482–97.

    PubMed  Google Scholar 

  60. 60.

    Yang S, Dong H, Lei F. Phylogeography of regional fauna on the Tibetan Plateau: a review. Prog Nat Sci. 2009;19(7):789–99.

    CAS  Google Scholar 

  61. 61.

    Zhang JW, Nie ZL, Wen J, Sun H. Molecular phylogeny and biogeography of three closely related genera, Soroseris, Stebbinsia, and Syncalathium (Asteraceae, Cichorieae), endemic to the Tibetan Plateau. SW China Taxon. 2011;60:15–26.

    Google Scholar 

  62. 62.

    Schönswetter P, Tribsch A. Vicariance and dispersal in the alpine perennial Bupleurum stellatum L. (Apiaceae). Taxon. 2005;54:725–32.

    Google Scholar 

  63. 63.

    Bryja J, Mikula O, Patzenhauerová H, Oguge NO, Šumbera R, Verheyen E. The role of dispersal and vicariance in the Pleistocene history of an East African mountain rodent, Praomys delectorum. J Biogeogr. 2014;41(1):196–208.

    Google Scholar 

  64. 64.

    Xu TT, Abbott RJ, Milne RI, Mao KS, Du FK, Wu GL, et al. Phylogeography and allopatric divergence of cypress species (Cupressus L.) in the Qinghai-Tibetan Plateau and adjacent regions. BMC Evol Biol. 2010;10(1):194.

    PubMed Central  PubMed  Google Scholar 

  65. 65.

    Wang YJ, Liu JQ, Miehe G. Phylogenetic origins of the Himalayan endemic Dolomiaea, Diplazoptilon and Xanthopappus (Asteraceae: Cardueae) based on three DNA regions. Ann Bot. 2007;99(2):311–22.

    PubMed Central  CAS  PubMed  Google Scholar 

  66. 66.

    Li Y, Stocks M, Hemmila S, Källman T, Zhu H, Zhou Y, et al. Demographic histories of four spruce (Picea) species of the Qinghai-Tibetan Plateau and neighboring areas inferred from multiple nuclear loci. Mol Biol Evol. 2010;27:1001–14.

    PubMed  Google Scholar 

  67. 67.

    Liu JQ, Wang YJ, Wang AL, Hideaki O, Abbott RJ. Radiation and diversification within the Ligularia–Cremanthodium–Parasenecio complex (Asteraceae) triggered by uplift of the Qinghai-Tibetan Plateau. Mol Phylogenet Evol. 2006;38(1):31–49.

    CAS  PubMed  Google Scholar 

  68. 68.

    Jia DR, Abbott RJ, Liu TL, Mao KS, Bartish IV, Liu JQ. Out of the Qinghai–Tibet Plateau: evidence for the origin and dispersal of Eurasian temperate plants from a phylogeographic study of Hippophaë rhamnoides (Elaeagnaceae). New Phytol. 2012;194(4):1123–33.

    PubMed  Google Scholar 

  69. 69.

    Wiens JJ, Donoghue MJ. Historical biogeography, ecology and species richness. Trends Ecol Evol. 2004;19(12):639–44.

    PubMed  Google Scholar 

  70. 70.

    Stephens PR, Wiens JJ. Bridging the gap between community ecology and historical biogeography: niche conservatism and community structure in emydid turtles. Mol Ecol. 2009;18(22):4664–79.

    CAS  PubMed  Google Scholar 

  71. 71.

    Kozak KH, Wiens JJ. Niche conservatism drives elevational diversity patterns in Appalachian salamanders. Am Nat. 2010;176(1):40–54.

    PubMed  Google Scholar 

  72. 72.

    Hewitt GM. Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc B Biol Sci. 2004;359(1442):183–95.

    CAS  Google Scholar 

  73. 73.

    Ren GP, Abbott RJ, Zhou YF, Zhang LR, Peng YL, Liu JQ. Genetic divergence, range expansion and possible homoploid hybrid speciation among pine species in Northeast China. Heredity. 2012;108(5):552–62.

    PubMed Central  CAS  PubMed  Google Scholar 

  74. 74.

    Polechova J, Barton N, Marion G. Species’ range: adaptation in space and time. Am Nat. 2009;174:186–204.

    Google Scholar 

  75. 75.

    Alleaume‐Benharira M, Pen IR, Ronce O. Geographical patterns of adaptation within a species’ range: interactions between drift and gene flow. J Evol Biol. 2006;19(1):203–15.

    PubMed  Google Scholar 

  76. 76.

    Kawecki TJ. Adaptation to marginal habitats. Ann Rev Ecol Evol Syst. 2008;39:321–42.

    Google Scholar 

  77. 77.

    Lavergne S, Mouquet N, Thuiller W, Ronce O. Biodiversity and climate change: integrating evolutionary and ecological responses of species and communities. Ann Rev Ecol Evol Syst. 2010;41:321–50.

    Google Scholar 

  78. 78.

    Wiens JJ, Graham CH, Moen DS, Smith SA, Reeder TW. Evolutionary and ecological causes of the latitudinal diversity gradient in hylid frogs: treefrog trees unearth the roots of high tropical diversity. Am Nat. 2006;168(5):579–96.

    PubMed  Google Scholar 

  79. 79.

    Angert AL. The niche, limits to species’ distributions, and spatiotemporal variation in demography across the elevation ranges of two monkeyflowers. Proc Natl Acad Sci USA. 2009;106:19693–8.

    PubMed Central  CAS  PubMed  Google Scholar 

  80. 80.

    Kostikova A, Litsios G, Burgy S, Milani L, Pearman PB, Salamin N. Scale‐dependent adaptive evolution and morphological convergence to climatic niche in Californian eriogonoids (Polygonaceae). J Biogeogr. 2014;41(7):1326–37.

    Google Scholar 

  81. 81.

    Hoskin CJ, Higgie M, McDonald KR, Moritz C. Reinforcement drives rapid allopatric speciation. Nature. 2005;437:1353–6.

    CAS  PubMed  Google Scholar 

  82. 82.

    Li DZ, Gao LM, Li HT, Wang H, Ge XJ, Liu JQ, et al. Comparative analysis of a large dataset indicates that internal transcribed spacer (ITS) should be incorporated into the core barcode for seed plants. Proc Natl Acad Sci USA. 2011;108:19641–6.

    PubMed Central  CAS  PubMed  Google Scholar 

  83. 83.

    Katoh K, Toh H. Recent developments in the MAFFT multiple sequence alignment program. Brief Bioinform. 2008;9:286–98.

    CAS  PubMed  Google Scholar 

  84. 84.

    Fuertes AJ, Rosselló JA, Nieto FG. Molecular evidence for the compilospecies model of reticulate evolution in Armeria (Plumbaginaceae). Syst Biol. 1999;48:735–54.

    Google Scholar 

  85. 85.

    Fuertes AJ, Nieto FG. Additive polymorphisms and reticulation in an ITS phylogeny of thrifts (Armeria, Plumbaginaceae). Mol Phylogenet Evol. 2003;28:430–47.

    Google Scholar 

  86. 86.

    Vaidya G, Lohman DJ, Meier R. SequenceMatrix: concatenation software for the fast assembly of multi-gene datasets with character set and codon information. Cladistics. 2011;27:171–80.

    Google Scholar 

  87. 87.

    Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.

    CAS  PubMed  Google Scholar 

  88. 88.

    Guindon S, Delsuc F, Dufayard JF, Gascuel O. Estimating maximum likelihood phylogenies with PhyML. Methods Mol Biol. 2009;537:113–37.

    CAS  PubMed  Google Scholar 

  89. 89.

    Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19:1572–74.

    CAS  PubMed  Google Scholar 

  90. 90.

    Sauquet H, Ho SYW, Gandolfo MA, Jordan GJ, Wilf P, Cantrill DJ, et al. Testing the impact of calibration on molecular divergence times using a fossil-rich group: the case of Nothofagus (Fagales). Syst Biol. 2012;61(2):289–313.

    PubMed  Google Scholar 

  91. 91.

    Drummond A, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.

    PubMed Central  PubMed  Google Scholar 

  92. 92.

    Schneeweiss GM, Schönswetter P, Kelso S, Niklfeld H. Complex biogeographic patterns in Androsace (Primulaceae) anstrd related genera: evidence from phylogenetic analyses of nuclear internal transcribed spacer and plastid trnL-F sequences. Syst Biol. 2004;53:856–76.

    PubMed  Google Scholar 

  93. 93.

    Boucher FC, Thuiller W, Roquet C, Douzet R, Aubert S, Alvarez N, et al. Reconstructing the origins of high-alpine niches and cushion life form in the genus Androsace s.l. (Primulaceae). Evolution. 2012;66:1255–68.

    PubMed Central  PubMed  Google Scholar 

  94. 94.

    Czaja A. Paläokarpologische Untersuchungen von Taphozönosen des Unter- und Mittelmiozäns aus dem Braunkohlentagebau Berzdorf/Oberlausitz (Sachsen). Palaeontogr Abt B. 2003;265:1–148.

    Google Scholar 

  95. 95.

    Yu Y, Harris AJ, He XJ. S-DIVA (Statistical Dispersal-Vicariance Analysis): a tool for inferring biogeographic histories. Mol Phylogenet Evol. 2010;56:848–50.

    PubMed  Google Scholar 

  96. 96.

    Yu Y, Harris AJ, He XJ. RASP (Reconstruct Ancestral State in Phylogenies): a tool for historical biogeography. Mol Phylogenet Evol. 2015;87:46–9.

    PubMed  Google Scholar 

  97. 97.

    Xie Y, MacKinnon J, Li D. Study on biogeographical divisions of China. Biodivers Conserv. 2004;13(7):1391–417.

    Google Scholar 

  98. 98.

    Roquet C, Boucher FC, Thuiller W, Lavergne S. Replicated radiations of the alpine genus Androsace (Primulaceae) driven by range expansion and convergent key innovations. J Biogeogr. 2013;40(10):1874–86.

    PubMed Central  PubMed  Google Scholar 

  99. 99.

    Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25(15):1965–78.

    Google Scholar 

  100. 100.

    Hijmans R, van Etten J. raster: geographic analysis and modeling with raster data. 2012.

  101. 101.

    Venables WN, Ripley BD. Modern Applied Statistics with S. 4th ed. New York: Springer; 2002.

    Google Scholar 

  102. 102.

    Beaulieu JM, Jhwueng DC, Boettiger C, O’Meara BC. Modeling stabilizing selection: expanding the Ornstein-Uhlenbeck model of adaptive evolution. Evolution. 2012;66:2369–83.

    PubMed  Google Scholar 

  103. 103.

    Burnham KP, Anderson DR. Model selection and multimodel inference:a practical information theoretic approach. New York: Springer; 2002.

Download references


We would like to thank Prof. JQ Liu for supporting the fieldwork, C. Yang for help in the field, Dr. B. Tian for providing sequences of Primula zambalensis, the Harvard University herbaria for providing samples of P. pamirica and D. Silvestro, G. Litsios, A. Kostikova & M. Serrano for their great help in analyzing data. We also profited from insightful discussions with M. Nowak. We are grateful for his help in the lab work and the nuclear primers he provided. The paper benefited greatly from the comments received from two anonymous reviewers. This work was funded by the University of Lausanne research fund, the grant 31003A_138282 from the Swiss National Science Foundation to NS and the China Scholarship Council (award to GPR for four year’s PhD study abroad at the University of Lausanne). We received support for computational work from the Vital-IT facilities from the Swiss Institute of Bioinformatics.

Author information



Corresponding author

Correspondence to Nicolas Salamin.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

NS, GPR and EC conceived the study. GPR carried out the sampling and the lab work, performed the molecular analysis and drafted the manuscript. NS supervised GPR, coordinated the project and helped to draft the manuscript. EC participated in the design and revising of the manuscript. All authors read, reviewed and approved the final manuscript.

Additional files

Additional file 1:

The 19 bioclimatic variables for each samples of the three groups used for the niche models. (DOCX 153 kb)

Additional file 2:

The PC1 and PC2 values summarized from the 19 bioclimatic variables used for the niche models. (DOCX 79 kb)

Additional file 3:

Average AICc weights for each niche evolutionary model. The values are averages across estimations done on 100 trees with 10 stochastic maps. Quantiles (reported in brackets below each average AICc weight) are calculated as 2.5 % and 97.5 % from the distribution of AICc weights based on 100 trees, with 10 stochastic maps. The bold values show the best-fit models for each set and each PC axis. (DOCX 86 kb)

Additional file 4:

Model fit and estimates parameters of supported OUwie models for the four group-sets. Parameter estimates are averages across estimations done on 100 trees with 10 stochastic maps. Quantiles (reported in brackets) are calculated as 2.5 % and 97.5 % from the distribution of AICc weights based on 100 trees, with 10 stochastic maps. F1, F2 and NT are three groups defined based on the chloroplast tree (see Fig. 4). (DOCX 119 kb)

Additional file 5:

List of samples used in the present study. (DOCX 125 kb)

Additional file 6:

GenBank accession numbers for DNA sequence data of the 13 taxa in the family Primulaceae chloroplast DNA dataset that was used to provide a secondary calibration for the section Armerina dataset. (DOCX 81 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ren, G., Conti, E. & Salamin, N. Phylogeny and biogeography of Primula sect. Armerina: implications for plant evolution under climate change and the uplift of the Qinghai-Tibet Plateau. BMC Evol Biol 15, 161 (2015).

Download citation


  • Incomplete Lineage Sorting
  • Biogeographic Pattern
  • Hengduan Mountain
  • Outgroup Species
  • Bioclimatic Variable