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

Climatic fluctuations during the Quaternary significantly affect many species in their intraspecific divergence and population structure across northwest China. In order to investigate the impact of climate change on herbaceous plants, we studied Panzerina lanata (Lamiaceae), a widely distributed species. Sequences of two chloroplast DNA (cpDNA) intergenic spacers (trnH-psbA and rpoB-trnC) and a nuclear ribosomal region (nrDNA, ITS) were generated from 27 populations of Panzerina lanata and resulted in the identification of seven chloroplast haplotypes and thirty-two nuclear haplotypes. We applied AMOVA, neutrality test and mismatch distribution analysis to estimate genetic differentiation and demographic characteristics. The divergence times of the seven cpDNA haplotypes were estimated using BEAST. Our results revealed high levels of genetic diversity (cpDNA: Hcp = 0.6691, HT = 0.673; nrDNA: Hnr = 0.5668, HT = 0.577). High level of genetic differentiation (GST = 0.950) among populations was observed in the cpDNA sequences, while the genetic differentiation values (GST = 0.348) were low in nuclear sequences. AMOVA results revealed major genetic variation among the three groups: northern, central, and eastern group. However, the genetic differentiation in ITS data was not found. The species distribution modeling and demographic analysis indicated that P. lanata had not experienced recent range expansion. The occurrence of divergence between seven cpDNA haplotypes, probably during Pleistocene, coincides with aridification and expansion of the desert across northwest China that resulted in species diversification and habitat fragmentation. In addition, we discovered that the deserts and the Helan Mountains acted as effective geographic barriers that promoting the intraspecific diversity of P. lanata.


INTRODUCTION
In China, studies on plant phylogeography have been mainly focusing on four regions: Qinghai-Tibet Plateau and Southwest China, West China, North and Northeast China, South and Southeast China (Liu et al., 2012;Qiu, Fu & Comes, 2011). The Qinghai-Tibet Plateau and southwest China are among the hotspots of biodiversity research in the world, which were not directly affected by ice shield. Most studies have shown that phylogeographic patterns were mainly influenced by climate fluctuations during the Quaternary, resulting in intraspecific divergences and regional range expansion (Guo et al., 2010;Liu et al., 2012). However, recent research on plant phylogeography in northwest China has expanded (Meng et al., 2015;Wang et al., 2016;Zhang et al., 2017), as its arid area includes not only Xinjiang but also Hexi Corridor, Qaidam Basin and western Helan Mountains (Dang & Pan, 2001). In these areas, the genetic structure and phylogeography of plants were mainly influenced by Quaternary climatic fluctuations (Meng et al., 2015;Wang et al., 2013). Climate change during the Pleistocene led to the low temperatures and aridification in northwest China, which promoted the desert expansion. Aridification and desert expansion have significant impacts on the phylogeography of many species in northwest China (Meng et al., 2015;Wang et al., 2016), such as allopatric divergences, speciation and habitat fragmentation of desert plants (Ma, Zhang & Sanderson, 2012;Meng & Zhang, 2011;Wang et al., 2016;Xu & Zhang, 2015a). In addition, geographical barriers (deserts, mountains, etc.) have resulted in segregation of species, limited seeds dispersal, and therefore rare genetic communication between populations, which lead to isolation and differentiation of genetic lineages of species (Cun & Wang, 2010;Liu et al., 2012;Wang et al., 2016). Research in these regions mainly concentrated on shrubs rather than herbaceous plants (Shi & Zhang, 2015;Su, Lu & Zhang, 2016;Su et al., 2015;Xu & Zhang, 2015a) that are more sensitive to climate oscillation. Here, we selected Panzerina lanata (Lamiaceae) as a suitable model to study the genetic structure of desert species in arid northwest China and its response to Quaternary climatic fluctuations.
The genus Panzerina (Lamioideae, Lamieae) contains two species (Panzerina canescens and Panzerina lanata ) that are mainly distributed in the desert and desert grassland areas of central Asia, and have been described in the Flora of China (Li & Hedge, 1994). Panzerina lanata is a perennial medicinal herb and widely distributed in the sandy desert steppes of Inner Mongolia, Gansu, Ningxia, and Shanxi. By contrast, Panzerina canescens, with a small population, is only found in the dry stony area of Xinjiang. Previous studies on P. lanata have mainly focused on biological characteristics (Cheryomushkina & Astashenkov, 2014;Li & Zhao, 2000), chromosome research (Li, Cao & Zhao, 1999), plant taxonomy and floristic analyses (Zhao et al., 1998;Zhao & Liu, 1997), other than its genetic diversity and phylogeography pattern.
Here, the phylogeographic structure of P. lanata was inferred by investigating two chloroplast DNA (cpDNA) intergenic spacers (trnH-psbA and rpoB-trnC) and a nuclear ribosomal region (nrDNA, ITS). Our objectives are to determine (1) the genetic diversity and structure of P. lanata and (2) how climatic fluctuations and geographical obstruction contribute to lineage differentiation of P. lanata.

Population sampling
Across northwest China, there are approximately 27 natural populations of P. lanata, from which we collected 269 individuals (Table 1) 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).

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: intergroups, 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.

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 Haplotype , TATTTAGATATGTT; , TTTTCT; •, TATTTATCTATTCTATTTTCA. 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.  (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

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

DISCUSSION 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 withinpopulation 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.

CONCLUSIONS
Our results show that the P. lanata recently has not experienced range expansion. Both chloroplast data and nuclear gene data showed high genetic diversity. Because of the different pollination mechanisms, the gene flow of cpDNA is mainly mediated by seeds while the nuclear DNA by both seeds and pollen. The chloroplast data indicated that the variation occurred mainly among three groups. Aridification and geographical isolation (Deserts and Helan Mountains) limited gene flow among populations and played critical roles in affecting genetic diversity and genetic differentiation of P. lanata. According to our nuclear gene data, the geographical isolation such as desert and Mountains had less influence on the differentiation of P. lanata due to the ability of pollens for longer distance traveling.