Comparison of genetic variation between northern and southern populations of Lilium cernuum (Liliaceae): Implications for Pleistocene refugia

Abstract

The so-called “Baekdudaegan” (BDDG), a mountain range that stretches along the Korean Peninsula, has been recently proposed as a major “southern” glacial refugium for boreal or temperate plant species based on palaeoecological and, especially, genetic data. Genetic studies comparing genetic variation between population occurring on the BDDG and more northern ones (i.e. in NE China and/or in Russian Far East) are, however, still too few to draw firm conclusions on the role of the BDDG as a refugium and a source for possible northward post-glacial recolonizations. In order to fill this gap, we selected a boreal/temperate herb, Lilium cernuum, and compared levels of allozyme-based genetic diversity of five populations from NE China with five populations from South Korea (home of its hypothesized refuge areas). As a complementary tool, we used the maximum entropy algorithm implemented in MaxEnt to infer the species’ potential distribution for the present time, which was projected to different past climate scenarios for the Last Glacial Maximum (LGM). Permutation tests revealed that Korean populations harbored significantly higher levels of within-population genetic variation than those from NE China (expected heterozygosity = 0.173 vs. 0.095, respectively). Our results suggest that the lowered levels of genetic diversity in NE Chinese populations might be due to founder effects associated with post-glacial migration from southern regions. Congruent with genetic data, past distribution models showed higher probability of occurrence in southern ranges than in northern ones during the LGM. In addition, a positive correlation was detected between the expected heterozygosity and environmental LGM suitability. From a conservation perspective, our results further suggest that the southern populations in South Korea may be particularly worthy of protection.

About the Authors:

Mi Yoon Chung

Contributed equally to this work with: Mi Yoon Chung, Son Hai Vu

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

¶‡ Co-first authors.

Affiliation: Department of Biology and the Research Institute of Natural Science, Gyeongsang National University, Jinju, Republic of Korea

Son Hai Vu

Contributed equally to this work with: Mi Yoon Chung, Son Hai Vu

Roles Conceptualization, Data curation, Validation, Visualization, Writing – original draft

¶‡ Co-first authors.

Affiliation: Department of Biology and the Research Institute of Natural Science, Gyeongsang National University, Jinju, Republic of Korea

Jordi López-Pujol

Roles Conceptualization, Formal analysis, Methodology, Writing – original draft, Writing – review & editing

Affiliation: BioC-GReB, Instituto Botánico de Barcelona, IBB, CSIC-ICUB, Barcelona, Spain

Sonia Herrando-Moraira

Roles Formal analysis, Methodology, Writing – review & editing

Affiliation: BioC-GReB, Instituto Botánico de Barcelona, IBB, CSIC-ICUB, Barcelona, Spain

Sungwon Son

Roles Conceptualization, Data curation, Investigation, Project administration, Resources

Affiliation: Plant Conservation Division, Korea National Arboretum, Pocheon, Republic of Korea

Gang Uk Suh

Roles Conceptualization, Data curation, Investigation, Project administration, Resources

Affiliation: Plant Conservation Division, Korea National Arboretum, Pocheon, Republic of Korea

Hoa Thi Quynh Le

Roles Data curation, Formal analysis, Methodology

Affiliation: Department of Biology and the Research Institute of Natural Science, Gyeongsang National University, Jinju, Republic of Korea

Myong Gi Chung

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Department of Biology and the Research Institute of Natural Science, Gyeongsang National University, Jinju, Republic of Korea

ORCID logo http://orcid.org/0000-0002-1283-3574

Introduction

Unglaciated regions and/or regions that provided relatively stable environmental conditions during the Last Glacial Maximum (LGM; ca. 21,000 yr BP) are generally regarded as candidates for glacial refugia for plant species, at least in the Northern Hemisphere [1–4]. Therefore, it can be anticipated that populations located within refugia would show high levels of genetic diversity, as plant species would have endured the Pleistocene harsh climatic conditions by maintaining relative large populations there. Furthermore, plant populations in refugial areas would have had a longer demographic history than those that are result of post-glacial colonization because the latter would have experienced recurrent bottlenecks/founder effects [5]. In the formerly glaciated areas, on the contrary, low levels of genetic diversity are expected, as populations are of recent origin (originated through recolonization following deglaciation and, thus, suffering from founder effects).

In agreement with these expectations, a large body of studies—especially conducted in Europe—have revealed a common pattern of latitudinal decrease of genetic variation, which is often referred to as “southern richness vs. northern purity” [2,6,7,8,9]. Such studies are also abundant in North America (e.g. [10,11]), especially in the eastern coast of the United States where there is one of the most studied refugia of the world, the Southern Appalachians (see [12] for a review). There are many examples of higher levels of genetic diversity in plant populations of the Southern Appalachians compared to more northern conspecific ones, with several cases showing a south-to-north decrease of genetic diversity (e.g. [13,14]). In addition, several congeneric comparisons have consistently revealed lower genetic variability in northern species compared to their southern congeners [15–18].

On the Korean Peninsula, the so-called “Baekdudaegan” (hereafter the “BDDG”) is undeniably one of the most important mountain systems in East Asia in terms of length, ecological value, and sacredness [5,12]. More specifically, the BDDG is the longest mountain chain of NE Asia, stretching ca. 1625 km in a north-to-south direction, and it is the major center of plant diversity of the Korean Peninsula, harboring over 1500 taxa. In addition, it contains sites holy to all the major religions present in the country, namely Shamanism, Buddhism, Daoism, Neo-Confucianism, and even Christianism [5,12]. More significantly, it has been recently proposed as a major Pleistocene refugium for the boreal and temperate flora of NE Asia based on both palaeoecological and genetic evidence [5]. First, all the high-resolution pollen records available for the peninsula unambiguously suggest that most of the BDDG was covered by boreal or mixed boreal/temperate forests, a scenario also supported by several palaeovegetation reconstructions ([19–22]; but see [23,24]). Second, a series of allozyme genetic diversity studies conducted in boreal and/or temperate plant species mostly centered in the southern part of the BDDG has revealed a pattern of high within-population and low to moderate among-population genetic variability (Table 1 in [5]), suggesting that the BDDG acted as a refugium for these plant species. However, there are still very few comparisons between populations of boreal and/or temperate plant species occurring on the BDDG and populations farther north (i.e. in NE China and/or in Russian Far East).

Broadly, two scenarios could be proposed in relation to the levels of within-population genetic variation of boreal/temperate plant species that are native to the BDDG but that also occur in more northern latitudes. If a given species vanished from northern areas during the LGM and then recolonized the area through post-glacial migration from its putative refugia in the BDDG (see below), we may expect diminished levels of genetic diversity in populations located at more northerly latitudes, due to founder effects. Alternatively, if the species would have also endured the LGM in situ in northern areas, we may expect similar levels of genetic diversity in northern and southern populations.

In this study, we tested these alternative hypotheses with the boreal/temperate and shade-tolerant herb Lilium cernuum Kom. (Liliaceae). We believe that this is a good study system because of the following reasons. First, the species is distributed from Primorsky Krai in Russian Far East, Jilin and Liaoning provinces in NE China to the south of the Korean Peninsula ([25]; Fig 1); that is, the species’ range covers a latitudinal stretch of nearly 10°. And, second, a previous study conducted in the southern section of the BDDG and in its main branch (the so-called “Nakdongjeongmaek”) reported high within-population genetic variation [26,27], suggesting that this region was a putative glacial refugium for the species. We have selected a set of populations from the BDDG and a set of populations located much further north, in Jilin (NE China) and compared their genetic variability. As a complementary tool to the genetic study, we have used ecological niche modeling (ENM) to get insights into the palaeodistribution (at the LGM) of L. cernuum and to see whether there is congruence between ENM and genetic data.

[Figure omitted. See PDF.]

Fig 1. Location of sampled populations of Lilium cernuum (black circles, LC-1 to LC-5 from NE China; LC-6 to LC-10 from South Korea).

Grey dashed line indicates the approximate location of the ridge of the main mountain range of the country, the so-called “Baekdudaegan” (BDDG), which runs north to south along the Korean Peninsula. The blue circles are precise occurrence records (which correspond to those used in ecological niche modeling; see Fig 5), whereas the red circles are approximate localities of L. cernuum. Presence records, regardless of their precision, have been mostly obtained from herbarium specimens but also from databases, scientific papers, and grey literature.

https://doi.org/10.1371/journal.pone.0190520.g001

Materials and methods

Study plant

Lilium cernuum is a herb 30–80 cm tall, with scattered, narrowly linear leaves. Pale purple-red flowers (ca. 4.0 cm long) are solitary or arranged in racemes of 2–7, nodding, and fragrant, and open in July. Fruits (capsules) are globose to obovoid, 1.5 cm long. The breeding system and pollinators are unknown [27]. It is a diploid species (2n = 24; http://ccdb.tau.ac.il/Angiosperms/Liliaceae/Lilium/Lilium%20cernuum%20Ko) that occurs in thickets and sunny grassy slopes in the northern part of its range (NE China, Russian Far East; [25,28]) but under Quercus mongolica-dominated temperate deciduous forests in Korea [29]. It is regarded as threatened both in Russia [30] and China [31] but not in Korea.

Lilium cernuum is considered endangered in China [it is listed as “VU” (“Vulnerable”) in the Red List of 2013] but, to our knowledge, it is not protected because it was not included in the laws of 1984 (the so-called “National List of Rare and Endangered Plant Species”) and 1999 (the “Catalogue of the National Protected Key Wild Plants”). In South Korea, L. cernuum was delisted as an endangered species in 2012. Now it has been listed as “LC” (“Least Concern”), i.e., not considered threatened.

We selected five populations from NE China (n = 207) and five populations from South Korea (n = 217), and these were sampled from areas of similar size to prevent area effects (Fig 1). None of these populations were inside nature reserves. In each population, we collected samples randomly from adults (flowering individuals). To minimize the damage to these lilies, we collected only one leaf per individual. Korean populations (LC-6 to LC-10) correspond to the LC-2 to LC-6 populations of the former study of 2014 [27], although the sample sizes of the present study are larger (mean n = 43 vs. 34). In addition, the resolved allozymes are not the same (Mdh-1 and Mdh-2 are lacking in the present study as the resolution of banding patterns was not optimal for the NE China populations).

Enzyme electrophoresis

For enzyme extraction, we crushed leaf samples using chilled mortars and pestles and an extraction buffer [32] to solubilize and stabilize the enzymes. We absorbed enzyme extracts onto 4 × 6 mm paper wicks (Whatman 3MM chromatography paper) and we stored them at −70°C until needed. We conducted electrophoresis on 13% starch gels, with two buffer systems (Table 1). We stained gels for eight enzyme systems, which were alcohol dehydrogenase (ADH, E.C.1.1.1.1), diaphorase (DIA, E.C.1.6.99.1), fluorescent esterase (FE, E.C.3.1.1.1), isocitrate dehydrogenase (IDH, E.C.1.1.1.42), phosphoglucoisomerase (PGI, E.C.5.3.1.9), 6-phosphogluconate dehydrogenase (6PGD, E.C.1.1.1.44), phosphoglucomutase (PGM, E.C.2.7.5.1), and triosephosphate isomerase (TPI, E.C.5.3.1.1). We utilized stain recipes from [33], except for DIA [34]. We designated putative loci designated sequentially, with the most anodally migrating isozyme designated as 1 and the next 2. We also designated different alleles per locus sequentially by alphabetical order.

[Figure omitted. See PDF.]

Table 1. The allozyme loci assayed for Lilium cernuum and the buffer systems used to resolve thema.

Buffer system no. 34 was a modification [35] of the system 6 of [33]. The morpholine-citrate system (MC, pH 6.1) was from [36].

https://doi.org/10.1371/journal.pone.0190520.t001

Data analyses

To estimate genetic diversity and structure, we considered that a locus was polymorphic when two or more alleles were observed, regardless of their frequencies. We estimated the following genetic diversity parameters using the programs POPGENE [37] and FSTAT [38]: percent polymorphic loci (%P), mean number of alleles per locus (A), allelic richness (AR) using a rarefaction method that compensates uneven population sample sizes [39,40], observed heterozygosity (Ho), and Nei [41]’s unbiased gene diversity or Hardy–Weinberg (H–W) expected heterozygosity (He). Hereafter, the subscript “s” indicates species’ (or pooled samples) values, while the subscript “p” indicates population means.

To test for differences between populations in South Korea and those in NE China for observed statistics, OSx (AR, Ho, and He), we used a permutation scheme (999 replicates) by randomly allocating whole samples to the different groups, keeping the number of samples in each group constant and calculating differences between South Korean and northern Chinese populations for randomized statistics, RSx. We then obtained the P value of the test as the proportion of randomized data sets giving higher values for RSx than for OSx. These calculations were performed using FSTAT [38].

To test for recent decreases in effective population size (bottlenecks), we evaluated for individual loci the difference between the H–W He and the equilibrium heterozygosity (Heq) expected assuming mutation–drift equilibrium. These differences were evaluated using a sign test and a Wilcoxon sign-rank test conducted across loci under an infinite allele model using the program BOTTLENECK [42]. Since allelic diversity is generally lost more rapidly than He [43], recently bottlenecked populations are expected to exhibit an excess of H–W equilibrium He relative to Heq [44,45].

We estimated population-level FIS (inbreeding) and calculated its significance level (P values) by gene permutation tests (999 replicates) under the null hypothesis (FIS = 0) using the program SPAGeDi [46]. We also calculated Wright [47]’s FIS and FST over loci following [48]. These fixation indices measure the average deviation from H–W equilibrium of individuals relative to their local populations (FIS, a measure of local inbreeding) and local populations relative to the total population (FST, a measure of differentiation between local populations). The significance of multi-population FIS and FST estimates was determined by a permutation test (999 randomizations of alleles between individuals within samples and 999 randomizations of genotypes between populations, respectively). These calculations were performed using FSTAT [38]. Statistical significance of differences in FIS and FST between populations in southern Korea and those in NE China was determined as outlined above for AR, Ho, and He. To test for the influence of individuals within populations, populations within regions, and regions on the observed genetic variation, we conducted an analysis of molecular variance (AMOVA) using the program GenAlEx 6.5 [49].

A UPGMA (unweighted pair-group method using arithmetic averages) phenogram was generated from Nei et al. [50] genetic distance matrix with branch support produced by 1000 bootstrapping over loci, utilizing Populations 1.2.30 [51] and TreeView 1.6 [52]. We assessed the genetic structure by means of the Bayesian algorithm implemented in STRUCTURE 2.3.4 [53]. The program estimates the likelihood of the individuals being structured in a given number of genetic clusters (or genetic populations, K). The admixture ancestry model with correlated allele frequencies was selected as an appropriate option for the analysis, as it is the best combination for dealing with species with presumably complex dynamics; individuals with mixed ancestry, events of migration, and populations with shared ancestry are presumably to occur in L. cernuum. The burn-in period and Markov Chain Monte Carlo (MCMC) were set to 50,000 and 500,000 iterations, respectively, and 10 replicates per K were run, with a priori grouping of individuals into populations. The most likely value of K was determined both by choosing the smallest K after the log probability of data [ln Pr(X|K)] values reached a plateau [53] and by the ΔK statistic of [54], with the aid of STRUCTURE HARVESTER [55]. For the most likely K, Clumpp 1.1.2b [56] was used to combine the results of the 10 replicates of the best K. To plot the output result produced by Clumpp, we used the program Distruct 1.1 [57].

Finally, we investigated whether there was isolation-by-distance by conducting a Mantel test between all pairwise FST/(1 –FST) values (FST was calculated following [48]) and the corresponding logarithm of pairwise geographical distances [58]. A linear regression model under the null hypothesis of no spatial genetic structure (regression slope, b = 0) was tested (by making 999 replicates) using the program Permute![59].

Ecological niche modeling (ENM)

Ecological niche modeling (ENM) was performed to evaluate the potential distribution of L. cernuum under present climatic conditions and to project it to the LGM. We employed the maximum entropy algorithm, as implemented in MaxEnt 3.3 [60]. The current distribution information for L. cernuum was obtained from: (1) major databases and information systems, including the Global Biodiversity Information Facility (GBIF; http://www.gbif.org/) and the Chinese Virtual Herbarium (http://www.cvh.ac.cn); (2) research articles, books, and grey literature (i.e. theses (e.g. [61]) and unpublished reports); (3) citizen science projects, such as iNaturalist (http://www.inaturalist.org), Chinese Field Herbarium (CFH; http://www.cfh.ac.cn), and Plant Photo Bank of China (PPBC; http://www.plantphoto.cn); and (4) the sampling sites of this study (Table 2). In total, after removing duplicate records within each pixel (2.5 arc-min, ca. 5 km), we obtained 58 presences. A set of 19 bioclimatic variables at 2.5 arc-min resolution covering the species distribution range (and neighboring areas) under current conditions (1950–2000) were downloaded from the WorldClim website (http://www.worldclim.org; [62]). After a pairwise correlation analysis, carried with the SDM Toolbox extension for ArcMap [63], we selected a smaller set of six (relatively) uncorrelated variables: mean diurnal range (bio2), isothermality (bio3), mean temperature of the wettest quarter (bio8), mean temperature of the coldest quarter (bio11), precipitation of the wettest month (bio13), and precipitation of the driest month (bio14). The selection of variables from pairs or groups of highly correlated (r ≥ |0.8|) ones was done on the basis of their relative contribution to the model (percent contribution, jackknife tests of variable importance) and the shape of their response curves.

[Figure omitted. See PDF.]

Table 2. Summary of genetic diversity measures and mean fixation values (FIS) for five Chinese populations and five Korean populations of Lilium cernuuma.

https://doi.org/10.1371/journal.pone.0190520.t002

The distribution model under current conditions was projected to the LGM using palaeoclimatic layers simulated by both the Community Climate System Model Version 4 (CCSM4; [64]), the Model for Interdisciplinary Research on Climate Earth System Model (MIROC-ESM; [65]), and the New Earth System Model of the Max Planck Institute for Meteorology (MPI-ESM-P: http://www.mpimet.mpg.de/en/science/models/mpi-esm/). Replicate runs (20) of MaxEnt using the “subsample” method were performed to ensure reliable results. Model performance was assessed using the area under the curve (AUC) of the receiver operating characteristic plot, with 25% of the localities randomly selected to test the model. AUC scores may range between 0.5 (randomness) and 1 (exact match), with those above 0.9 indicating a good performance of the model [66]. The MaxEnt jackknife analysis was used to evaluate the relative importance of the six bioclimatic variables employed, based on their gain values when used in isolation. In a simple way to reduce the uncertainty associated to the use of different global climate models (GCMs), we produced a LGM-ensemble map that identified all areas predicted as suitable that were common to the three LGM projections. To do this, firstly we converted the three continuous projection maps into binary (presence/absence) distribution maps, applying the maximum sensitivity plus specificity logistic threshold, which is very robust with all types of data [67]. Secondly, we overlapped the binary output maps with the Intersect Tool of ArcGIS. All ENM predictions were visualized in ArcGIS 10.2 (ESRI, Redlands, CA, USA).

Comparison between genetic diversity and habitat suitability

In biogeography, ENM is used as a complementary tool to the genetic markers in order to reconstruct species palaeodistributions, often showing good agreement (e.g. [68,69]). Herein, following a methodology that is based on [70], we analyzed the relationship between He and habitat suitability (measured as MaxEnt’s logistic output, a value that ranges from 0 to 1 and that can be cautiously interpreted as the predicted probability of presence; see [71]) by means of Pearson’s correlations, for all ten studied populations of L. cernuum. Correlations were performed separately for the present time and for the LGM; for this latter, the values of habitat suitability were obtained by averaging the values for each of the three palaeoclimatic scenarios (CCSM4, MIROC-ESM, and MPI-ESM-P).

Results

Genetic diversity in northern and southern populations of Lilium cernuum

Of the 12 putative loci surveyed for L. cernuum (Table 2), nine (Adh, Dia-1, Fe, Idh, 6Pgd-1, 6Pgd-2, Pgi-2, Pgm, and Tpi-1) were polymorphic across 10 populations, resulting in high levels of genetic variation in pooled samples (n = 410, %Ps = 75.0, As = 2.33, and Hes = 0.152; Table 2). Considerably higher levels of genetic variation were found within southern populations (South Korea) compared to northern populations (NE China): %Pp = 56.7 vs. 43.3, Ap = 1.84 vs. 1.63, AR = 1.79 vs. 1.60, and Hep = 0.173 vs. 0.095 (Table 2). Accordingly, we found significant differences for AR, Hop, and Hep between South Korean and NE Chinese populations (one-sided P-values were 0.043, 0.002, and 0.006, respectively). Regarding the allelic composition, we identified 22 alleles in the northern populations and 28 alleles in the southern ones; six alleles were exclusive to the Korean populations and none of them to the Chinese ones (i.e. the Chinese populations could be considered a subset of the allelic diversity of the Korean ones; S1 Table). Among the 10 studied populations, there were no statistically significant signals of recent bottlenecks (S2 Table).

Inbreeding and population genetic structure

We found a significant deficiency of heterozygotes (at the 0.05 level) relative to H–W expectations in all but LC-4 and LC-7 populations (Table 2). These results, as well as the significant multi-population-level FIS (FIS = 0.185, P = 0.001; Table 2), indicated a substantial deficit of heterozygotes within populations. The values of pooled multi-population FIS for northern and southern populations were significantly greater than zero (FIS = 0.280, P = 0.001 vs. FIS = 0.136, P = 0.001; Table 2), but the former was significantly larger than the latter (P = 0.001). Deviation from H–W expectations due to allele frequency differences between populations were, despite low, significantly different from zero for the two regions (FST = 0.108, P = 0.001 for northern populations vs. FST = 0.094, P = 0.001 for southern populations). However, we found no statistically significant differences between the two groups of populations (P = 0.571). Across 10 populations, the value of FST was 0.133 (P = 0.001). Consistent with these results, the AMOVA indicated that most of the variation resided within populations (84%), with variation attributable to differences between regions (South Korea vs. NE China) and among populations within regions being only 7% and 9%, respectively (Table 3).

[Figure omitted. See PDF.]

Table 3. Analysis of molecular variance (AMOVA) of Lilium cernuum populationsa.

https://doi.org/10.1371/journal.pone.0190520.t003

The UPGMA phenogram (Fig 2) revealed a clear genetic separation between northern and southern populations. The best clustering scheme of STRUCTURE (K = 2, according to both the ln Pr(X|K) and the ΔK statistic; S1 Fig), agreed with the UPGMA phenogram, although some populations, such as LC-1, and especially LC-10, showed a high degree of admixture (Fig 3).

[Figure omitted. See PDF.]

Fig 2. UPGMA phenogram based on Nei et al. [50] genetic distances between populations of Lilium cernuum.

Note that LC-1 to LC-5 are from NE China, whereas LC-6 to LC-10 are from South Korea. Bootstrap values greater than 50% from 1000 replicates were found only for three branches.

https://doi.org/10.1371/journal.pone.0190520.g002

[Figure omitted. See PDF.]

Fig 3. Results of STRUCTURE for all Lilium cernuum studied individuals.

Assignation of individuals to genetic clusters is at K = 2.

https://doi.org/10.1371/journal.pone.0190520.g003

Finally, we found a significant positive linear relationship between pairwise FST/(1 –FST; see Table 4 for pairwise FST estimates) and logarithm of pairwise linear geographic distances (r = 0.312, P = 0.041; Fig 4).

[Figure omitted. See PDF.]

Table 4. Pairwise FST values between 10 populations of Lilium cernuum.

The significant FST values are in bold.

https://doi.org/10.1371/journal.pone.0190520.t004

[Figure omitted. See PDF.]

Fig 4. Differentiation between populations of Lilium cernuum.

Multilocus estimates of pairwise differentiation of FST/(1 − FST) are plotted against logarithm (ln) of pairwise geographic distances in meters. A significant relationship was found for L. cernuum (r = 0.312, P = 0.041).

https://doi.org/10.1371/journal.pone.0190520.g004

Ecological niche modeling and congruence with genetic data

The AUC scores averaged across 20 runs were high (> 0.93) for the four models built (present only, and present projected to CCSM4, MIROC-ESM, and MPI-ESM-P), which supported their predictive power. According to the MaxEnt jackknife tests of variable importance, the most informative variables for all models were bio11, bio14, and bio8 (in this order). The present day distributional predictions for L. cernuum were largely congruent with the known occurrences; the large areas of North Korea without occurrences are mostly due to the lack of floristic research since the 1950s and not to the actual absence of the plant there. Projections of the species niche to the LGM climate produced considerably different maps of probability of occurrence compared to the present time. All three models, although variable, showed a similar pattern of southwards displacement of the suitable areas; such areas included the southern portion of North Korea, most of South Korea, a large part the exposed East China Sea shelf (including the Yellow Sea), large parts of the Korea Strait, and southern Japan (Fig 5). Indeed, the LGM-ensemble map hardly showed suitable areas for the presence of L. cernuum north of 40°N (S2 Fig).

[Figure omitted. See PDF.]

Fig 5. A comparison of potential habitats for Lilium cernuum under the present climate and three climatic scenarios for the Last Glacial Maximum (LGM, ca. 21,000 yr BP).

(A) Potential habitats under the present climate; (B) Potential habitats expected for the LGM using the CCSM4 model, (C) the MIROC-ESM model, and (D) the MPI-ESM-P model. The darker color indicates a higher probability of occurrence. Maps were generated using the software ArcGIS 10.2. The reconstructed LGM coastlines are represented as dashed lines. Blue circles are the precise occurrence records of L. cernuum (at 2.5 arc-min resolution) used to build the models.

https://doi.org/10.1371/journal.pone.0190520.g005

As anticipated, ENM and genetic data showed good agreement, as there was a positive and statistically significant correlation between He and habitat suitability for the studied populations of L. cernuum, for both the present time and the LGM (Fig 6). Notably, there was no any northern population (NE China) that showed a higher probability of occurrence relative to southern ones (South Korea) for both time frames (Fig 6).

[Figure omitted. See PDF.]

Fig 6. Pearson’s correlations between expected the expected heterozygosity (He) and habitat suitability (measured as MaxEnt’s logistic output) for all Lilium cernuum studied individuals, for (A) the present time, and (B) the LGM.

https://doi.org/10.1371/journal.pone.0190520.g006

Discussion

Patterns of genetic diversity in northern and southern populations: Implications for glacial/post-glacial plant dynamics

Levels of within-population genetic diversity found in L. cernuum are in agreement with our prediction that the southern populations (those from South Korea) harbor greater genetic variation than northern ones (NE China) (%Pp = 56.7 vs. 43.3, Ap = 1.84 vs. 1.63, AR = 1.79 vs. 1.60, and Hep = 0.173 vs. 0.095; Table 2). In addition, the Chinese populations contain a subset of the allelic diversity of the Korean ones; no allele is exclusive to the Chinese populations (S1 Table). Our results suggest that L. cernuum was probably not present in situ in NE China during the LGM and, thus, the current northern populations are the result of post-glacial recolonization from southern glacial refugia. The higher genetic polymorphism detected at the studied populations from South Korea indicates that such refugia would have been located along the BDDG, which mostly constitutes the current southern range of the species. This scenario of survival in southern refugia followed by northwards migration is consistent with the reconstructed palaeodistribution of L. cernuum for the LGM, a time when the only suitable habitats for the species would have been limited to areas below 40°N or so (Fig 5 and S2 Fig). Most of the current range of the species in South Korea (Fig 1) is included within areas of medium to high probability of occurrence in the LGM models, with the relative exception of MIROC-ESM (Fig 5). Therefore, L. cernuum would have had the chance to persist in situ in South Korea during the LGM and, perhaps, at other glacial maxima. In contrast, the current northern range of the species (the Chinese provinces of Jilin and Liaoning, and Primorsky Krai in Russian Far East) would not have harbored climatically suitable areas for L. cernuum during the Pleistocene glacial periods. Given the positive relationship between He and habitat suitability in L. cernuum, ENM might constitute a simple and cost-effective surrogate of genetic diversity studies. Although not explicitly tested, some studies with boreal/temperate plant species that are also native to our study region (Korea plus NE China and Russian Far East) also show an apparent congruence between genetic diversity levels and habitat suitability, identifying the mountains of South Korea as a putative glacial refugium [72,73]. Of these, it is especially relevant the study carried out with Quercus mongolica [73], the dominant species of the plant communities in which L. cernuum usually occurs.

The results from the present study are in agreement with a recently published review of the region’s Pleistocene biogeography [5], in which the authors demonstrated that the BDDG served as a Pleistocene refugium for (but not exclusively) the boreal and temperate flora of NE Asia; such statement is based on both palaeoecological and genetic evidences. In this review [5], the authors summarized the extant palaeovegetation reconstructions for NE Asia, which mostly suggest that the vegetation belts suffered major range displacements, nearly comparable to those occurred in other regions of the world even though ice sheets were not present in the vast majority of NE Asia during the LGM (see Fig 1 in [5]). In NE China and in the southern parts of the Russian Far East (i.e. Primorsky Krai), LGM vegetation was either treeless or a mixture of pure boreal forests and non-forest vegetation (tundra, shrub-tundra, or polar deserts) [20–24,74,75]. A few pollen records available for the areas that correspond to the northern range of L. cernuum also indicate that the vegetation there adapted to much colder conditions than today. For instance, in a couple of sites at latitude ca. 43°N, in the southernmost tip of Primorsky Krai, boreal forests dominated by Betula (accompanied by Larix, Picea, and Pinus) would have occurred [76]. Similarly, in the contiguous Jilin Province of China (a site located in a hilly area at ca. 100 km of the Korean border and ca. 250 km west from LC-2), vegetation probably dominated by cold deciduous and needle-leaved forests (Betula and Picea with the presence of Abies and Pinus) would have occurred [76]. Taken all together, it is highly likely that boreal/temperate herbaceous elements like L. cernuum would have not been present in situ in areas such as NE China or Russian Far East during the LGM.

More favorable habitats for this lily, instead, would have occurred further south. According to the tentative LGM vegetation map for the Korean Peninsula provided by [5], non-arboreal vegetation would have been limited to the northern half of North Korea, whereas the rest of the peninsula would have been covered by boreal and mixed (boreal/temperate) forests, with temperate ones at the southernmost tip. Although pollen records within or around the BDDG are very scarce, it is likely that the central to northern sections would have sustained boreal forests, whereas mixed forests were probably ubiquitous in the central to southern parts of this mountain range (see [5] and references therein for more details).

Although the results must be interpreted with caution because of the limited sampling in the northern range of L. cernuum (neither populations from Liaoning Province nor Russian Far East have been studied), there is a clear pattern of higher genetic diversity for the Korean populations compared to the Chinese ones; such results are fully compatible with the genetic studies conducted for species endemic or native to the region, which have been recently compiled by [5]. First, as noted in the Introduction, the meta-analysis of the allozyme literature published for plant species occurring on the BDDG revealed a pattern of high within-population genetic diversity (Hep = 0.159; [5]), much larger than the “classical” values of reference (e.g. [77,78]). The value of Hep for the Korean populations of L. cernuum (0.173) is close to the mean value reported by [5], whereas that for the Chinese populations is much lower (Hep = 0.095). Second, the still small but growing body of studies on genetic diversity and phylogeography at regional level including populations located in southern refugia on the Korean Peninsula and populations located further north, are mostly consistent with the “southern richness vs. northern purity” pattern. That is, Korean populations harbor more genetic diversity, ancestral haplotypes, and/or significant amounts of unique haplotypes/alleles (see [5] and references therein). Very similar to our study, for example, [79] examined the genetic variation of the boreal tree Pinus koraiensis across the NE Asian continent (South Korea—with three of the four populations located in the BDDG, NE China and Russian Far East) by means of allozyme and RAPD markers. The authors found that levels of within-population genetic variation significantly decreased with increasing latitude. More recently, [72] found a similar pattern of latitudinal decrease of both genetic diversity and genetic singularity across NE Asia for the boreal/temperate tree Acer mono; both cpDNA and microsatellites suggested that the BDDG was the origin for the recolonization of NE China, a scenario also recovered from ENM [72].

Levels of genetic diversity and structure of L. cernuum and implications for conservation

The Korean populations of L. cernuum maintain high levels of genetic variation when compared to other herbaceous perennials mainly occurring on the BDDG (Table 1 in [5]), as also reported in a previous study [27]. Therefore, Korean populations seem not to be of conservation concern, although it should be taken into account that several populations of L. cernuum in South Korea are small and discontinuously distributed along the ridge of the BDDG (M. Y. Chung and M. G. Chung, pers. obs.). However, other populations, such as LC-7 (Mt. Deokhang), would enjoy large effective population sizes (Ne) [26]. The low FST values both at local (0.019 at a scale of <1 km; [26]), country (0.119 at a scale of 15–300 km; [27]) and regional level (0.094 at a scale of 14–730 km; present study) suggest high historical gene flow between populations. This result is an expected outcome given the role of the BDDG as Pleistocene refugium but also the species’ high potential for seed dispersal [26].

Even being relatively isolated geographically at the southern tip of the species’ distribution, the Korean population LC-10, that occurs at a high-elevation mountain within the BDDG (Mt. Sobaek, ca. 1220 m), has the lowest amount of genetic diversity among the five southern populations (but still higher than the northern populations). On the other hand, the populations LC-7 (Mt. Deokhang) and LC-8 (Yongyeon Cave), located approximately at the center of the current range for the species within South Korea, maintain the highest levels of genetic diversity among the studied populations. As mentioned above, populations around LC-7 maintain large effective population sizes on the order of hundreds [26]. These findings are in accordance with the “center-periphery” model (reviewed in [80]) which proposes that marginal populations are genetically less diverse than core ones; such model is based in the assumption that both effective population sizes (Ne) and rates of gene flow (m) are higher at the range center compared to range margins [80]. Historical processes are increasingly recognized as major factors creating such pattern (e.g. the influence of the Pleistocene glacial/interglacial cycles on species’ demography). The low levels of polymorphism detected for the population LC-10, as also reported for the population of Mt. Gaji for the same species [27], fit well the genetic pattern expected for the “rear-edge populations” hypothesis of [81]. According to these authors, populations located at the low-latitude margins of species ranges are expected to show reduced within-population genetic variation (see also [82]) and, thus, may merit specific conservation measures.

In contrast to some populations distributed along the ridge of the BDDG, populations from NE China and Russian Far East that are close to the border with North Korea are relatively large and continuously distributed (S. Son and G. U. Suh, pers. obs.). Given present-day population sizes, population genetics theory predicts that genetic diversity should be higher in the populations from NE China compared to populations from South Korea. However, our results show the opposite pattern (most likely due to the northwards post-glacial expansion of L. cernuum), stressing the importance of empirical genetic studies in plant species [83]. In spite of their lower polymorphism compared to their southern counterparts, the studied NE Chinese populations of L. cernuum should not be regarded as genetically extremely depauperate. Reference values for plants in general (Hep = 0.116; [78]), plants regionally-distributed (Hep = 0.118; [77]), and short-lived herbaceous perennials (Hep = 0.096; [77]) are only slightly higher than the observed values for the Chinese populations of L. cernuum (mean Hep = 0.095).

If we use the equilibrium equation between FIS and t (outcrossing rates) or s (selfing rates) (FIS = s/(2 − s) or FIS = (1 − t)/(1 + t); [84,85]), we obtain a t of 0.688 for L. cernuum, suggesting a mixed mating system in the species. The occurrence of selfing, however, might not be the only factor explaining the observed heterozygote deficiency at most populations, as other factors, including non-random mating with relatives (biparental inbreeding) and Wahlund effect (population subdivision), may account. Because NE Chinese populations of L. cernuum have been observed to be large and continuous (S. Son and G. U. Suh, pers. obs.), the higher levels of inbreeding for NE Chinese (FIS = 0.280) compared to South Korean (FIS = 0.136) populations are probably indicating different reproductive behavior between the two groups instead of population subdivision. Detailed field work would be necessary to see whether there is reduced pollinator service in the northern populations compared to the southern ones. The significantly larger FIS estimates found in the NE Chinese populations compared to those in the southern ones may have conservation implications, as inbreeding may lead to reduced fitness [86].

The loss or reduction of genetic diversity is of particular concern for conservationists and managers because it may lead directly to the decrease of biological fitness of populations. Thus, maintenance of genetic diversity is often viewed as a prerequisite for adaptation to environmental changes [86] and, thus, coping with the Sixth Mass Extinction. Understanding the genetic make-up of populations occurring on the BDDG is crucial because they have been regarded as a reservoir and at the same time the genetic stock for post-glacial recolonization for temperate and boreal plant species in NE Asia [5]. Our results are consistent with such scenario, emphasizing the conservation value of southern populations, which could be used for future translocation or reinforcement actions.

Conclusions

Our results from genetic and ENM analyses support the hypothesis that L. cernuum vanished from the Chinese provinces of Jilin and Liaoning, and Primorsky Krai in Russian Far East during the LGM, and that the species was only able to persist in southern glacial refugia, presumably located along the BDDG in South Korea. Our results are in concordance with previous studies showing that a pattern of survival in southern refugia followed by post-glacial recolonizations would have been common for the boreal and temperate flora of NE Asia [5,12]. Our results further suggest that the populations of L. cernuum in South Korea are particularly worthy of protection.

Supporting information

[Figure omitted. See PDF.]

S1 Fig. Estimation of the most likely value of K, as determined by the ln Pr(X|K) and the ΔK statistics (left and right, respectively), with the aid of Structure Harvester.

https://doi.org/10.1371/journal.pone.0190520.s001

(DOCX)

S2 Fig. LGM-ensemble map obtained from the projections with the three global climate models employed (CCSM4, MIROC-ESM, and MPI-ESM-P).

https://doi.org/10.1371/journal.pone.0190520.s002

(DOCX)

S1 Table. Allele frequencies of nine polymorphic loci of the 10 populations of Lilium cernuum.

https://doi.org/10.1371/journal.pone.0190520.s003

(XLS)

S2 Table. Results of statistical tests for evidence of recent population bottlenecks in Lilium cernuum.

Numbers reported are P values of sign and Wilcoxon sign-rank tests conducted using the program BOTTLENECK.

https://doi.org/10.1371/journal.pone.0190520.s004

(DOCX)

Acknowledgments

The authors thank Beom Jin Shim, Myeong Soon Park for laboratory assistances; Hoa Thi Quynh Le for statistical analyses; Jae Min Chung and Jin Seok Kim for helping us in locating populations of Lilium cernuum. The authors also appreciate the tuition-fee exemption program by Gyeongsang National University Graduate School during S. H. Vu’s study as a master student at M.G. Chung’s laboratory from September 2015 to July 2017.

Citation: Chung MY, Vu SH, López-Pujol J, Herrando-Moraira S, Son S, Suh GU, et al. (2018) Comparison of genetic variation between northern and southern populations of Lilium cernuum (Liliaceae): Implications for Pleistocene refugia. PLoS ONE 13(1): e0190520. https://doi.org/10.1371/journal.pone.0190520

References

1. Hewitt GM. Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc. 1996; 58: 247–276.

2. Hewitt GM. The genetic legacy of the Quaternary ice ages. Nature. 2000; 405: 907–913. pmid:10879524

3. Comes HP, Kadereit JW. The effect of quaternary climatic changes on plant distribution and evolution. Trends Plant Sci. 1998; 3: 1360–1385.

4. Gavin DG, Fitzpatrick MC, Gugger PF, Heath KD, Rodríguez-Sánchez F, Dobrowski SZ, et al. Climate refugia: joint inference from fossil records, species distribution models and phylogeography. New Phytol. 2014; 204: 37–54. pmid:25039238

5. Chung MY, López-Pujol J, Chung MG. The role of the Baekdudaegan (Korean Peninsula) as a major glacial refugium for plant species: A priority for conservation. Biol Conserv. 2017; 206: 236–248.

6. Taberlet P, Fumagalli L, Wust-Saucy A-G, Cosson J-F. Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998; 7: 453–464. pmid:9628000

7. Hewitt GM. Genetic consequence of climatic oscillations in the Quaternary. Philos Trans R Soc B-Biol Sci. 2004; 359: 183–195.

8. Hu FS, Hampe A, Petit RJ. Paleoecology meets genetics: deciphering past vegetational dynamics. Front Ecol Environ. 2009; 7: 371–379.

9. Nieto-Feliner G. Patterns and processes in plant phylogeography in the Mediterranean Basin. Perspect Plant Ecol Evol Syst. 2014; 16: 265–278.

10. Soltis DE, Gitzendanner MA, Strenge DD, Soltis PS. Chloroplast DNA intraspecific phylogeography of plant from the Pacific Northwest of North America. Plant Syst Evol. 1997; 206: 353–373.

11. Soltis DE, Morris AB, McLachlan JS, Manos PS, Soltis PS. Comparative phylogeography of unglaciated eastern North America. Mol Ecol. 2006; 15: 4261–4293. pmid:17107465

12. Chung MY, López-Pujol J, Chung MG. Is the Baekdudaegan “the Southern Appalachians of the East”? A comparison between these mountain systems, focusing on their role as glacial refugia. Korean J Plant Taxon. 2016; 46: 337–347.

13. Godt MJW, Hamrick JL, Bratton S. Genetic diversity in a threatened wetland species, Helonias bullata (Liliaceae). Conserv Biol. 1995; 9: 596–604.

14. Broyles SB. Postglacial migration and the loss of allozyme variation in northern populations of Asclepias exaltata (Asclepiadaceae). Am J Bot. 1998; 85: 1091–1097. pmid:21684994

15. Bayer RJ, La Duke JC, Crawford DJ. Isozyme variation in Trillium nivale (Liliaceae). Can J Bot. 1987; 65: 2250–2254.

16. Whitkus R, Bryan FA, Les DH, Tyrell LE. Genetic structure in a heterocyanic population of Trillium sessile (Liliaceae). Plant Species Biol. 1987; 2: 67–73.

17. Lewis PO, Crawford DJ. Pleistocene refugium endemics exhibit greater allozyme diversity than widespread congeners in the genus Polygonella (Polygonaceae). Am J Bot. 1995; 82: 141–149.

18. Hannan GL, Orick MW. Isozyme diversity in Iris cristata and the threatened glacial endemic I. lacustris (Iridaceae). Am J Bot. 2000; 87: 293–301. pmid:10718990

19. Petit-Maire N, Bouysse P. Geological records of the recent past, a key to the near future world environments. Episodes. 2000; 23: 230–246.

20. Hope G, Kershaw AP, van der Kaars S, Sun X, Liew P-M, Heusser LE, et al. History of vegetation and habitat change in the Austral-Asian region. Quatern Int. 2004; 118–119: 103–126.

21. Prentice IC, Harrison SP, Bartlein PJ. Global vegetation and terrestrial carbon cycle changes after the last ice age. New Phytol. 2011; 189: 988–998. pmid:21288244

22. Woillez M-N, Kageyama M, Krinner G, de Noblet-Ducoudré N, Viovy N, Mancip M. Impact of CO2 and climate on the Last Glacial Maximum vegetation: results from the ORCHIDEE/IPSL models. Clim Past. 2011; 7: 557–577.

23. Adams JM, Faure H. Review and atlas of palaeovegetation: preliminary land ecosystem maps of the world since the Last Glacial Maximum. Oak Ridge: Oak Ridge National Laboratory. 1997. http://www.esd.ornl.gov/projects/qen/adams1.html

24. Zheng YQ, Qian ZC, He HR, Liu HP, Zeng XM, Yu G. Simulations of water resource environmental changes in China during the last 20,000 years by a regional climate model. Glob Planet Chang. 2007; 55: 284–300.

25. Liang S, Tamura MN. Lilium. In: Wu ZY, Raven PH, Hong DY, editors. Flora of China. Vol. 24 (Liliaceae). Beijing and St. Louis: Science Press and Missouri Botanical Garden Press; 2000. pp. 135–149.

26. Chung MY, Chung MG. Large effective population size and high levels of gene flow between subpopulations of Lilium cernuum (Liliaceae). Biochem Syst Ecol. 2014; 54: 354–361.

27. Chung MY, Chung MG, López-Pujol J, Ren M-X, Zhang Z-Y, Park SJ. Were the main mountain ranges in the Korean peninsula a glacial refugium for plants? Insights from the congeneric pair Lilium cernuum–Lilium amabile. Biochem Syst Ecol. 2014; 53: 36–45.

28. Komarov VL, editor. Flora of the USSR. Volume IV. Liliiflorae and Microspermae (Translated version). Jerusalem: Israel Program for Scientific Translations; 1968.

29. Ku YB, Oh HK, Shim KY, Kim MS, Lee SM. Growth characteristics, genetic diversity and conservation of endangered plants (I): The case of Astilboides tabularis, Euchresta japonica, Echinosophora koreensis and Lilium cernuum. NIER 2006–54–836. Incheon: National Institute of Environmental Research; 2006 [in Korean].

30. Trutnev YP, editor. Red Book of Russian Federation (plants and fungi). Moscow: KMK Scientific Press; 2008 [In Russian].

31. MEP–CAS (Ministry of Environmental Protection–Chinese Academy of Sciences). China Red List of higher plants–Evaluation’s report. Beijing: Ministry of Environmental Protection of the People’s Republic of China & Chinese Academy of Sciences; 2013 [in Chinese].

32. Mitton JB, Linhart YB, Sturgeon KB, Hamrick JL. Allozyme polymorphisms detected in mature needle tissue of ponderosa pine. J Hered. 1979; 70: 86–89.

33. Soltis DE, Haufler CH, Darrow DC, Gastony GJ. Starch gel electrophoresis of ferns: a compilation of grinding buffers, gel and electrode buffers, and staining schedules. Am Fern J. 1983; 73: 9–27.

34. Cheliak WM, Pitel JP. Technique for starch gel electrophoresis of enzyme from forest tree species. Information Report PI-X-42. Chalk River: Petawawa National Forestry Institute; 1984.

35. Haufler CH. Enzyme variability and modes of evolution in Bommeria (Pteridaceae). Syst Bot. 1985; 10: 92–104.

36. Clayton JW, Tretiak DN. Amine citrate buffers for pH control in starch gel electrophoresis. J Fish Res Board Can. 1972; 29: 1169–1172.

37. Yeh FC, Yang RC, Boyle TBJ. POPGENE version 1.31-Microsoft Windows-based freeware for population genetic analysis. Quick users’ guide. Edmonton: University of Alberta; 1999.

38. Goudet J. FSTAT (Version 1.2): A computer program to calculate F-statistics. J Hered. 1995; 86: 485–486.

39. Hurlbert SH. The nonconcept of species diversity: a critique and alternative parameters. Ecology. 1971; 52: 577–586. pmid:28973811

40. El Mousadik A, Petit RJ. High level of genetic differentiation for allelic richness among populations of the argan tree [Argania spinosa (L.) Skeels] endemic to Morocco. Theor Appl Genet. 1996; 92: 832–839. pmid:24166548

41. Nei M. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics. 1978; 89: 583–590. pmid:17248844

42. Piry S, Luikart G, Cornuet JM. BOTTLENECK: A computer program for detecting recent reductions in the effective population size using allele frequency data. J Hered. 1999; 90: 502–503.

43. Nei M, Maruyama T, Chakraborty R. The bottleneck effect and genetic variability in populations. Evolution. 1975; 29: 1–10. pmid:28563291

44. Cornuet JM, Luikart G. Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics. 1996; 144: 2001–2014. pmid:8978083

45. Luikart G, Allendorf FW, Cornuet JM, Sherwin WB. Distortion of allele frequency distributions provides a test for recent population bottlenecks. J Hered. 1998; 89: 238–247. pmid:9656466

46. Hardy OJ, Vekemans X. SPAGeDi: a versatile computer program to analyze spatial genetic structure at the individual or population level. Mol Ecol Notes. 2002; 2: 618–620.

47. Wright S. The interpretation of population structure by F-statistics with special regard to systems of mating. Evolution. 1965; 19: 395–420.

48. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984; 38: 1358–1370. pmid:28563791

49. Peakall R, Smouse PE. GenAlEx 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes. 2006; 6: 288–295.

50. Nei M, Tajima F, Tateno Y. Accuracy of estimated phylogenetic trees from molecular data. II. Gene frequency data. J Mol Evol. 1983; 19: 153–170. pmid:6571220

51. Langella O. Populations version 1.2.30. 1999. http://bioinformatics.org/~tryphon/populations/

52. Page RDM. TreeView: an application to display phylogenetic trees on personal computers. Comput Appl Biosci. 1996; 12: 357–358. pmid:8902363

53. Pritchard JK, Wen X, Falush D. Documentation for structure software: Version 2.3. Chicago: Department of Human Genetics, University of Chicago. 2010. http://pritch.bsd.uchicago.edu/structure_software/release_versions/v2.3.4/structure_doc.pdf

54. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005; 14: 2611–2620. pmid:15969739

55. Earl DA, vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012; 4: 359–361.

56. Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007; 23: 1801–1806. pmid:17485429

57. Rosenberg NA. DISTRUCT: A program for the graphical display of population structure. Mol Ecol Notes. 2004; 4: 137–138.

58. Rousset F. Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics. 1997; 145: 1219–1228. pmid:9093870

59. Legendre P, Lapointe F-J, Casgrain P. Modeling brain evolution from behavior: a permutational regression approach. Evolution. 1994; 48: 1487–1499. pmid:28568410

60. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006; 190: 231–259.

61. Jeong J-H. A study of growth environment and morphological characters of Lilium cernuum Kom. M.Sc. Thesis, Daejin University. 2011.

62. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005; 25: 1965–1978.

63. Brown JL. SDMtoolbox: a python-based GIS toolkit for landscape genetic, biogeographic, and species distribution model analyses. Methods Ecol Evol. 2014; 5: 694–700.

64. Gent PR, Danabasoglu G, Donner LJ, Holland MM, Hunke EC, Jayne SR, et al. The Community Climate System Model Version 4. J Clim. 2011; 24: 4973–4991.

65. Watanabe S, Hajima T, Sudo K, Nagashima T, Takemura T, Okajima H, et al. MIROC-ESM 2010: Model description and basic results of CMIP5-20c3m experiments. Geosci Model Dev. 2011; 4: 845–872.

66. Swets JA. Measuring the accuracy of diagnostic systems. Science. 1988; 240: 1285–1293. pmid:3287615

67. Liu C, Newell G, White M. On the selection of thresholds for predicting species occurrence with presence-only data. Ecol Evol. 2016; 6: 337–348. pmid:26811797

68. Huang P, Schaal BA. Association between the geographic distribution during the last glacial maximum of Asian wild rice, Oryza rufipogon (Poaceae), and its current genetic variation. Am J Bot. 2012; 99: 1866–1874. pmid:23125436

69. Barnard-Kubow KB, Debban CL, Galloway LF. Multiple glacial refugia lead to genetic structuring and the potential for reproductive isolation in a herbaceous plant. Am J Bot. 2015; 102: 1842–1853. pmid:26542847

70. Diniz-Filho JAF, Rodrigues H, Telles MPDC, Oliveira GD, Terribile LC, Soares TN, et al. Correlation between genetic diversity and environmental suitability: taking uncertainty from ecological niche models into account. Mol Ecol Res. 2015; 15: 1059–1066.

71. Merow C, Smith MJ, Silander JA. A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. Ecography. 2013; 36: 1058–1069.

72. Guo X-D, Wang H-F, Bao L, Wang T-M, Bai W-N, Ye J-W, et al. Evolutionary history of a widespread tree species Acer mono in East Asia. Ecol. Evol. 2014; 4: 4332–4345. pmid:25540694

73. Zeng Y-F, Wang W-T, Liao W-J, Wang H-F, Zhang D-Y. Multiple glacial refugia for cool-temperate deciduous trees in northern East Asia: the Mongolian oak as a case study. Mol Ecol. 2015; 24: 5676–5691. pmid:26439083

74. Harrison SP, Yu G, Takahara H, Prentice IC. Diversity of temperate plants in east Asia. Nature. 2001; 413: 129–130. pmid:11557970

75. Wang S, Xu X, Shrestha N, Zimmermann NE, Tang Z, Wang Z. Response of spatial vegetation distribution in China to climate changes since the Last Glacial Maximum (LGM). PLoS ONE. 2017; 12: e0175742. pmid:28426780

76. Cao X, Herzschuh U, Ni J, Zhao Y, Böhmer T. Spatial and temporal distributions of major tree taxa in eastern continental Asia during the last 22,000 years. Holocene. 2015; 25: 79–91.

77. Hamrick JL, Godt MJW. Allozyme diversity in plant species. In: Brown AHD, Clegg MT, Weir BS, editors. Plant population genetics. Breeding and genetic resources. Sunderland: Sinauer Associates; 1989. pp. 43–63.

78. Godt MJW, Hamrick JL. Allozyme diversity in the grasses. In: Cheplick GP, editor. Population biology in grasses. Cambridge: Cambridge University Press; 1998. pp. 11–29.

79. Kim ZS, Hwang JW, Lee SW, Yang C, Gorovoy PG. Genetic variation of Korean pine (Pinus koraiensis Sieb. et Zucc.) at allozyme and RAPD markers in Korea, China and Russia. Silvae Genet. 2005; 54: 235–246.

80. Eckert CG, Samis KE, Lougheed C. Genetic variation across species’ geographical ranges: the central-marginal hypothesis and beyond. Mol Ecol. 2008; 17: 1170–1188. pmid:18302683

81. Hampe A, Petit RJ. Conserving biodiversity under climate change: the rear edge matters. Ecol Lett. 2005; 8: 461–467. pmid:21352449

82. Petit RJ, Aguinagalde I, de Beaulieu J-L, Bittkau C, Brewer S, Cheddadi R, et al. Glacial refugia: hotspots but not melting pots of genetic diversity. Science. 2003; 300: 1563–1565. pmid:12791991

83. Godt MJW, Hamrick JL. Genetic diversity in rare southeastern plants. Natl Areas J. 2001; 21: 61–70.

84. Hedrick PW, Cockerham CC. Partial inbreeding: equilibrium heterozygosity and the heterozygosity paradox. Evolution. 1986; 40: 856–861. pmid:28556171

85. Weir BS. Genetic data analysis II: Methods for discrete population genetic data. Sunderland: Sinauer Associates; 1996.

86. Reed DH, Frankham R. Correlation between fitness and genetic diversity. Conserv Biol. 2003; 17: 230–237.