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, 2000; Qiu 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.

TABLE 1

Table 1. Details of sample location, sample size, and genetic diversity of chloroplast DNA (cpDNA), nuclear single copy gene (PEPC) and nuclear microsatellites (nSSRs) of Schisandra chinensis.

Materials and Methods

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.

Chloroplast and Nuclear DNA Sequence Analyses

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¹, 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² 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).

Microsatellite Data Analysis

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⁵ initial burn-in, followed by 1 × 10⁶ 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³ (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⁴ 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.

Results

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 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).

FIGURE 1

Figure 1. The most parsimonious network (a) and geographic distribution (b) of haplotypes based on four Schisandra chinensis chloroplast fragment, matK, ndhA, trnL-trnF and trnS-trnG. Outgroups represents Kadsura longipedunculata. Circle sizes are proportional to sample size of each population.

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).

FIGURE 2

Figure 2. The most parsimonious network (a) and geographic distribution (b) of haplotypes based on one Schisandra chinensis nuclear single copy gene, PEPC. Outgroups represents K. longipedunculata. Circle sizes are proportional to sample size of each population.

Neither significant phylogeographic structure (N_ST = 0.323 and G_ST = 0.311 in cpDNA and N_ST = 0.100 and G_ST = 0.065 in nSCG, P > 0.05) nor correlations between HR and latitude (r = 0.21, P = 0.38 in cpDNA and r = 0.21, P = 0.37 in nSCG) 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).

FIGURE 3

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.

FIGURE 4

Figure 4. Isolation by distance in Schisandra chinensis. Pairwised genetic distance (as measured by F_ST/(1 – F_ST)) is regressed onto the natural logarithm of pairwised geographic distance between all populations. Light gray dots represent population with sample size lower than 10.

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).

FIGURE 5

Figure 5. (A) Color-coded grouping of the 20 Schisandra chinensis populations according to the most likely K = 3 in Structure analysis. (B) Histogram of the Structure assignment test for all populations at the likely K = 2 and 3. Three genetic clusters (west, north and east) were also shown. Circle sizes are proportional to sample size of each population.

In the DIYABC analysis, the most possible scenario was that the west cluster diverged before the other two clusters (Supplementary Figure S5, PP = 0.90). The population divergence times between west and north clusters (t₂) and among all three clusters (t₁) were estimated as 3.97 × 10³ generations with a 95% highest-probability-density interval (HPD) of 0.79 × 10³ – 1.20 × 10³) and 7.09 × 10³ generations (95% HPD: 1.64 × 10³ – 12.40 × 10³), respectively. The generation time of S. chinensis is assumed to be five years. Therefore, the absolute t₂ was 19.85 kya (95% HPD: 3.95 – 60.00 kya) and the absolute t₁ was 35.45 kya (95% HPD: 8.20 – 62.00 kya). The mutation rate was estimated as 3.91 × 10⁻⁵/locus/generation. The estimated effective population size in west, north, and east clusters were 4.31 × 10⁴, 8.18 × 10⁴, and 7.13 × 10⁵, respectively (Table 2, Supplementary Figure S5). In Migrate analysis, almost all the estimated gene flows (4Nm) were very high (Table 3).

Ecological Niche Modeling

High ROC values (0.958 ± 0.012) indicated good accuracy of model predictions. During LGM, both models predicted range 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).

FIGURE 6

Figure 6. Potential species distribution of S. chinensis based on ENM at present (a), during the last glacial maximum using the MIROC 3.2 model (b) and in the CCSM4 model (c), and during the last interglacial (d).

TABLE 2

Table 2. Posterior median estimation and 95% highest posterior density interval (HPDI) for demographic parameters in the seventh divergence scenario of S. chinensis in DIYABC.

Discussion

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.

TABLE 3

Table 3. Estimates of gene flow (4Nm) among three groups of S. chinensis based on nuclear microsatellite structure division.

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

Data Availability

DNA sequences: deposited in Genbank under accessions MH580160–MH580199. All the chloroplast and nuclear sequences were deposited in TreeBASE (https://www.treebase.org/treebase-web/home.html) under submission number 23007 (http://purl.org/phylo/treebase/phylows/study/TB2:S23007).

Author Contributions

H-FW, LB, and J-PG conceived the study. Z-KZ contributed to the sampling. Z-KZ and J-WY collected and analyzed the data. J-WY, LB, and H-FW wrote the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (31210103911, 31770410, 31600301, and 31570381).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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₁, N₂, and N₃ 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₁ and t₂.

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₂, divergence time between west and north Structure clusters (b), t₁, 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.

TABLE S2 | Nuclear microsatellite locus, motif, allele size range, and annealing temperate in Schisandra chinensis.

TABLE S3 | Prior distributions for model parameters used in the seven divergence scenarios of Schisandra chinensis in DIYABC.

TABLE S4 | Selected bioclimatic variables with low correlations (r < 0.8) used in ecological niche modeling for Schisandra chinensis.

TABLE S5 | Haplotypes derived from four chloroplast DNA fragments in Schisandra chinensis.

TABLE S6 | Haplotypes derived from PEPC sequence in Schisandra chinensis.

TABLE S7 | Genetic diversity and genetic differentiation of eight nuclear microsatellite loci in Schisandra chinensis.

Footnotes

1. ^ http://www.codoncode.com/aligner/
2. ^ http://www.pierroton.inra.fr/genetics/labo/Software/Permut/
3. ^ http://data.gbif.org
4. ^ http://pmip2.lsce.ipsl.fr/

References

Ai, J., Wang, Y. P., Li, C. Y., Yan, T. L., and Guo, X. W. (2007). Pollen shape and pollination character of Schisandra chinensis. J. Jilin Agric. Univ. 29, 293–297.

Google Scholar

Antao, T., Lopes, A., Lopes, R. J., Beja-Pereira, A., and Luikart, G. (2008). LOSITAN: a workbench to detect molecular adaptation based on a Fst-outlier method. BMC Bioinformatics 9:323. doi: 10.1186/1471-2105-9-323

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, W. N., Liao, W. J., and Zhang, D. Y. (2010). Nuclear and chloroplast DNA phylogeography reveal two refuge areas with asymmetrical gene flow in a temperate walnut tree from East Asia. New Phytol. 188, 892–901. doi: 10.1111/j.1469-8137.2010.03407.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bandelt, H. J., Forster, P., and Röhl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036

PubMed Abstract | CrossRef Full Text | Google Scholar

Bao, L., Ayijiamali, K., Bai, W. N., Chen, R. Z., Wang, T. M., Wang, H. F., et al. (2015). Contributions of multiple refugia during the last glacial period to current mainland populations of Korean pine (Pinus koraiensis). Sci. Rep. 5:18608. doi: 10.1038/srep18608

PubMed Abstract | CrossRef Full Text | Google Scholar

Beaumont, M. A., Zhang, W. Y., and Balding, D. J. (2002). Approximate Bayesian computation in population genetics. Genetics 162, 2025–2035. doi: 10.1515/sagmb-2012-0050

PubMed Abstract | CrossRef Full Text | Google Scholar

Beerli, P. (2006). Comparison of Bayesian and maximum-likelihood inference of population genetic parameters. Bioinformatics 22, 341–345. doi: 10.1093/bioinformatics/bti803

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, X. Y., Herzschuh, U., Ni, J., Zhao, Y., and Boehmer, T. (2015). Spatial and temporal distributions of major tree taxa in eastern continental Asia during the last 22,000 years. Holocene 25, 79–91. doi: 10.1111/oik.01525

CrossRef Full Text | Google Scholar

Clark, P. U., and Mccabe, A. M. (2009). The last glacial maximum. Science 325, 710–714. doi: 10.1126/science.1172873

PubMed Abstract | CrossRef Full Text | Google Scholar

Cornuet, J. M., and Luikart, G. (1996). Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics 144, 2001–2014. doi: 10.3892/ijo_00000551

PubMed Abstract | CrossRef Full Text | Google Scholar

Cornuet, J. M., Pudlo, P., Veyssier, J., Dehnegarcia, A., Gautier, M., Leblois, R., et al. (2014). DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30, 1187–1189. doi: 10.1093/bioinformatics/btt763

PubMed Abstract | CrossRef Full Text | Google Scholar

Evanno, G., Regnaut, S., and Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Excoffier, L., Foll, M., and Petit, R. J. (2009). Genetic consequences of range expansions. Annu. Rev. Ecol. Syst. 40, 481–501. doi: 10.1146/annurev.ecolsys.39.110707.173414

CrossRef Full Text | Google Scholar

Falush, D., Stephens, M., and Pritchard, J. K. (2007). Inference of population structure using multilocus genotype data: dominant markers and null alleles. Mol. Ecol. Notes 7, 574–578. doi: 10.1111/j.1471-8286.2007.01758.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fawcett, T. (2006). An introduction to ROC analysis. Pattern Recognit. Lett. 27, 861–874. doi: 10.1016/j.patrec.2005.10.010

CrossRef Full Text | Google Scholar

Gavin, D. G., Fitzpatrick, M. C., Gugger, P. F., Heath, K. D., Francisco, R. S., Dobrowski, S. Z., et al. (2014). Climate refugia: joint inference from fossil records, species distribution models and phylogeography. New Phytol. 204, 37–54. doi: 10.1111/nph.12929

PubMed Abstract | CrossRef Full Text | Google Scholar

Goudet, J. (2001). FSTAT, a Program to Estimate and Test Gene Diversities and Fixation Indices (Version 2.9.3). Lausanne: Lausanne University.

Google Scholar

Guo, X. D., Wang, H. F., Bao, L., Wang, T. M., Bai, W. N., Ye, J. W., et al. (2014). Evolutionary history of a widespread tree species Acer mono in East Asia. Ecol. Evol. 4, 4332–4345. doi: 10.1002/ece3.1278

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, S. P., Yu, G., Takahara, H., and Prentice, I. C. (2001). Palaeovegetation (Communications arising): diversity of temperate plants in East Asia. Nature 413, 129–130. doi: 10.1038/35093166

PubMed Abstract | CrossRef Full Text | Google Scholar

Hewitt, G. M. (1996). Some genetic consequences of ice ages, and their role in divergence and speciation. Biol. J. Linn. Soc. 58, 247–276. doi: 10.1111/j.1095-8312.1996.tb01434.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hewitt, G. M. (2000). The genetic legacy of the Quaternary ice ages. Nature 405, 907–913. doi: 10.1038/35016000

PubMed Abstract | CrossRef Full Text | Google Scholar

Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., and Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. doi: 10.1002/joc.1276

CrossRef Full Text | Google Scholar

Jiang, Z. Y., Peng, Y. L., Hu, X. X., Zhou, Y. F., and Liu, J. Q. (2011). Cytoplasmic DNA variation in and genetic delimitation of Abies nephrolepis and Abies holophylla in northeastern China. Can. J. For. Res. 41, 1555–1561. doi: 10.1139/x11-069

CrossRef Full Text | Google Scholar

Kalinowski, S. T. (2005). hp-rare 1.0: a computer program for performing rarefaction on measures of allelic richness. Mol. Ecol. Notes 5, 187–189. doi: 10.1111/j.1471-8286.2004.00845.x

CrossRef Full Text | Google Scholar

Lafontaine, G. D., Ducousso, A., Lefèvre, S., Magnanou, E., and Petit, R. J. (2013). Stronger spatial genetic structure in recolonized areas than in refugia in the European beech. Mol. Ecol. 22, 4397–4412. doi: 10.1111/mec.12403

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y. X., and Zhang, S. J. (2017). Cultivation technology research of Aralia mandshrica and Schisandra chinensis. For. Sci. Technol. 42, 36–38.

Google Scholar

Liu, C. P., Tsuda, Y., Shen, H. L., Hu, L. J., Saito, Y., and Ide, Y. (2014). Genetic structure and hierarchical population divergence history of Acer mono var. Mono in south and northeast china. PLoS One 9:e87187. doi: 10.1371/journal.pone.0087187

PubMed Abstract | CrossRef Full Text | Google Scholar

Peakall, R., and Smouse, P. E. (2012). GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics 28, 2537–2539. doi: 10.1093/bioinformatics/bts460

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, S. J., Anderson, R. P., and Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259. doi: 10.1016/j.ecolmodel.2005.03.026

CrossRef Full Text | Google Scholar

Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959. doi: 10.1111/j.1471-8286.2007.01758.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Provan, J., and Bennett, K. D. (2008). Phylogeographic insights into cryptic glacial refugia. Trends Ecol. Evol. 23, 564–571. doi: 10.1016/j.tree.2008.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, Y. X., Fu, C. X., and Comes, H. P. (2011). Plant molecular phylogeography in China and adjacent regions: tracing the genetic imprints of Quaternary climate and environmental change in the world’s most diverse temperate flora. Mol. Phylogenet. Evol. 59, 225–244. doi: 10.1016/j.ympev.2011.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team (2013). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

Google Scholar

Rousset, F. (1996). Equilibrium values of measures of population subdivision for stepwise mutation processes. Genetics 142, 1357–1362. doi: 10.1111/j.1432-1033.1974.tb03647.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rozas, J., Sánchez-Delbarrio, J. C., Messeguer, X., and Rozas, R. (2003). DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics 19, 2496–2497. doi: 10.1079/9780851994758.0139

CrossRef Full Text | Google Scholar

Shields, C. A., Bailey, D. A., Danabasoglu, G., Jochum, M., Kiehl, J. T., Levis, S., et al. (2012). The low-resolution CCSM4. J. Clim. 25, 3993–4014. doi: 10.1175/JCLI-D-11-00260.1

CrossRef Full Text | Google Scholar

Takezaki, N., Nei, M., and Tamura, K. (2010). POPTREE2: software for constructing population trees from allele frequency data and computing other population statistics with Windows interface. Mol. Biol. Evol. 27, 747–752. doi: 10.1093/molbev/msp312

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L. J., Bao, L., Wang, H. F., and Ge, J. P. (2014). Study on the genetic diversity of Phyllodendron amurense. China Sciencepaper 9, 1397–1401. doi: 10.1002/ece3.1133

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S. H., Bao, L., Wang, T. M., Wang, H. F., and Ge, J. P. (2016). Contrasting genetic patterns between two coexisting Eleutherococcus species in northern China. Ecol. Evol. 6, 3311–3324. doi: 10.1002/ece3.2118

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W. T., Xu, B., Zhang, D. Y., and Bai, W. N. (2016). Phylogeography of postglacial range expansion in Juglans mandshurica (Juglandaceae) reveals no evidence of bottleneck, loss of genetic diversity, or isolation by distance in the leading-edge populations. Mol. Phylogenet. Evol. 102, 255–264. doi: 10.1016/j.ympev.2016.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Waters, J. M., Fraser, C. I., and Hewitt, G. M. (2013). Founder takes all: density-dependent processes structure biodiversity. Trends Ecol. Evol. 28, 78–85. doi: 10.1016/j.tree.2012.08.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Weir, B. S., and Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution 38, 1358–1370. doi: 10.1111/j.1558-5646.1984.tb05657.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wright, S. (1943). Isolation by distance. Genetics 28, 114–138. doi: 10.1007/BF02982775

CrossRef Full Text | Google Scholar

Wu, Z. Y. (1980). Chinese Vegetation. Beijing: Science Press.

Google Scholar

Ye, J. W., Jiang, T., Wang, H. F., Wang, T. M., Bao, L., and Ge, J. P. (2018). Repeated expansions and fragmentations linked to Pleistocene climate changes shaped the genetic structure of a woody climber, Actinidia arguta (Actinidiaceae). Botany 96, 19–31. doi: 10.1139/cjb-2017-0058

CrossRef Full Text | Google Scholar

Ye, J. W., Yuan, Y. G., Cai, L., and Wang, X. J. (2017). Research progress of phylogeographic studies of plant species in temperate coniferous and broadleaf mixed forests in Northeastern China. Biodiveris. Sci. 25, 1339–1349. doi: 10.17520/biods.2017265

CrossRef Full Text | Google Scholar

Yuan, L. C. (2007). Pollination Biology and Seed Dispersal of Several Species in Schisandraceae. Ph.D. dissertation, Institute of Botany, Chinese Academy of Sciencess, Beijing.

Google Scholar

Zeng, Y. F., Wang, W. T., Liao, W. J., Wang, H. F., and Zhang, D. Y. (2015). Multiple glacial refugia for cool-temperate deciduous trees in northern East Asia: the Mongolian Oak as a case study. Mol. Ecol. 24, 5676–5691. doi: 10.1111/mec.13408

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Niu, Y. B., Yan, L., Cui, Z. J., Liu, C. C., and Mu, K. H. (2008). Late Pleistocene glaciation of the Changbai Mountains in northeastern China. Chin. Sci. Bull. 53, 2672–2684. doi: 10.1007/s11434-008-0347-9

CrossRef Full Text | Google Scholar

Zhang, X. S. (2007). Vegetation Map of China and Its Geographic Pattern. Beijing: Geological Publishing House.

Google Scholar

Zhao, C., Wang, C. B., Ma, X. G., Liang, Q. L., and He, X. J. (2013). Phylogeographic analysis of a temperate-deciduous forest restricted plant (Bupleurum longiradiatum Turcz.) reveals two refuge areas in China with subsequent refugial isolation promoting speciation. Mol. Phylogenet. Evol. 68, 628–643. doi: 10.1016/j.ympev.2013.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar