Phylogeography of Schisandra chinensis (Magnoliaceae) Reveal Multiple Refugia With Ample Gene Flow in Northeast China

Temperate conifers and broadleaved mixed forests in northeast China are ideal to investigate the genetic consequences of climate changes during the last glacial maximum (LGM), 29 – 16 kya. As previous studies were focused on tree species with long generation time; here, the evolutionary history of Schisandra chinensis, a climber species with a generation time of five years, was investigated using chloroplast DNA (cpDNA), nuclear single copy gene (nSCG), and nuclear single sequence repeats (nSSRs, i.e., microsatellite) markers, along with ecological niche modeling (ENM), which predicted a suitable habitat in Korea Peninsula (KP) during the LGM. Private haplotypes and high genetic diversity of both cpDNA and nSCG were mainly found in KP and Changbai Mt. (CB). Although no significant phylogeographic structure was detected in the cpDNA and nSCG, three nSSRs clusters roughly distributed in west (CB and KP), east (north China), and north (Xiaoxing’an Range, XR) regions were found in Structure analysis. The approximate Bayesian computation analysis showed the west cluster diverged at 35.45 kya, and the other two clusters at 19.85 kya. The genetic diversity calculated for each of the three markers showed no significant correlation with latitude. Genetic differentiation of nSSRs was also not correlated with geographic distance. Migrate analysis estimated extensive gene flow between almost all genetic cluster pairs and BOTTLENECK analysis showed that few populations experienced severe bottlenecks. Overall, results indicate that S. chinensis survived the LGM in situ in multiple refugia, which likely include two macrorefugia (KP and CB) and two microrefugia (XR and north China). Extensive postglacial gene flow among the three nSSRs clusters led to uniformly distributed genetic diversity and low genetic differentiation.


INTRODUCTION
Quaternary climate changes, especially the last glacial maximum (LGM, 29-16 kya; Clark and Mccabe, 2009), have greatly influenced the distribution of temperate forests in the Northern Hemisphere (Hewitt, 1996(Hewitt, , 2000Qiu et al., 2011). Contractions during glaciations and expansions between or after glaciation periods of species distributions have left imprints on the genetic structure and distribution of genetic diversity and differentiation (Provan and Bennett, 2008;Excoffier et al., 2009;Waters et al., 2013). In northeast China, the temperate conifer and broadleaved mixed forests (hereafter mixed forests) between 40 and 50 • N are the contact zone between the southern warm-temperate forests and the northern cool-temperate forests (Supplementary Figure S1; Wu, 1980). Because the mixed forests are sensitive to climate cooling and warming, it is an ideal system to investigate the genetic consequences of the LGM (Wu, 1980;Zhang, 2007).
Vegetation reconstructions using fossil pollen data showed that the mixed forests retreated southward to 25-30 • N during the LGM (Harrison et al., 2001;Cao et al., 2015). However, phylogeographic studies provided different scenarios. Two scenarios have been proposed, which are the single refugium scenario and the multiple refugia scenario (Supplementary Figure S1). The single refugium scenario suggests that Changbai Mt. (CB), located at about 40-42 • N, is the only macrorefugia for conifer trees within the genus Abies (Jiang et al., 2011), the broadleaf tree Juglans mandshurica (Bai et al., 2010), and the perennial herb Buplerum longiradiatum (Zhao et al., 2013). The multiple refugia scenario suggests the Korea Peninsula (KP) as another macrorefugia in addition to CB (Guo et al., 2014;Wang et al., 2014;Bao et al., 2015;Zeng et al., 2015). Other microrefugia were also detected, with the northernmost ones located at the northern mixed forests margin [Xiaoxing'an Range (XR) and Russian Far East] (Bao et al., 2015;Zeng et al., 2015;Wang S.H. et al., 2016). Most species examined were tree and shrub species. The only climber species that has been studied is Actinidia arguta, but its limited genetic variation in chloroplast DNA (cpDNA) hindered the precise discrimination of the two different scenarios (Ye et al., 2018). Thus, whether climber species would follow the single refugium or multiple refugia scenario remains unknown.
Genetic diversity and genetic differentiation distribution patterns in mixed forests are also complex (Ye et al., 2017). Genetic diversity is expected to show significant declines due to genetic drift and bottlenecks during the northward expansion from a single source (such as the CB macrorefugia) (Excoffier et al., 2009;Waters et al., 2013). The reduction of genetic diversity in nuclear microsatellites (nuclear single sequence repeats, nSSRs) of Acer mono with increasing latitude agrees with this prediction (Guo et al., 2014;Liu et al., 2014). However, J. mandshurica shows an inconsistent pattern. Gradual expansion with gravity-mediated seed dispersal combined with extensive wind-mediated pollen gene flow have prevented the decrease of its genetic diversity in the direction of range expansions (Wang W.T. et al., 2016). In addition, if there was substantial postglacial southward expansion from northern microrefugia (such as XR), the genetic diversity would also become uniformly distributed (Bao et al., 2015;Zeng et al., 2015). Genetic differentiation can also show different patterns. An increased genetic differentiation with increased distance (i.e., isolation by distance, IBD) pattern in re-colonized areas is predicted under rapid range expansion, while ample gene flow among source and colonized populations would erase the IBD pattern (Lafontaine et al., 2013). It should be noticed that the above-mentioned studies were performed on tree species with long generation time. Thus, the distributions of genetic diversity and genetic differentiation in plants with short generation time need additional investigation.
Schisandra chinensis (Turcz.) Baill, a common climber species with short generation time (about five years) (Li and Zhang, 2017), was chosen for the present study. Its seeds are dispersed by birds (Yuan, 2007), while its pollen is mainly dispersed by insects and, occasionally, by wind (Ai et al., 2007). The genetic patterns in cpDNA, nuclear single copy gene (nSCG), and nSSRs were combined with ecological niche modeling (ENM) to investigate the evolutionary history of S. chinensis. Firstly, genetic structure and potential habitat predictions were used to infer potential refugia. Potential divergence histories were compared and parameters of the most possible scenario were estimated using approximate Bayesian computation (ABC). Then, the correlations between genetic diversity and latitude or IBD were calculated. Historical gene flow was estimated using the maximum-likelihood (ML) algorithm, and whether populations have experienced genetic bottleneck was also estimated.
To further understand the influence of the LGM on mixed forests, the present study aimed to (1) determine whether S. chinensis conforms to the single refugium or multiple refugia scenario, and (2) reveal the distribution patterns of genetic diversity and genetic differentiation of S. chinensis and their potential causes.

Sampling
We collected 355 individuals of S. chinensis from 20 populations ( Table 1) and four individuals of Kadsura longipedunculata were collected as outgroups. There were at least 30 m between any two individuals in each population. Silica gels was used to desiccate and preserve all leaf samples. Voucher specimens were deposited at the Beijing Normal University Herbaria (BNU), Beijing, China.

DNA Extraction and Sequencing
Total genomic DNA was extracted using a Plant Genomic DNA Kit (DP305-03; Tiangen, Beijing, China). We amplified and sequenced four chloroplast gene fragments, namely matK, ndhA, trnL-trnF, and trnS-trnG (Supplementary Table S1) for 148 individuals and the nuclear gene PEPC for 167 individuals ( Table 1). The PCR amplification was performed as described by Guo et al. (2014). The amplicons were sequenced from both directions at the Beijing Genomics Institute (Beijing, China). All electropherograms of DNA sequences were visually analyzed before they were assembled and read in CodonCode Aligner 3.6.1 (CodonCode Corporation 1 , Centerville, MA, United States).

Genetic Diversity and Construction of the Most Parsimonious Network
CodonCode Aligner 3.6.1 with the CLUSTAL module was used to align all sequences. The alignment was verified visually and low-quality sections at the beginning and end of sequences were deleted. The PHASE function in DnaSP 5.10.01 (Rozas et al., 2003) was used to determine heterozygous sequences in the nSCG. Then, DnaSP 5.10.01 was used to determine haplotypes, without considering indels, and to calculate nucleotide diversity (π) and haplotype richness (HR). Permut 1.0 2 with 1000 permutations was used to compare two forms of genetic differentiation (G ST and N ST ) in all populations. The most parsimonious network was inferred using Network 4.6.1.1 (Bandelt et al., 1999). Pearson correlations ("corr.test" function) between HR and latitude were performed in R 3.2.3 (R Core Team, 2013).

Amplification and Genotyping
Eight nSSRs loci were used to determine the genotypes of all sampled individuals in each population (Supplementary Table S2). The PCR procedures were the same as those for amplifying nuclear and cpDNA sequences, except for the annealing temperature (Supplementary Table S2). Amplicons were loaded onto a 3730XL Automated Genetic Analyzer (Applied Biosystems, Foster City, CA, United States) and scored using GeneMapper 4.0 (Applied Biosystems). To reduce score error, alleles were independently read by two people, and any disputes (<5% of calls) were decided by a third person.

Genetic Diversity and Bottleneck
The neutral test in Lositan (Antao et al., 2008) indicated that no loci violated the neutral hypothesis; thus, all eight loci were used for further analyses. For each locus and population, standard genetic diversity statistics were calculated. The deviation of the fixation index (F IS ) from zero was used to test the deviation from Hardy-Weinberg equilibrium in each population. Genotypic disequilibrium was tested for all loci pairs in all populations by randomization, and the obtained P-values (=0.05) were adjusted using the Bonferroni correction. All the above calculations were performed in Fstat 2.9.3 (Goudet, 2001). Allele richness (R S ) and private allele richness (P AR ) in all populations were calculated in hp-rare 1.0 (Kalinowski, 2005) using rarefaction with a sample number of 10. Pearson correlations ("corr.test" function) between nSSRs genetic diversity (R S , P AR , and H e , expected heterozygosity) and latitude were performed in R 3.2.3. Severe population size decreases were detected in BOTTLENECK (Cornuet and Luikart, 1996), applying the two-phase mutation (TPM) with 70% stepwise mutation model (SMM) and 30% multistep mutations.

Population Structure
The genetic differentiation index F ST (Weir and Cockerham, 1984) was determined for the eight loci using Fstat 2.9.3. Genetic distance, Nei's Da, was used to create a neighbor-joining (NJ) tree using Poptree2 (Takezaki et al., 2010) with 1000 bootstrap pseudoreplicates. Structure 2.3.4 (Falush et al., 2007) was used to detect potential population structure without using population location as prior information. Ten independent runs were performed for each K (K = 1-20) with 1 × 10 5 initial burnin, followed by 1 × 10 6 Markov Chain Monte Carlo (MCMC) steps using an admixture model with correlated allele frequencies.
Potential clusters (K) were determined by LnP(D), the change in log-likelihood of the data for each run (Pritchard et al., 2000), and by K, the second-order rate of change of LnP(D) between successive K values (Evanno et al., 2005). Isolation by distance (Wright, 1943) was evaluated according to Rousset (1996) using a Mantel test between genetic differentiation in terms of F ST /(1 -F ST ) and the natural logarithm of geographic distance. The Mantel test was performed with 9999 permutations in GenAlEx 6.5 (Peakall and Smouse, 2012).

Population Divergence History and Gene Flow
To detect population divergence history, all populations were assigned to one of the three gene pools corresponding to its largest proportion of ancestry; population LW ( Table 1) was the exception as its proportion of ancestry was approximately identical for the three clusters (see section "Results"). Population divergence history was modeled using the ABC procedure (Beaumont et al., 2002) and performed in DIYABC 2.0 (Cornuet et al., 2014). Three possible divergence history scenarios were compared among the three genetic clusters (see section "Results"; Supplementary Figure S2 and Supplementary Table S3). The simulations were summarized based on the following statistics: mean number of alleles, mean genetic diversity and F ST for each lineage, mean classification index, and shared allele distance between pairs of lineages. The simulations were repeated 3,000,000 times and, after logit transformation, local linear regression was applied to choose the 1% simulated data sets that were closest to the observation. Using these data sets, we compared the posterior probability (PP) of the three scenarios and estimated the parameters of the most possible scenario. The historical gene flow among the nSSRs clusters was estimated using the ML algorithm in Migrate 3.6.8 (Beerli, 2006). We estimated θ = 4Nµ (N, the effective population size; µ, the mutation rate per locus per generation) and M = m/µ (m, the migration rate per generation) to calculate the effective number of migrants (Nm). We used 10 short chains (10,000 trees), three long chains (100,000 trees), 20 genealogies, and burn-in of 10,000 initial trees. The starting values for θ and M were estimated from F ST values. All individuals were included and the analyses were run three times, independently.

Ecological Niche Modeling
The maximum entropy modeling technique (MAXENT) (Phillips et al., 2006) was utilized to predict the current potential distribution of S. chinensis as well as those during the LGM and last interglacial (LIG). Fifty-five presence records were obtained from this study (20 sample sites) and from the Global Biodiversity Information Facility 3 (35 occurrence locations). Seven climate variables with low correlation (<0.8) were used to model the niche (Supplementary Table S4) (Hijmans et al., 2005). Model validation was performed 10 times independently, with 20% randomly chosen data. The accuracy of model predictions was evaluated based on the area under the receiver operating characteristic (ROC) curve (AUC) (Fawcett, 2006).
The established model was projected to both MIROC 3.2 4 and CCSM4 (Shields et al., 2012) models to predict potential distributions during LGM and LIG distribution. All climate layers were prepared using a 2.5 arc-minute resolution. The paleocoastlines during the LGM were estimated assuming that the sea level was 130 m lower than the current sea level.

cpDNA and nSCG Sequences
In cpDNA, seven variable sites, including two singleton variable sites and five parsimony informative sites, were detected (Supplementary Table S5). Six haplotypes were identified with a total alignment length of 3,507 bp (Supplementary Table S5). Haplotype 2 (H2), located in the center of the network, was likely the ancestral haplotype. It was also the most abundant haplotype, being distributed in almost all populations. Only three populations had private 4 http://pmip2.lsce.ipsl.fr/ haplotypes: FY (H4), in northernmost Changbai Mt., and ZY (H6) and JW (H5) in Korea Peninsula. These three populations also had high HR (Table 1); other populations with high HR were mostly distributed in Changbai Mt. (DS, BS, and LW) and in the northern mixed forests margin (KD and RH) (Figure 1).
In the nSCG, twelve variable sites, including three singleton variable sites and nine parsimony informative sites, were detected. Twelve haplotypes were identified with a total alignment length of 809 bp (Supplementary Table S6). Because nuclear haplotype 6 (N6) was located in the center of the network, it was likely the ancestral haplotype. It was distributed in north China (LM, QS, and CB) and Changbai Mt. (CB). Populations ZY (N11 and N12), JW (N8 to N10), and CH (N7) in northernmost Changbai Mt. had private haplotypes. Private haplotypes in Korea Peninsula (N8, N9, and N10 in JW, N11 and N12 in ZY) populations were not closely linked (Figure 2). All populations except XR (HR = 0.250) had high HR ranging from 0.556 to 0.867 (Table 1)  were found. All the newly obtained sequences were uploaded to GenBank (Accessions MH580160-MH580199).

Nuclear Microsatellites
The genetic diversity of the eight nSSRs loci and 20 populations are shown in Table 1 and Supplementary Table S7. For each population, F IS showed no significant deviation from zero ( Table 1), suggesting no violations of the Hardy-Weinberg equilibrium assumptions. No significant genotypic disequilibrium was observed among the 28 loci pairs in any population. No significant correlations between genetic diversity and latitude were found when all populations were considered (H E , r = 0.35, P = 0.13; R S , r = 0.20, P = 0.40; P AR , r = -0.30, P = 0.19) or when populations with sample size < 10 were excluded (H E , r = -0.30, P = 0.32; R S , r = -0.45, P = 0.12; P AR , r = -0.23, P = 0.45) (Figure 3). An IBD pattern was not detected (r = 0.22, P = 0.06; Figure 4). Two populations (KY and MJ) with limited sample size and population HN experienced bottlenecks ( Table 1).
In Structure analysis, K was high for K = 2 or K = 3 (Supplementary Figure S3), indicating that populations can be clustered into two or three groups. However, LnP(D) and K were both higher at K = 3 than at K = 2, indicating that three clusters were more likely than two clusters. The three clusters were roughly distributed in the east, west and north regions ( Figure 5). No population structure was detected using the NJ tree (Supplementary Figure S4).

Ecological Niche Modeling
High ROC values (0.958 ± 0.012) indicated good accuracy of model predictions. During LGM, both models predicted range FIGURE 3 | Correlations between genetic diversity, H E , expected heterozygosity (A); R S , allelic richness (B); and P AR , private allelic richness (C) in eight nuclear microsatellites, with latitude in Schisandra chinensis. Light gray dots represent population with sample size lower than ten.
contractions. The CCSM4 predicted that the most suitable habitat was located at KP, while MIROC 3.2 indicated the Yellow Sea land bridge and adjacent regions and KP as the most suitable locations. During the LIG, suitable habitats were predicted to be more limited than that at present (Figure 6).

In situ Glacial Survival in Multiple Refugia
Previous phylogeographic studies in mixed forests provide limited information for climber species (Ye et al., 2017). Here, we investigated the evolutionary history of S. chinensis, a common climber species in mixed forests, in detail based on multiple genetic markers and potential habitat predictions. Contrary to the southern contraction of mixed forests during the LGM inferred by vegetation reconstructions (Harrison et al., 2001;Cao et al., 2015), the current investigation clearly showed that S. chinensis survived the unsuitable climate during the LGM in situ in multiple refugia, with the northernmost contacting the northern mixed forests margin.
Low genetic differentiation (G ST = 0.311 in cpDNA, G ST = 0.065 in nSCG, and F ST = 0.130 in nSSRs) was revealed by the different genetic markers and no phylogeographic structure was revealed in either cpDNA and nSCG and in the NJ tree of nSSRs. However, nSSRs Structure analysis indicated that populations could be divided into three clusters (west, north, and east) roughly corresponding to their geographic ranges. The DIYABC analysis showed the west cluster diverged before the other two clusters (35.45 kya) and before the LGM, and persisted through the LGM. The divergence time between the north and east clusters (19.85 kya) indicated their divergence was likely triggered by LGM and these two clusters have persisted in situ since the LGM. Thus, the three S. chinensis clusters probably had multiple LGM refugia.
The macrorefugia of S. chinensis can be inferred through potential habitat prediction. The Korea Peninsula was suggested as one of the macrorefugia for this species, as it provided a suitable habitat during the LGM in both climate models in the   N 1 ( × 10 4 ) N 2 ( × 10 4 ) N 3 ( × 10 5 ) t 2 ( × 10 3 ) t 1 ( × 10 3 ) µ ( × 10 −5 ) P N 1 , N 2 , N 3 : The current population size of east cluster, north cluster, west cluster and ancestral populations at t 2 , respectively; t: divergence time; µ: mutation rate (per generation per locus); P, the proportion of multiple step mutations in the generalized stepwise model, GSM.  ENM. Further, KP populations showed high proportion of private alleles in all three markers indicating long-term persistence. Although the ENM failed to detect a potential suitable habitat in CB during the LGM, populations in CB showed high proportion of private alleles and high genetic diversity in cpDNA and nSCG data. Based on these data and on previous phylogeographic studies (Ye et al., 2017), CB might be another macrorefugia. The Changbai Mountains harbor several forest types, and glacial advances in late Pleistocene only took place at elevations above 2000 m above sea level (Zhang et al., 2008). Thus, low elevation regions with varied vegetation types could serve as a refuge for S. chinensis. This hypothesis also explains the existence of the west cluster, which is mainly distributed in KP and CB, and has persisted for longer and much higher effective population size than other two clusters as inferred by the DIYABC analysis. Potential microrefugia cannot be effectively detected by ENM, because it uses spatially and temporally smoothed climate data, thus being unable to capture climatic variance and the effects of topography on microclimate (Gavin et al., 2014). As proposed by Zeng et al. (2015), our nSSRs also revealed the existence of northern microrefugia. The north nSSRs cluster was mostly distributed in the XR and it diverged during the LGM (19.85 kya) indicating that XR, at the northern mixed forests margin, might be a microrefugium. Moreover, the nSSRs cluster distributed in the east indicated another possible refugium in north China. The ancestral nSCG haplotype N6 distributed mainly in north China also supports the existence of a refugium in this region. Thus, contrasting to the mixed forests' retreat to 30 • N proposed by vegetation reconstructions using fossil pollen data (Harrison et al., 2001;Cao et al., 2015), S. chinensis is proposed to have survived the LGM in situ in multiple refugia, i.e., in KP, CB, XR, and north China.
Our study on S. chinensis provides information on the evolution of another climber species in mixed forests. Ye et al. (2018) suggested repeated expansions and fragmentations in A. arguta linked to Pleistocene climate changes have shaped its genetic structure across its distribution from northeast China to subtropical China, although limited cpDNA variation was found. The present study on S. chinensis using multiple datasets provides a more detailed evolutionary history of a climber species and a substantial complement to previous mixed forests studies. A brief evolutionary history of mixed forests can be thus be inferred. The mixed forests have survived in situ during the LGM using KP and CB as the most important macrorefugia (Jiang et al., 2011;Guo et al., 2014;Wang et al., 2014;Zeng et al., 2015) and microrefugia located in other southern or northern regions were also used. Southern microrefugia can be found in north China (Wang S.H. et al., 2016;Wang W.T. et al., 2016), Shandong Province (Ye et al., 2018), and northern microrefugia can be found in the northern mixed forests margin, such as XR (Bao et al., 2015;Wang S.H. et al., 2016) or the Russian Far East (Zeng et al., 2015).

Uniformly Distributed Genetic Diversity and Genetic Differentiation
Assuming that the two southern macrorefugia (KP and CB) as sources for northward range expansion, S. chinensis is certainly expected to show significant genetic changes along latitude (Excoffier et al., 2009;Waters et al., 2013). However, no significant decreases in genetic diversity, allelic richness, or allelic privacy in organelle and nuclear markers were detected with latitude increase. The same genetic diversity distribution pattern was revealed for Pinus koraiensis (Bao et al., 2015). The ABC procedure estimated significant gene flow between southern macrorefugia (KP or CB) and the northern microrefugium XR in P. koraiensis. Bao et al. (2015) suggested that this ample gene flow resulted in uniformly distributed genetic diversity and shared dominant haplotypes. In Quercus mongolica, Zeng et al. (2015) suggested that if northern microrefugia (such as the Russian Far East) contributed little to post-glacial expansions, the genetic diversity of this species would also show a significant decline. In agreement with that found for P. koraiensis (Bao et al., 2015), the present study revealed a northern macrorefugium for S. chinensis at XR, and Migrate estimations showed extensive gene flow among the three Structure clusters, the greatest being found from the north to the west cluster (4Nm north→west = 113.68, Table 3). Thus, after the in situ LGM persistence of the different S. chinensis nSSRs clusters, postglacial seed movement by birds (Yuan, 2007), pollen movement by wind (Ai et al., 2007), and the species short generation time led to ample gene flow among the different clusters. Most of the genetic diversity present in macrorefugia (CB and KP) or in the microrefugium (XR) has been preserved in populations in the contact zone among different refugia, and private cpDNA or nSCG haplotypes can only be found in a few S. chinensis populations. Therefore, no significant negative correlations with latitude were found for genetic diversity indices using the three independent genetic markers. The homogenizing gene flow among all populations also resulted in the absence of IBD.

CONCLUSION
A detailed evolutionary history was revealed for S. chinensis using multiple genetic markers and potential habitat predictions. Divergence times among the three nSSRs clusters showed that S. chinensis persisted in situ through and since the LGM. Two macrorefugia, at KP and CB, were revealed by ENM and high allelic privacy and genetic diversity distributions. Haplotype and nSSRs cluster distributions indicated two possible microrefugia located at XR and north China. Extensive gene flow among clusters resulted in the uniform distribution of genetic diversity and in the absence of an IBD pattern. Short generation time, and seed and pollen movement have facilitated this extensive gene flow. The present study using a common climber species with short generation time is an important complement of previous studies that mainly used tree species with long generation times.

ACKNOWLEDGMENTS
We are grateful to Da-Yong Zhang, Shou-Hsien Li, and Victoria L. Sork for useful discussions and insightful comments and to Sheng-Hong Wang, Xi-Di Guo, and Tao Jiang who assisted with sample collections.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00199/ full#supplementary-material FIGURE S1 | Location of the macrorefugia and microrefugia in the temperate conifers and broadleaved mixed forests inferred by Ye et al. (2017). The blue region indicates the range of the mixed forests following Wu (1980). The 35 occurrence records (black dots) obtained from the Global Biodiversity Information Facility database and from the 20 sampled populations (red dots) are shown.
FIGURE S2 | Illustration of the three scenarios (a-c) proposed for the divergence history of the three clusters of Schisandra chinensis obtained in DIYABC using eight nuclear microsatellites. N 1 , N 2 , and N 3 represent the effective population sizes of the east, north, and west genetic clusters inferred by Bayesian clustering. Posterior probability (PP) of the different scenarios is shown. Divergence times for the depicted events are labeled as t 1 and t 2 .
FIGURE S3 | K and LnP(D) obtained in the Structure analysis of the 20 Schisandra chinensis populations conducted for predefined group numbers (K = 1-10). The standard deviations of LnP(D) obtained from 10 independent runs for each predefined group size are also shown.
FIGURE S4 | Neighbor-joining tree of Schisandra chinensis populations using Nei's D a as the genetic distance of nSSRs. Bootstrap values above 70% are presented above nodes.
FIGURE S5 | Prior (red line) and posterior (green line) distribution of P, the proportion of multiple step mutations in the generalized stepwise model (a), t 2 , divergence time between west and north Structure clusters (b), t 1 , divergence time of all three Structure clusters (c), µ, mutation rate (d), and effective population sizes of west (e), north (f), and east (g) Structure clusters.
TABLE S1 | Primers used for chloroplast DNA and nuclear single copy gene sequencings in Schisandra chinensis. An annealing temperature of 54 • C was used in each case.