Intraspecific divergences and phylogeography of Panzerina lanata (Lamiaceae) in northwest China

Population sampling

Across northwest China, there are approximately 27 natural populations of P. lanata, from which we collected 269 individuals (Table 1). All individuals and 200 individuals from 27 natural populations were used for chloroplast analysis and nuclear gene analysis, respectively. In each population, 4–14 individuals were collected. To avoid duplicated sampling, we set 30 m as the minimum distance between individuals within each population. Fresh leaves were harvested and dehydrated using silica gel for later DNA isolation. Two specimens were collected from each population, and the voucher specimens were deposited in the Herbarium of Xinjiang Institute of Ecology and Geography, Chinese Academy of Science (XJBI). Leonurus turkestanicus and Lagochilus ilicifolius were included as outgroups in this analysis (Meng & Zhang, 2013).

DNA isolation, PCR amplification, and sequencing

CTAB method was used for genomic DNA isolation (Doyle & Doyle, 1987). The cpDNA intergenic spacers, trnH-psbA and rpoB-trnC (Sang, Crawford & Stuessy, 1997; Shaw, Lickey & Beck, 2005), and a nuclear ribosomal region (ITS1-ITS4) were amplified. The PCR reactions (25 µL) contained 10 ×PCR buffer (2.5 µL), MgCl₂ (25 mM, 2.5 µL), dNTP mixture (2.5 mM, 2.0 µL), forward primer (1 µL) and reverse primer (1 µL), Taq polymerase (0.125 µL) (Takara, Kusatsu, Japan), and genomic DNA (1 µL) as template. The amplification cycles were as following: initial denaturation (95 °C, 4 min), 36 cycles of denaturation (94 °C, 30s), annealing (52 °C, 30s) and extension (72 °C, 1 min), and a final extension (72 °C, 10 min). The PCR products were examined by electrophoresis using agarose gel (1%), isolated with a QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany), and sequenced on ABI Prism 3730 Genetic Analyzer (Sangon, Shanghai, China).

For DNA sequence editing, we used SeqMan (Lasergene, USA), and used BioEdit (Hall, 1999) for alignment. All sequences are deposited in GenBank under accession numbers MK299519–MK299550 for ITS, MK299551–MK299556 for trnH-psbA, MK299557–MK299559 for rpoB-trnC.

Genetic diversity and population structure

Arlequin 3.5 was applied to determine haplotype diversity (H) and nucleotide diversity (π) (Excoffier & Lischer, 2010). For subdivision of their geographical structure, we used SAMOVA v1.0 (Dupanloup, Schneider & Excoffier, 2002). We set 2 ≤ K ≤ 12 until the F_CT values reached the maximum and, when a single population was clustered into one group, the combination was excluded (Beatty, Provan & Comes, 2013; Iwasaki et al., 2012). With Arlequin 3.5, AMOVA estimated genetic variation on three different levels: inter-groups, inter-populations within groups, and intra-populations. Calculation of genetic differentiation was determined using a significance test based on 1,000 permutations. To identify genetic differentiation (G_ST, N_ST), average heterozygosity within populations (H_S) and across total populations (H_T), we used Permut v1.0 (Pons & Petit, 1996) for 1,000 permutation tests. H_T and H_S were used to estimate genetic diversity. These two parameters (N_ST, G_ST) were used to test if phylogeographic structure existed. By applying the median-joining algorithm, we also implemented the phylogenetic relationship among the haplotypes Network v5.0 (Bandelt, Forster & Röhl, 1999).

Demographic history and divergence time analyses

In order to test whether all populations and groups of P. lanata divided by the SAMOVA experienced demographic expansion, we used Arlequin 3.5 to determine Tajima’s D (Tajima, 1989) and Fu’s F_S (Fu, 1997) with 1,000 permutation tests. We further applied a mismatch distribution analysis to measure population expansion for all groups. In general, unimodal distribution patterns showed a recent expansion event, while the distribution pattern was stable according to bimodal and multimodal results. We also estimated the sum of squared deviations (SSD) and the raggedness index of Harpending (HRag) (Harpending, 1994) between the observed and expected mismatches in Arlequin 3.5 with 1,000 permutation tests. The significance of P was examined if the populations experienced expansion.

We estimated the divergence time between different lineages of P. lanata using Beast v1.6.1 (Drummond & Rambaut, 2007). As most angiosperms have a cpDNA nucleotide substitution rate of 1.0–3. 0 × 10⁻⁹ s/s/y, we use 2. 0 × 10⁻⁹ s/s/y with a SD of 6. 080 × 10⁻¹⁰ s/s/y to estimate the divergence time of these two cpDNA regions (trnH-psbA and rpoB-trnC) (Jia et al., 2012; Wang et al., 2016). We ran the GTR+G substitution model and the MCMC chains for 10,000,000 generations and sampled every 1,000 generations. The effective population sizes (ESS > 200) were checked in Trace v1.5. We used TreeAnnotator v1.6.1 to implement the maximum clade credibility tree with a burn-in of 1,000 trees. For tree editing, we applied FigTree v1.3.1.

Species distribution inference

For estimation of the past (Last Glacial Maximum, 21,000 before present) and present distribution ranges of P. lanata, we used the ecological niche model to infer the potential distribution area. We simulated the potential distribution with the maximum entropy method in MAXENT 3.3.1 (Phillips, Anderson & Schapire, 2006). Geographical information of P. lanata was gleaned from field collection records and the Chinese Virtual Herbarium (http://www.cvh.org.cn/). We used 36 points in these modeling analyses (Fig. S1). For estimation of current potential distribution, 19 bioclimatic variables in the WorldClim database (http://www.worldclim.org/download) were applied. To simulate the distribution during Last Glacial Maximum, we followed the Model for Interdisciplinary Research on Climate (Hasumi & Emori, 2004) and the Community Climate System Model (Collins et al., 2006). For the model and the training data, we used the AUC value to estimate the goodness of fit.

Sequence differences and haplotype patterns of cpDNA

After alignment, the rpoB-trnC and trnH-psbA sequences were 1,145 bp and 367 bp respectively, and 1,512 bp in total. In the combined data, we identified seven haplotypes (H1–H7) (Table 1) and 17 polymorphic sites (14 substitutions and three indels) from 269 individuals collected from 27 populations (Table 2). The total haplotype diversity (H_cp) was 0.6691 with a within-population variation of 0–0.5333 while the overall nucleotide diversity (πcp) was 0.0077 with a range of 0–0.0130. The BYT2 population had higher H and π values than others (Table 1).

Among the 27 populations examined, only three of them had two haplotypes and all others a single haplotype. In the seven haplotypes identified among P. lanata samples, H1 was widely distributed and was the only haplotype in most populations. H5 was found in four populations, three of which were “H5-only” ones. H2, H3, and H7 were exclusively detected in three different populations. H4 and H6 were specific haplotypes found in BYHT and BYT2, respectively (Table 1; Fig. 1).

Sequence differences and haplotype patterns of ITS

The attempts for amplification and sequencing of ribosomal DNA from several individuals in each of the 27 populations failed. Two hundred individuals were successfully sequenced. The total length was 635 bp, and thirty-two haplotypes (C1-C32) were obtained (Table 1). The total haplotype diversity (H_nr) was 0.5668, with a range of 0.0769–0.8918, whereas the nucleotide diversity (πnr) was 0.0018, with a range of 0.0001–0.0059. Populations JLT, SZ, BYT1, BYHT, XN, EGB, BYT2, XJW, ELT and DB had high H and π values. Most populations had Haplotype C1 or C2.

Genetic diversity and population structure

SAMOVA analysis divided all populations of P. lanata into three groups based on the cpDNA analysis: (1) the northern group (populations 1–3), (2) the central group (populations 4–11), and (3) the eastern group (populations 12–27).

Among all of the populations, the genetic diversity was high (cpDNA: H cp = 0.6691 ± 0.0259, H_T = 0.673 ±0.0869; nrDNA: H nr = 0.5668 ± 0.0277, H_T = 0.577 ± 0.0685), while the average within-population genetic diversity was relatively low in the cpDNA (H_S = 0.033 ±  0.0214) (Table 3). The cpDNA showed high levels of inter-population genetic diversity (G_ST = 0.950 ± 0.0308, N_ST = 0.922 ± 0.0559), suggesting that genetic differentiation mainly occurred among populations. However, as indicated by the nrDNA sequences, the genetic differentiation among the populations (G_ST = 0.348 ±  0.0562) was low. In either cpDNA or nrDNA, N_ST and G_ST showed no significance regarding their difference (p > 0.05), thus no clear phylogeographic structure is present in P. lanata. In the central group, the genetic diversity was high (H cp = 0.6673 ±  0.0177, H_T = 0.753 ± 0.0476), with a low average intra-population genetic diversity (H_S = 0.021 ± 0.0208). The level of genetic differentiation among populations was high (G_ST = 0.972 ± 0.0263, N_ST = 0.755 ±  0.1630). By contrast, both the genetic diversity (H cp = 0.1364 ± 0.0373, H_T = 0.137 ±  0.1122) and the average within-population genetic diversity (H_S = 0.046 ± 0.0348) of the eastern group were low (Table 3). Of the total variation, 89.87% (p < 0.001) occurred in the cpDNA sequences among northern, central, and eastern groups, 5.22% was inter-populations variation within groups, and 4.91% was intra-populations variation (Table 4). This result of the hierarchical analysis was consistent with the F_ST (0.9509, p < 0.001) and F_CT (0.8987, p < 0.001) values. However, genetic variation in the nrDNA mainly occurs within populations (F_Sc = 0.30026; F_CT = 0.08700) (Table 4).

Population divergence and demography: cpDNA analysis

The parameters of Tajima’s D and Fu’s F_S (p > 0.05) were positive and insignificant, indicating that P. lanata has not experienced a recent expansion (Table 5; Fig. 2A). In the central group, the Tajima’s D and Fu’s F_S (p > 0.05) were insignificant and positive, indicating that the group has not experienced recent expansion. This conclusion was also confirmed by the P-values (p < 0.05) of SSD, Hrag, and the multimodal mismatch analysis (Table 5; Fig. 2B). In the eastern group, the parameter of Tajima’s D (−1.8814, p < 0.05 ) was significantly negative, with the P-values of SSD and Hrag were both greater than 0.05. In addition, the mismatch distribution was unimodal, indicating that the populations had experienced a regional-scale expansion (Table 5; Fig. 2C). In the northern group, H5 was the only haplotype, and therefore it was impossible to analyze its population expansion further.

Meanwhile, Beast analysis unraveled that the occurrence time of divergence between these seven haplotypes was from the early Pleistocene (1.6053; 95% HPD: 0.7215–2.8307) Mya to the late Pleistocene (0.0857; 95% HPD: 0.0033–0.2721) Mya (Fig. 3).

Species expansion trends: current and future distribution

The AUC values of P. lanata were 0.997/0.997 (the current model) and 0.997/0.997 (under the MIROC climate model). The higher AUC values indicated that the model was more suitable for the current distribution and the potential distribution during the LGM. The results revealed that the current distribution of P. lanata has not experienced recent expansion. The central region remained stable, whereas the eastern region showed an expanded distribution pattern (Figs. 4A, 4B).

Genetic diversity of P. lanata

We found the high genetic diversity in P. lanata. The value of genetic differentiation (G_ST = 0.950) was higher than the average reported for other angiosperms (G_ST = 0.637) (Petit et al., 2005), indicating a stronger inter-population differentiation of the cpDNA sequences in P. lanata. The total genetic diversity of P. lanata (H_T = 0.673) was similar to that of Allium mongolicum (H_T = 0.693) but lower than that of Lagochilus ilicifolius (H_T = 0.925), both of which are herbaceous plants in northwest China. The average within-population genetic diversity (H_S = 0.033) of P. lanata was lower than Allium mongolicum (H_S = 0.180) (Meng & Zhang, 2011; Zhang et al., 2017), mainly due to its morphology and the influence of gravity on seed dispersal. Since the seeds can only be spread over short distances, a severely restricted gene flow between populations was resulted (Meng & Zhang, 2011). Concerning the nuclear gene sequence, the genetic diversity (H_T = 0.577, H_S = 0.376) was higher due to the long distance of pollen transmission. High genetic variation and unique haplotypes are usually associated with centers of plant diversity or potential refugia, whereas regions of recent colonization have low levels of genetic variation (Li et al., 2010; Meng & Zhang, 2011; Stewart et al., 2010). The central group had the highest genetic diversity. The Helan Mountains is a diversity center for species, and a high level of genetic diversity is expected. Such diversity centers get support from other species distributed in the area (Meng & Zhang, 2011; Shi & Zhang, 2015). In the eastern group, haplotype H1 was widely distributed in most populations and the most haplotypes for all populations except HQH and BYT2. During the range expansion process, the founder effect should be responsible for less genetic diversity, often causing a single prevailing haplotype (Hewitt, 2000; Xu & Zhang, 2015b; Zhang, Volis & Sun, 2010). Therefore, the low genetic diversity in the eastern group should be attributed to the founder effect.

Intraspecific divergence across DNA sequences

All populations of P. lanata were divided into three groups. In the eastern group, the haplotype H1 was dominant, and H2, H5, and H6 were rare haplotypes. The central group was dominated by the haplotypes H2, H3, and H7, and H4 was rare. The northern group only had the haplotype H5. The haplotypes were unique among the three groups except for H2 and H5, suggesting a division between the three groups. AMOVA analyses showed that the variation among the three groups contributed to the total variation, indicating the presence of inter-group restricted gene flow. The occurrence of divergence among P. lanata populations was in the early to late Pleistocene, as the aridification and desert expansion also occurred during Pleistocene in northwest China. Therefore, we speculate that this species showed a diversification distribution, which may be attributed to the aridification in the early Pleistocene. In addition, the northern and eastern groups are separated by the desert (Hobq Desert, Mu Us Sandy Land, and Ulan Buh Desert), the geographical barriers that may limit the gene flow among the three groups, causing genetic differentiation within the species. The central group is located around the Helan Mountains. Previous studies have shown that the Helan Mountains acted as migration corridors for recolonization after ice ages as well as a geographical barrier (Meng et al., 2015; Zhang et al., 2017). Therefore, we speculate that deserts and the Helan Mountains may act as geographical barriers restricting the long-distance dispersal of P. lanata seeds, eventually resulting in isolation and differentiation of the populations in these three groups (northern, eastern, and central groups). Nuclear gene data showed that low genetic differentiation among the three groups. Because pollen transmission was far away, geographical barriers (deserts) did not hinder the gene flow among the three groups.

Demographic history of P. lanata

According to the neutrality test and mismatch analysis, no expansion event has occurred in this species recently, except the eastern group. Previous studies have shown that plants have lower genetic diversity and a single, widely distributed haplotype in regions with rapid expansion (Hewitt, 1996; Shi & Zhang, 2015; Zhang & Zhang, 2012). The eastern group has a low number of haplotypes and limited genetic diversity. H1 is widely distributed and is the only haplotype for most populations in the eastern group, indicating that the region has experienced an expansion. The mismatch analysis and neutrality test also confirmed these results. In comparison to the LGM distribution of P. lanata, the central region showed a stable distribution pattern, whereas the eastern region showed an expanding distribution. This may be due to the high annual precipitation in the central region and eastern region allowing for more suitable habitat for P. lanata.