Population sampling
During the summers of 2008 through 2009, we collected the samples throughout the range of S.tetraptera.
Fresh leaves were collected from 35 populations and, with few exceptions, 4–12 individuals that were at least 50 m apart from each other were sampled from each population (Table 1, Fig. 1). We measured the longitude, latitude and altitude of each collection location using an Etrex Global Positioning System (Garmin). In total, 301 individuals were sampled and leaves were dried with silica gel. All samples were identified by Professor Guoying Zhou. Voucher specimens of all populations were deposited at the herbarium of the QTPMB (Qinghai-Tibetan Plateau Museum of biology), Xining, Qinghai Province, China.
Table 1
Geographic distributions, gene diversity, nucleotide diversity, and haplotype frequencies of cpDNA sequences for Swertia tetraptera
|
Population
code
|
Longitude
|
latitude
|
Sample
number
|
Hd
|
Pi
|
Haplotype
(no. of individuals)
|
Qilian, Qinghai
|
QL
|
100.27ºE
|
38.16ºN
|
8
|
0.67857
|
0.00068
|
H1(4), H4(3), H6(1)
|
Mengyuanqingshizui, Qinghai
|
MY1
|
101.41ºE
|
37.46ºN
|
7
|
0.80952
|
0.00134
|
H3(2), H4(3), H7(1), H10(1)
|
Mengyuanxianmi, Qinghai
|
MY2
|
102.01ºE
|
37.28ºN
|
4
|
0.50000
|
0.00083
|
H24(1), H54(3)
|
Mengyuanxianmi, Qinghai
|
MY3
|
102.02ºE
|
37.38ºN
|
11
|
0.80000
|
0.00321
|
H7(2), H16(3), H53(4), H54(2)
|
Huzhu, Qingahi
|
HZ
|
102.85 ºE
|
37.03ºN
|
11
|
0.18182
|
0.00015
|
H3(1), H4(10)
|
Gangcha, Qinghai
|
GC
|
101.18ºE
|
37.25ºN
|
7
|
0.71429
|
0.00291
|
H4(4), H16(1), H32(1), H42(1)
|
Datong, Qinghai
|
DT
|
101.53ºE
|
37.22ºN
|
9
|
0.80556
|
0.00101
|
H1(2), H2(1), H3 (1), H4 (4), H5(1)
|
Huangyuan, Qinghai
|
HY
|
101.37ºE
|
36.73ºN
|
10
|
0.86667
|
0.00106
|
H1(1), H3 (1), H4 (3), H6(3), H7(1), H8(1)
|
Huangzhong, Qinghai
|
HZH
|
101.70ºE
|
36.30ºN
|
7
|
0.66667
|
0.00087
|
H1(4), H4 (1), H6(2)
|
Pingan, Qinghai
|
PA
|
101.91ºE
|
36.32ºN
|
12
|
0.86364
|
0.00170
|
H1(1),H2 (1),H3 (4),H4(3),H7(1),H9(1),H10(1)
|
Hualong, Qinghai
|
HL
|
101.97ºE
|
36.27ºN
|
11
|
0.87273
|
0.00183
|
H1(4),H2(2),H6(1),H7(1),H11(1),H12(1),H13(1)
|
Ledu, Qinghai
|
LD
|
102.39ºE
|
36.66ºN
|
8
|
0.25000
|
0.00083
|
H1(7),H14(1)
|
Minghe, Qinghai
|
MH
|
102.69ºE
|
36.09ºN
|
9
|
0.94444
|
0.00381
|
H1(1),H4(4),H8(1),H13(2),H15(1),H16(2),H17(1)
|
Xunhua, Qinghai
|
XH
|
102.31ºE
|
35.74ºN
|
10
|
0.77778
|
0.00184
|
H1(1),H15(5),H16(1),H18(1),H19(1),H20(1)
|
Tongre, Qinghai
|
TR
|
101.94ºE
|
35.31ºN
|
4
|
0.50000
|
0.00041
|
H4(1),H6(3)
|
Zeku, Qinghai
|
ZK
|
101.87ºE
|
35.25ºN
|
9
|
0.69444
|
0.00073
|
H3(1),H4(1),H6(2),H7(5)
|
Henan, Qinghai
|
HN
|
101.55ºE
|
34.60ºN
|
6
|
0.80000
|
0.00138
|
H3(1),H4(3),H6(1),H21(1)
|
Maqing, Qinghai
|
MQ
|
100.21ºE
|
34.47ºN
|
9
|
0.83333
|
0.00128
|
H1(1),H3(1),H4(4),H6(1),H7(1),H9(1)
|
Jiuzhi, Qinghai
|
JZ
|
101.37ºE
|
33.39ºN
|
10
|
0.73333
|
0.00202
|
H4(5),H5(2),H22(2),H23(1)
|
Dari, Qinghai
|
DR
|
99.60ºE
|
33.75ºN
|
10
|
0.66667
|
0.00066
|
H1(1),H4(6),H6(1),H24(1),H25(1)
|
Banma, Qinghai
|
BM
|
100.71ºE
|
33.05ºN
|
10
|
0.88889
|
0.00229
|
H1(2),H3(3),H4(2),H24(1),H26(1),H27(1)
|
Yushu, Qinghai
|
YS
|
96.73ºE
|
33.13ºN
|
7
|
0.90476
|
0.00173
|
H3(2),H7(1),H28(1),H29(2),H30(1)
|
Chenduo, Qinghai
|
CD
|
97.44ºE
|
33.11ºN
|
8
|
0.96429
|
0.00357
|
H7(1),H12(1),H15(1),H19(1),H22(1),H31(1),H32(2)
|
Nangqian, Qinghai
|
NQ
|
96.14ºE
|
32.44ºN
|
11
|
0.98182
|
0.00316
|
H3(1),H6(1),H9(1),H12(1),H13(1),H32(1),H33(1),
H34(2),H35(1),H36(1)
|
Tianzhu, Gansu
|
TZ
|
102.60ºE
|
36.94ºN
|
8
|
0.92857
|
0.00177
|
H3(2),H7(1),H30(1),H37(1),H38(2),H39(1)
|
Xiahe, Gansu
|
XHE
|
102.47ºE
|
35.15ºN
|
11
|
0.96364
|
0.00201
|
H1(2),H3(1),H4(2),H7(1),H13(1),H40(1),H41(1),
H42(1),H43(1)
|
Hezuo, Gansu
|
HZU
|
102.91ºE
|
35.06ºN
|
9
|
0.91667
|
0.00142
|
H1(2),H2(1),H3 (2),H4(2),H7(1),H44(1)
|
Maqu, Gansu
|
MQU
|
102.64ºE
|
34.10ºN
|
10
|
0.95556
|
0.00209
|
H4(2),H5(1),H24(1),H30(1),H35(1),H45(1),H46(2),
H47(1)
|
Diebu, Gansu
|
DB
|
103.89ºE
|
34.23ºN
|
9
|
0.97222
|
0.00362
|
H1(1),H4(1),H7(1),H14(1),H16(2),H26(1),H31(1),
H48(1)
|
Shiqu, Sichuan
|
SQ
|
97.52ºE
|
33.15ºN
|
9
|
0.69444
|
0.00092
|
H1(2),H4(5),H24(1),H49(1)
|
Ruoergaijiangzha, Sichuan
|
REG1
|
103.18ºE
|
33.60ºN
|
10
|
0.86667
|
0.00198
|
H3(3),H4(3),H11(1),H21(1),H34(1),H50(1)
|
Ruoergaibasi, Sichuan
|
REG2
|
102.78ºE
|
34.14ºN
|
9
|
0.91667
|
0.00229
|
H1(1),H4(3),H7(1),H36(1),H42(1),H51(1),H52(1)
|
Ruoergaibanyou, Sichuan
|
REG3
|
103.10ºE
|
33.59ºN
|
9
|
0.86111
|
0.00335
|
H1(2),H2(2),H16(3),H44(1),H46(1)
|
Aba, Sichuan
|
AB
|
101.85ºE
|
32.92ºN
|
5
|
0.70000
|
0.00149
|
H14(1),H16(1),H53(3)
|
Hongyuan, Sichuan
|
HOY
|
102.33ºE
|
32.65ºN
|
4
|
0.50000
|
0.00083
|
H24(1),H54(3)
|
Total
|
|
|
|
301
|
0.90800
|
0.00257
|
|
DNA extraction, amplification and sequencing
Total genomic DNA was extracted from silica gel-dried leaves using the modified CTAB method [31] and used as template in the polymerase chain reaction (PCR). A preliminary universal primer scans of chloroplast DNA genomes were performed on 10 individuals from 10 different populations using 5 pairs of primers. Primers a and f of Taberlet et al. (1991) [32] were used to amplify trnT-trnF region and sequenced with primers a, c and f. The other four regions psbB-psbH, rpl20–5′rps12, trnS-trnG and psbA-trnH were amplified and sequenced using primers described in Hamilton (1999) [33]. Primer used to amplify trnS-trnG, showed different sequences within the 10 individuals tested, and were used thereafter for the large-scale survey of haplotype variation within S.tetraptera. Polymerase chain reaction (PCR) was performed in a 25µLvolume, containing 10–40 ng plant DNA, 50 mMTris-HCl, 1.5 mM MgCl2, 250µg/mL BSA, 0.5 mM dNTPs, 2µMof each primer, and 1 unit of Taq polymerase. Initial template denaturation was programmed at 95°C for 5 min, followed by 35 cycles of 94°C for 1 min, 52°C for 1 min, and 72°C for 1 min plus a final extension of 72°C for 7 min. The BioDev Gel Extraction System B Kit (BioDev-Tech) was used to purify all successfully amplified DNA fragments. The PCR primers trnS and trnG were adopted to perform sequencing reactions using ABI Prism Bigdye™ Terminator Cycle Sequencing Ready Reaction Kit. The all trnL-trnF sequences of S.tetraptera used in this study were quoted from Yang et al.,(2011)[30].
Proofreading and Alignment of DNA, and Data Analysis
We used BioEdit v 7.0.9.0 [34] software for manual proofreading and examining the variable sites. And then we used MEGA v 7.0 [35] software to get rid of low quality sequences and we only used high quality sequences for subsequent analysis. The sequencematrix-1.8 was used to combine the trnL-trnF and trnS-trnG gene segments, and the combined cpDNA gene segment was used to carry out the subsequent analysis.
Genetic Variation and Genetic Structure Analysis
We used the program ARLEQUIN version 3.11 [36] to calculate the indices of unbiased genetic diversity (Hd) and nucleotide diversity (π) for each population. Differentiations among populations and within populations were calculated by analyses of molecular variance (AMOVA) using ARLEQUIN version 3.11.
The program PERMUT was used to calculate total gene diversity (HT), average gene diversity within populations (HS), GST (coefficient of genetic variation over all populations) and NST (equivalent coefficient taking into account sequences similarities between haplotypes) [37]. A comparison was made between GST and NST using a permutation test with 1000 permutations. The phylogeographic structure appear when the NST value is significantly larger than GST value.
Phylogenetic Analysis
We used DnaSP v5.10 to identify haplotypes of the cpDNA genes. Haplotype distribution was visualized by GenGIS [38]. The program NETWORK 4.5.1.0 (available at http:⁄ ⁄www.fluxus-engineering.com) was used to construct the haplotype network based on Median Joining (MJ) method [39] and MP calculation [40].
In this study, phylogenetic tree was constructed based on Bayesian inference (BI) [41] and maximum parsimony (MP). Mafft 7.205 software was used to compare the cpDNA sequences and remove the irregular sequences at both ends [42]. Before building BI tree, PAUP and MrModeltest are jointly run through MrMTgui. AkaikeInformation Criterion (AIC) results showed that the best model for BI analysis is GTR + I + G, with random tree as starting tree. Start with four Markov chains, namely three hot chains and one cold chain, save one tree every 100 generations, calculate 9 000 000 generations in total, discard the first 25% preheated (Burin-in) tree, and use the remaining tree to calculate the Bayesian posterior probability of the consistent tree and each branch (PP, Posteriorpssibility). PAUP 4.0b10 software was used to construct MP tree [43].
Demographic history
Mismatch distribution and neutral test for the existing populations of S.tetraptera were conducted with DnaSP Ver.5.10 [44] software. If the result of mismatch distribution is unimodal, it indicates that the population may have experienced recent expansion. If it is multiple peaks, it means that the population size is relatively stable and in individual equilibrium in a long time [45]. For the neutral test, two infinite mutation-site model indices such as Tajima’s D, Fu's Fs [46–48] were selected to predict the nature of sequence evolution and possible population history dynamics. Negative values of two indices indicate that the population may have undergone recent expansion or selective sweep. Positive values of two indices indicate that populations may have been geographically isolated for a long time and mutation differences between populations were accumulated or controlled by balanced selection [45, 49]. Arlequin Ver. 3.5.2 [36] was used to test the results of mismatch distribution analysis, among which SSD (sum of square deviation) was used to test whether to accept the hypothesis of population rapid expansion recently.
Species Distribution Modeling
A total of 97 sites were collected in this study, covering the known distribution areas of S.tetraptera. The geographic distribution data of S.tetraptera were obtained mainly by the following methods: (1) field investigation (for detailed information, see Table 1); (2) the network data, including the China digital plant herbarium (http://www.cvh.org.cn/), China plants subject database (http://www.plant.csdb.cn/) and the Chinese image library (http://www.plantphoto.cn/). (3) Literature review, including Chinese and English journals, flora of China, flora of local areas, investigation reports of nature reserves, etc. After removing duplicate records using ENMTools program from the same locality which can reduce the influences of autocorrelation, 82 sites of S.tetraptera were applied for the ecological niche modeling (ENM).
Potential distributions of the Last Glacial Maximum (LGM, 0.021–0.018 Ma), current and future (year 2050) of S.tetraptera was predicted using MaxEnt v3.4.1 [50–51]. Three general circulation models (BCC-CSM2-MR, CNRM-CM6-1 & MIROC-ES2L) and four shared socio-economic pathways (SSP1-2.6, SSP2-4.5, SSP3-7.0 & SSP5-8.5) were selected for 12 sets of climate simulation data for the future decade. As for as LGM period, paleo climatic layers simulated by the New Earth System Model of the Max Planck Institute for Meteorology, the Community Climate System Model Version Version 4 (CCSM4; [52]), and the Model for Interdisciplinary Research on Climate Earth System Model (MIROC-ESM; [53]). We downloaded layers for 19 bioclimatic variables (Table S2) of these models plus for the current time (1970–2000) at 2.5 arc-min resolution from the WorldClim website (www.worldclim.org) for the study area. There is a certain correlation between 19 bioclimatic variables (Peterson et al., 2011). Pearson correlation coefficient matrix among 19 contemporary bioclimatic variables was calculated using ENMTools version 1.4.4 [53], and 0.80 was used as the threshold to determine whether the correlation was significant. MaxEnt software version 3.4.1was used for pre-modeling to obtain the percentage contribution of each factor to the model and the analysis result of the Jackknife. By retaining the relatively important factor among the significant correlation factors, the combination of climate factors for model construction was determined [54]. The eventual chose variables that were applied to find changes in the distribution ranges of S.tetraptera were bio4, bio10, bio11, bio14, and bio15.
Feature class and regularization multiplier optimization realized using R program and Kuenm package [55]. Select the minimum value of AICc as the optimal setting and establish the model [54, 56]. LQ and 0.1 were selected as the optimal setting for feature class and regularization multiplier.
The distribution data of S.tetraptera and corresponding bioclimatic variables in each period were imported into Maxent V3.4.1. In parameter setting, the Maximum iterations is set to 5000, and the Subsample method is adopted to run 10 times repeatedly. The default parameters of other parameters were used for MaxEnt, and the final distribution model is obtained by taking average values. The output format of model analysis results is ASCII grid layer, and the value of fitness index is between 0 and 1. Receiver Operating characteristic Curve (ROC) was used to analyze the prediction accuracy. The greater the Area under the ROC curve (AUC), the higher the prediction accuracy of the model [57]. The calculation results of MaxEnt V3.4.1 were imported into DIVA-GIS V7.5, and the Mongolian map layer made by Chinese map was used to limit the analysis scope to China. As threshold rule, we selected the maximum test sensitivity plus specificity logistic threshold, which is very robust with all types of data [58].