Westwards and northwards dispersal of Triosteum himalayanum (Caprifoliaceae) from the Hengduan Mountains region based on chloroplast DNA phylogeography

Material sampling

Field expeditions were conducted between 2006 and 2015 to cover all of the possible distribution localities of T. himalayanum according to specimen records in the herbaria. However, populations of T. himalayanum that were previously distributed in central China (according to the specimen study) were difficult to collect due to land use changes. Fresh leaves of 4–25 individuals were collected for each population, separated at least 20 m apart. A total of 238 individuals from 20 populations were collected. Among these populations, 19 were from the HHM, and only one population (P20) was not from this region (Table 1; Fig. 1). Leaf materials were dried in silica gel and stored at room temperature. Voucher specimens of all populations are deposited in the Herbarium of Northwest Institute of Plateau Biology (HNWP), Xining, Qinghai, China.

Lonicera trichosantha Bureau & Franchet (Caprifoliaceae), has been chosen as out group for haplotype phylogenetic analyses.

DNA extraction, polymerase chain reaction amplification, and sequencing

Total DNA was extracted from the silica-gel dried leaf materials following a modified CTAB (Cetyltrimethyl Ammonium Bromide; Doyle & Doyle, 1987). Samples were genetically screened using the universal cpDNA primers described by Hamilton (1999) and Dong et al. (2012). After preliminary screening (pilot screening) of more than 10 intergenic spacers, rbcL-accD, rps15-ycf1, and trnH-psbA were chosen for further analysis, since these three intergenic spacers revealed sequence polymorphism. Polymerase chain reactions (PCRs) were performed in 25 μL volume containing 0.6 μL of template DNA (approximately 20 ng), 2.5 μL of 10 × PCR buffer (with Mg²⁺), 1.0 μL of 10 mM dNTPs, 0.5 μL of 5 pM of each primer, and 0.2 μL (1.5 units) of Taq polymerase (Takara, Dalian, China). Amplification conditions were: 5 min at 95 °C; followed by 35 cycles of 50 s at 95 °C, 1 min at 58 °C, and 1 min at 72 °C; with a final extension of 6 min at 72 °C. PCR products were purified using a CASpure PCR Purification Kit (CASarry, Shanghai, China) following the manufacturer’s protocols. Sequencing reactions and analysis were performed using an ABI 3730xl DNA sequencer (Applied Biosystems, Foster City, CA, USA). Primers used for sequencing were the same as those used for amplification.

Population genetic and phylogeographic analysis

The chromatograms of each sequence were checked visually using Chromas ver. 2.33 (http://www.technelysium.com.au). DNA sequences were aligned in Clustal_X (Thompson et al., 1997) and then checked by visual inspection with minor adjustment. The cpDNA haplotypes were identified by DnaSP ver. 5.0 (Librado & Rozas, 2009). All newly generated sequences from T. himalayanum were submitted to GenBank, and the access numbers are MH122760–MH122783. Three cpDNA fragments, rbcL-accD, rps15-ycf1, and trnH-psbA, were then concatenated for the subsequent analyses. In the following analysis, we did not code the indels; each indel was considered as a one-step mutation.

The spatial genetic structure of cpDNA haplotypes was analyzed to define groups of populations that were geographically homogeneous and maximally differentiated from each other using the program SAMOVA ver. 1.0 (Dupanloup, Schneider & Excoffier, 2002); the number of initial conditions was set to 100, and pairwise difference (DNA) was chosen as the molecular distance.

Estimates of average gene diversity within populations (H_S), total gene diversity (H_T), G_ST (Nei, 1987), and N_ST (Grivet & Petit, 2002) were calculated for the overall populations and for each regional group using PERMUT ver. 2.0 (Pons & Petit, 1996). G_ST and N_ST were used to estimate population differentiations; G_ST only takes haplotype frequency into account, whereas N_ST also takes sequence similarities between haplotypes into account. A permutation test with 1,000 permutations was used to compare G_ST and N_ST. Significantly higher values of N_ST than the estimated G_ST indicated the presence of a phylogeographical structure (Pons & Petit, 1996).

Analysis of molecular variance (AMOVA; Excoffier, Smouse & Quattro, 1992) for overall populations and each regional group of T. himalayanum was conducted using ARLEQUIN version 3.11 (Excoffier, Laval & Schneider, 2005). The F-statistic (F_ST) was calculated, and significance was tested.

Phylogenetic analyses

The appropriate model of DNA substitution for maximum likelihood (ML) analysis was determined using the AIC in jModelTest ver. 2.1.6 (Darriba et al., 2012; Guindon & Gascuel, 2003). Best-fit models for the cpDNA datasets were determined to be HKY+G. ML trees of cpDNA haplotypes of T. himalayanum were generated in PhyML ver. 3.1 (Guindon & Gascuel, 2003) using L. trichosantha as the outgroup. We selected BioNJ as the distance-based tree reconstruction method during tree searching. Subtree pruning and regrafting (SPR) moves were relied on to explore the space of tree topologies. The number of random starting trees was five. The tree branch support was estimated by non-parametric bootstrap analysis, and the number of bootstrap replicates was 100. Bayesian inference was conducted using MrBayes ver. 3.2 (Ronquist et al., 2012). The MCMC algorithm was run for 10,000,000 generations starting from random trees, and one in every 100 generations was sampled. The first 25% of the trees was discarded as burn-in followed by reconstruction of the majority consensus tree. Internodes with posterior probabilities >95% were considered to be statistically significant. When the average standard deviation of split frequencies between the two independent runs was smaller than 0.01, the analysis was stopped. Maximum parsimony (MP) analysis was implemented using Phylip ver. 3.69 (Felsenstein, 1985, 1988; Margush & McMorris, 1981). The number of bootstrap replicates was set to 1,000, and the random number seed was five.

To further detect genealogical relationships, a haplotype network was constructed among cpDNA haplotypes using the median‐joining method and MP calculations as implemented in Network ver. 5.0.0.0 (Bandelt, Forster & Röhl, 1999; Polzin & Daneschmand, 2003; http://www.fluxus-engineering.com).

Demographic analyses

For the overall populations and three regional groups, Tajima’s D and Fu’s Fs of neutrality tests with 1,000 simulated samples were performed using ARLEQUIN ver. 3.11 (Fu, 1997). Negative values are expected when there have been recent population expansions or population bottlenecks (Tajima, 1989). Pairwise mismatch distribution analysis (Rogers & Harpending, 1992) for overall populations and each group based on the cpDNA dataset was performed assuming historical population dynamics using DnaSP ver. 5.0 (Librado & Rozas, 2009). Unimodal distributions reflect rapid growth from small populations, while multimodal distributions reflect long-term population stability (Slatkin & Hudson, 1991; Rogers & Harpending, 1992).

Palaeodistributional reconstruction

Potential geographic ranges and the effects of past climatic oscillations for T. himalayanum has computed across current, Mid Holocene (6 ka), last glacial maximum (LGM; 22 ka), and last interglacial maximum (LIG; 130 ka, Otto-Bliesner et al., 2006). We simulated all the species distribution models, i.e., GBM (generalized boosted models; Ridgeway, 2006), SRE (surface range envelop; Busby, 1991), GLM (generalized linear models; Peter & John, 1989), CTA (classification tree analysis; Breiman, 2001, 2017), ANN (artificial neural network; Ripley, 1996), FDA (flexible discriminant analysis; Hastie, Tibshirani & Buja, 1994), MARS (multivariate adaptive regression splines; Friedman, 1991), RF (random forests; Breiman, 2001), and MAXENT (maximum entropy; Phillips, Dudík & Schapire, 2004) in R packages “Biomod 2” (Thuiller, Georges & Engler, 2014), supported with rworld map, rgdal, dismo, SDMTools. The effectiveness of all the models were evaluated using TSS (true skill statistics) and AUC (area under curve; curve-receiver operating characteristic curve) values >0.7. A total of nine bioclimatic variables (BIO4 = temperature seasonality, BIO5 = max temperature of warmest month, BIO6 = min temperature of coldest month, BIO7 = temperature annual range, BIO8 = mean temperature of wettest quarter, BIO9 = mean temperature of driest quarter, BIO10 = mean temperature of warmest quarter, BIO11 = mean temperature of coldest quarter, and BIO15 = precipitation seasonality) showed low correlation and high informativeness after a jackknife procedure on the 19 BIOCLIM variables downloaded from the WorldClim data set (Hijmans et al., 2005).

Haplotype variation and distribution

The combined cpDNA sequences of rbcL-accD, rps15-ycf1, and trnH-psbA of T. himalayanum ranged from 1,698 to 1,718 bp in length, with an alignment length of 1,727 bp. Nineteen haplotypes (H1–H19) were identified based on polymorphic sites across the concatenated cpDNA sequence data (Table 2). Among those polymorphic sites, we found 23 single-site mutations and eight indels. The rbcL-accD region was the most variable fragment (14 substitutions; three indels), followed by rps15–ycf1 (six substitutions; five indels), and trnH-psbA (three substitutions).

Twelve haplotypes (accounting for 63.2% of the total haplotypes) were relatively rare, each was restricted to a single population (H1, H3–H6, H9–H13, H15, and H17). Moreover, all of the twelve haplotypes were distributed in the HM group of five populations, indicating a very high level of molecular diversity (Table 1). Another seven haplotypes were shared by at least two populations, suggesting a substantial amount of population isolation at the cpDNA level.

The most common haplotype was H14, which occurred in nine populations with variable locations in the HM group (three populations) and NHM group (six populations). Haplotypes H7 and H16 both occurred in four populations but with locations not far from each other (Fig. 1).

Population variation and phylogeographic structure

Twelve populations (P1–P6, P13–P15, and P18–P20) were fixed for single haplotype, showing no variation within the population. However, the remaining eight populations (P7–P12, P16, and P17) in the HM and NHM groups were particularly diverse, with two to five haplotypes in each population (Table 1; Fig. 1). Populations P1–P6 were found to belong to the EH group; each contained only one haplotype, and all three haplotypes present in these six populations were endemic.

Population P19 was located far from P20, isolated by the Sichuan Basin. However, these two populations shared the sole haplotype H19.

SAMOVA failed to detect any meaningful geographical groups based on cpDNA haplotypes. K values increased from two to seven while F_CT values changed irregularly from 0.897 to 0.781. However, each newly defined group was represented by a single population. On the basis of the distribution ecoregions of the populations and the genetic structure of the T. himalayanum populations, we divided the populations into three groups: the East Himalayan group (EH group, including populations P1–P7), the Hengduan Mountains group (HM group, P8–P12), and the North Hengduan Mountains group (NHM group, including P13–P20). The division backgrounds which can be considered as temporal and biota differences between these three eco-regions are illustrated as follows: (1) distinctions between the Hengduan Mountains and the Himalaya. Mulch & Chamberlain (2006) suggests that the orogeny of the Hengduan Mountains occurred as a final propagation of the palaeoelevation history of the QTP uplift after 10 Ma, reflecting its much younger age than the Himalaya. In addition, the HM is more species-rich compared with the Himalaya, and biota diversification process of HM is distinguishing from that of Himalaya (Xing & Ree, 2017); (2) differences between the NHM group and HM group-even the NHM group occupies part of the Hengduan Mountains and adjoining highland, the genetic diversity of T. himalayanum and some other species (Du et al., 2017) in there is quite lower than that of core area of Hengduan Mountains (HM group). A similar grouping strategy can be found in many studies of plant species that are sympatric with T. himalayanum (Sun & Zhou, 1996; Zou et al., 2012; Liu et al., 2013; Yu et al., 2015; Gao et al., 2016; Du et al., 2017). The EH group was distributed along the Yarlung Zangbo River or Himalayan Mountains, ranging from the big bend gorge of Yalutsangpu to the Yadong valley. The HM group is located in the core region of the Hengduan Mountains with rich molecular diversity. The location of the NHM group is the northern part of the Hengduan Mountains region, which is also the eastern edge of the QTP.

For the overall populations, the N_ST value was 0.929, while the G_ST value was 0.781, and the difference between them was significant (p < 0.001). A significant difference between N_ST and G_ST was also found in the EH and NHM groups but not in the HM group. The values of total genetic diversity (H_T) were much higher than those of gene diversity within populations (H_S) for T. himalayanum in the overall populations (H_T = 0.878, H_S = 0.192), EH group (H_T = 0.714, H_S = 0.075), HM group (H_T = 0.922, H_S = 0.550), and NHM group (H_T = 0.498, H_S = 0.071, Table 3). AMOVA for the overall populations revealed that 89.43% of the total genetic variation occurred among populations, compared with 10.57% within populations; the pairwise F_ST value was 0.8943. Similar AMOVA results were obtained for both the EH group and NHM group. However, the HM group exhibited a higher level of within-population variation compared with among-population variation (Table 4).

Phylogenetic relationship

The relationships of the 19 haplotypes are shown in Figs. 2 and 3. Phylogenetic trees constructed using MP, ML, and Bayesian methods are not consistent with each other. They do not show clear relationships between cpDNA haplotypes. However, H19 was observed to occupy the radical position in all trees, and H5 and H6 clustered as one clade in the basal position but only in the Bayesian inference and ML trees. The brush structure observed in the Bayesian tree may be due to the rapid diversification of this species. In the network phylogenetic analysis, the number of mutations between H6 and H7 was seven, illustrating considerable differentiation between the EH group and HM group.

Demographic history

Tajima’s D values of the overall populations and the HM group were negative but not significant (p > 0.05), suggesting slight range expansion. Fu’s Fs statistics showed positive values or not significant negative values for the whole gene pool and three groups (Table 5). Mismatch distributions of the overall populations, the EH group and the NHM group were multimodal, indicating a demographic equilibrium; the HM was unimodal (Fig. 4). A star-like pattern network is considered to indicate a demographic expansion (Hudson, 1990), thus supporting the findings of the Tajima’s D. H7 and H14 are more likely to be the ancestral haplotypes because they occupy a relative central position in the network (Carbone & Kohn, 2001).

Paleaodistributional reconstruction

The Palaeodistributional reconstruction showed that during the last interglacial maximum T. himalayanum occupied large scope of Himalaya, almost the whole of Hengduan Mountains area, and adjoining region arounding the Sichuan Basin (Fig. 5). During the last glacial maximum, the species was restricted southwards with sharp reduction of distribution, and there were few distribution in East Himalaya and North Hengduan Mountains. From Mid Holocene up to the current the species is expanding Northwards and Westwards again. Climatic niche during three different times suggesting that this alpine herb went through repeated retreat and expansion, and affected by the severe climatic conditions of the QTP.

Phylogeographic structure of T. himalayanum and its correlation with heterogeneous habitats and Quaternary climatic oscillation

One characteristic of the geographic genetic structure of T. himalayanum is a high level of private haplotypes. Twelve out of 19 detected cpDNA haplotypes were private. High total diversity versus a relatively low within-population diversity is another feature of T. himalayanum, and this type of geographic genetic structure is found in many studied alpine species located in the HHM and adjoining regions (Yang et al., 2008; Wang et al., 2009; Zhang, Meng & Rao, 2014; Gao et al., 2016; Wan et al., 2016). The comparisons between the values of H_T and those of H_S in the overall populations and the three regional groups suggest high genetic diversity between populations. The AMOVA results also showed high levels of genetic differentiation between populations. These findings indicate a clear phylogeographic structure (Avise, 2004) except for the HM group (with an among-population variation of 39.66%). Factors promoting the among-population differentiation of T. himalayanum are likely to be complicated. Among those, fragmented distribution and population isolation caused by historical orogenesis and climatic oscillations could be the main factors. In heterogeneous habitats, diverse elevations, precipitation, temperatures, edaphic factors, wind systems, and ecological selection resulting from competitors, predators, parasites, and perhaps mutualists that might restrict other species’ ranges, could accelerate species differentiation (Renner, 2016; Chozas et al., 2017).

Orogenic activities cause a large geomorphologic adjustment, and hydrologic, and tectonic configurations are produced (Li & Fang, 1999). Studies on palaeoaltimetry of the QTP show a progressive uplift from south to north which is likely to have started during or shortly after the beginning of the collision (about 45–40 Ma; Li & Fang, 1999; Mulch & Chamberlain, 2006). Since then, the QTP has experienced shortening on its north–south axis accompanied by the formation of rough surface features and clockwise rotation in its eastern part (Favre et al., 2015). These result in the formation of north–south oriented valleys/rivers in the southeast plateau (including the Yangtze River, the Langcang River, and the Nujiang valleys), as well as west–east oriented valleys/rivers in the south (e.g., the Yarlung Zangbo River) surrounded by high peaks (Favre et al., 2015), which provide pathways for south–north and east–west migration (Wu, 1979; Yu et al., 2017). Due to its immense size and altitude, the QTP and adjacent highland is closely linked to development of the Asian monsoon system and consequently contributed to climatic changes in this region (Ruddiman & Kutzbach, 1991; Liu & Yin, 2002; Tang et al., 2013). During the glacial-interglacial cycle, many species experienced divergent or “contact” evolution during glacial contractions and postglacial expansions (Avise, 2000; Yu et al., 2015). Based on our survey, the distribution of diverse haplotypes in the HM group may have resulted from several contraction/expansion cycles following the repeated advance/retreat of glaciations. The contraction/expansion process may also result in the accumulation of rich haplotypes within the population. This region may play a diversification role and may even be the origin center of T. himalayanum. Successive founder events during postglacial recolonization would have led to reduced genomic variability and areas of genetic homogeneity along the recolonization route (Hewitt, 2004). This phenomenon has been frequently observed in the HHM biota (Opgenoorth et al., 2010; Gao et al., 2016; Yu et al., 2017, and references therein). Thus, the decline of gene diversity from the HM group (H_T = 0.922) to the NHM group (H_T = 0.498) could result from the south–north direction postglacial recolonization. This hypothesis was also best supported by the palaeodistributional reconstruction (Fig. 5). Furthermore, complicated topography produced by tectonic movement and Quaternary glacial-interglacial cycles could result in range fragmentation of T. himalayanum, leading to reduced gene flow and the initiation of intra-species allopatric divergence (Coyne, 1992; Rice & Hostert, 1993; Hewitt, 2000, 2004). This may explain that eight of the 12 haplotypes in the HM group and four of the five haplotypes in the EH group are private. Most studies on the alpine flora of the HHM and adjacent regions have highlighted the effect of the Quaternary period on intraspecific diversification (Zhang, Volis & Sun, 2010; Gao et al., 2012b; Zhang et al., 2012b; Wang et al., 2013; Li et al., 2013; Khan et al., 2015). Therefore, the high genetic variations among populations of T. himalayanum could result from strong bottlenecks and founder effects/genetic drift to fix different alleles in isolated populations of fragmented regions (Birky & William, 2001).

Demographic history of triosteum himalayanum in the HHM area

Even the values of Fu’s Fs and mismatch distributions for the overall populations (Table 5) reject the model of current extensive expansion, the star-like network together with the negative values of Tajima’s D for the overall populations indicate a demographic expansion. Additionally, the palaeodistributional reconstruction revealed a clear pattern of recolonization from the Mid Holocene up to current. Therefore, we prefer the indication that T. himalayanum did experience range expansion from the core region of Hengduan Mountains (HM group) to both the East Himalaya area (EH group) and the North Hengduan Mountains (NHM group) from Mid Holocene.

The phylogenic relationship of the star‐like network indicates that haplotypes of the EH group originated from Hengduan Mountains areas, and there was also a decline of gene diversity from the HM group (H_T = 0.922) to the EH group (H_T = 0.714). One of the migration corridors may be the Yarlung Zangbo Valley, which is the largest vapor channel of the QTP (Lin & Wu, 1990). Lower elevation and intense airflow create ideal circumstances for the dispersal of seed and pollen of T. himalayanum. The other important dispersal corridors could be the long, narrow mountain areas. Since the Himalayas are linked to the Hengduan Mountains, exchanges along the mountain ridges or slopes occur easily. This hypothesis has been assumed in many studies and reviews (Cun & Wang, 2010; Luebert & Weigend, 2014; Yu et al., 2015, 2017). As discussed above, the south–north direction of the Hengduan Mountains and the river valley provide channels for glacial contractions and postglacial expansions for plant and animals and facilitate communication between different ecoregions of the HM and NHM areas. Therefore, haplotype H14 is not only shared by both the HM group and the NHM group, but it is also the most dominant and widely distributed haplotype.

However, anthropogenic activities that reduce the distribution range and genetic diversity could be other factors responsible for lower gene diversity and within-population variation of the NHM group in addition to bottleneck or founder effects. Many previous habitats of T. himalayanum recorded in documents are now surrounded by folk houses, and we were not able to locate them during the field work conducted two years ago. We should pay more attention to this phenomenon and more effectively preserve this alpine herb for future study and utilization.