Rear‐edge, low‐diversity, and haplotypic uniformity in cold‐adapted Bupleurum euphorbioides interglacial refugia populations

Abstract The high genetic diversity of rear‐edge refugia populations is predicted to have resulted from species repeatedly migrating to low latitudes during glacial periods over the course of Quaternary climate change. However, several recent empirical studies of cold‐tolerant plants revealed the opposite pattern. We investigated whether current habitats of the cold‐adapted and range‐restricted Bupleurum euphorbioides in the Baekdudaegan, South Korea, and North Korea could be interglacial refugia, and documented how their rear‐edge populations differ genetically from those of typical temperate species. Phylogeographic analysis and ecological niche modeling (ENM) were used. The genetic structure was analyzed using microsatellite markers and chloroplast DNA sequences. The congener B. longiradiatum was included as a typical temperate plant species. Despite having almost identical life history traits, these congeneric species exhibited contrasting patterns of genetic diversity. ENM revealed an apparent maximum range contraction during the last interglacial. In contrast, its range expanded northward to the Russian Far East (Primorsky) during the Last Glacial Maximum. Thus, we hypothesize that B. euphorbioides retreated to its current refugia during interglacial periods. Unlike populations in the central region, the rear‐edge populations were genetically impoverished and uniform, both within populations and in pooled regional populations. The rear‐edge B. euphorbioides survived at least one past interglacial, contributing to the species’ genetic diversity. We believe that such genetic variation in the cold‐adapted B. euphorbioides gives the species the necessary adaptations to survive an upcoming favorable environment (the next glacial), unless there is artificial environmental change.


| INTRODUC TI ON
Genetic structure in extant plant populations is known to be affected by various factors, including the diversity of pollinators, reproductive systems, seed dispersal modes, and historical migration patterns; historical range change during Quaternary climatic oscillations is considered a primary factor (Petit et al., 2003). Although their relative importance may vary across time and space, the genetic features of populations in northeastern Asian temperate regions likely reflect historical, rather than current, levels of gene flow (Aizawa, Kim, & Yoshimaru, 2012;Bao et al., 2015;Chen et al., 2012;Tamaki et al., 2019). Quaternary glacial periods have repeatedly forced temperate plant species to retreat southward, where they remained for long periods (Hewitt, 1996(Hewitt, , 2000. Rear-edge populations are typically described as low-latitude glacial refugia, providing long-term stores of genetic diversity with high levels of variation and uniqueness (Provan & Maggs, 2011). Such genetic determinants, acquired via repeated adaptation processes, may play a key role in guaranteeing long-term survival (Gugger, González-Rodríguez, Rodríguez-Correa, Sugita, & Cavender-Bares, 2011;Lepais et al., 2013). Therefore, understanding rear-edge population genetic diversity is likely to become a conservation research priority; it will provide a framework for assessing the ability of natural populations to adapt to changing environmental conditions (Provan & Maggs, 2011;Scalfi, Piotti, Rossi, & Piovani, 2009).
Many studies on glacial refugia have considered the genetic diversity of temperate, usually warm-adapted, species (Petit, Hu, & Dick, 2008). This may be because there are relatively few coldadapted temperate plant species. Furthermore, the term "glacial refugia" may be too broad, as it does not adequately convey the wide range of behavioral patterns. Quaternary refugia are regions that species inhabit during the maximum range contraction in the glacial/ interglacial cycle (Stewart, Lister, Barnes, & Dalén, 2010). Several empirical studies in East Asia have explained reasonably well how this occurs, accounting for the time frame of the maximum range contraction during the Quaternary (e.g., for Rosa sericea: Gao, Zhang, Gao, & Zhu, 2015;Chrysanthemum indicum: Li, Wan, Guo, Abbott, & Rao, 2014;and Pinus kwangtungensis: Tian, López-Pujol, Wang, Ge, & Zhang, 2010). Nonetheless, our understanding of the genetic history and structure of rear-edge populations occurring in cryptic southern refugia remains to be elucidated.
Baekdudaegan, a main mountain range of the Korean Peninsula, is a key ecological axis of northeastern Asia; it is connected to the Sikhote-Alin mountains in the Russian Far East and the Lesser Khingan Range in China, spanning ca. 1,600 km from South Korea.
It is one of the longest mountain chains in Asia, with 14 subsidiary mountain ranges; it is home to many boreal and temperate species, as well as to widespread and range-restricted species (Chung, López-Pujol, & Chung, 2016). The range is also a well-known hotspot for biodiversity, particularly for threatened species (Chung et al., 2018).
Although most of it remains relatively pristine, some parts of Baekdudaegan in South Korea are threatened. Its high conservation value has made it the site of many genetic studies, which have provided valuable information for conserving threatened species Chang, Kim, Park, & Maunder, 2004;Jeong et al., 2010). Some studies have focused on the history of population establishment (Chung, Chung, López-Pujol, Park, & Chung, 2013;Chung et al., 2012). As an important refugium, the parts of Baekdudaegan in South Korea are thought to have provided a stable habitat for the boreal and temperate species of northeastern Asia during glacial periods (Chung, López-Pujol, & Chung, 2017).
Bupleurum euphorbioides Nakai (Apiaceae) is a cold-adapted annual or biennial herbal plant that is distributed on rocky slopes at the timberline (ca. 1,500 m above sea level); it is restricted to high-elevation peaks along the main Baekdudaegan ridge on the Korean Peninsula, both in South and North Korea (Figures 1 and 2).
Although the plants also occur further north than Baekdudaeganmostly around Mt. Baekdu in northeastern China and Vladivostok in the Russian Far East (notably, it occurs at lower altitudes in these regions)-most of its range occurs in Baekdudaegan. Its habitat is severely fragmented and threatened by human disturbance, such as increasing mountain tourism in South Korea, and by overgrazing by wild animals (Kim & Chang, 2005). Therefore, B. euphorbioides is internationally considered to be endangered ("EN" status) based on the F I G U R E 1 Photographs of study species taken by Dr. Soonku So on Mt. Seorak, South Korea. (a) Habitat on a rocky slope; (b) habitat under a temperate tree IUCN Red List of Threatened Species (Kim, Kim, & Son, 2016), and it is locally managed as vulnerable ("VU" status) in North Korea (Ju et al., 2016). In contrast, B. longiradiatum Turcz., a congener, is widespread at mid-elevations along the Baekdudaegan range and has a large distribution spanning from Mongolia to northeastern China, the Russian Far East, and Japan (Ohba, 1999;Sheh & Watson, 2005;Shishkin, 1950).
Plant genetic studies are often conducted by comparing congeneric species (Chung et al., 2013;Kim & Chang, 2005). This makes it possible to interpret results in an evolutionary framework.
Comparisons of congeneric species allow for partial control of "background variation" caused by differences in life history traits, particularly breeding systems and seed dispersal mechanisms (Chung et al., 2012). We speculate that the distribution of B. euphorbioides, which might have been continuous and extensive during glacial periods, now comprises isolated fragments as a result of postglacial vertical migration; similar patterns have been observed for other cold-tolerant plants (Gao et al., 2015;Tian et al., 2010). We therefore compared B. euphorbioides to B. longiradiatum, which has undergone range contraction during glacial periods (Zhao, Wang, Ma, Liang, & He, 2013).
In this study, we examined the possibility that the current habitats of B. euphorbioides are interglacial refugia, paying particular attention to the rear-edge populations. The genetic structures were analyzed using previously developed microsatellite markers and chloroplast DNA sequences. Ecological Niche Modeling (ENM) was employed to determine the maximum contraction stage of B. euphorbioides' geographic range. Furthermore, we characterized the genetic characteristics of the rear-edge populations to assess their contribution to the species' genetic diversity. Finally, we provide conservation guidelines for the recovery and management of threatened B. euphorbioides populations in the boreal and temperate climatic zones of the Baekdudaegan range.  1-8 km apart and are fragmented and isolated from one another. We regarded isolated populations between which seed dispersal could not occur as one independent populations. We were unable to find any B. euphorbioides individuals at a location on Mt. Odaesan where we had previously located and monitored them; that population is presumed to be locally extinct. To be able to compare patterns of genetic variation in B. euphorbioides and B. longiradiatum, we also collected 49 B. longiradiatum leaf samples (mean: 6.83) from the six populations in South Korea.

| DNA extraction, genotyping, and sequencing
Total genomic DNA was extracted from the samples using MG Plant Genomic DNA Extraction SV Miniprep Kit (Macrogen Inc., Seoul, Republic of Korea) according to the manufacturer's protocol. The extracted DNA was confirmed by electrophoresis, and the concentration was measured using a Nano-300 micro spectrophotometer, to obtain the same concentration of template DNA in each sample; this was then diluted to 15 ng of DNA for analysis.
For our preliminary cpDNA analysis, we sequenced three chloroplast non-coding regions from two individuals per population (a total of 16 individuals in the eight populations), using universal makers widely used in analyzing Bupleurum phylogeny (Ma, Zhao, Wang, Liang, & He, 2015). The preliminary analysis showed that trnL-trnF and rpl32-trnL, but not the rps16 region, were unsuitable because they contained PolyT or PolyA sequences, or had low levels of intraspecific genetic differentiation. Thus, we designed two primer pairs, after selecting regions that exhibit polymorphism within the species, using the cpDNA genomic data of B. euphorbioides and Taq polymerase, and sufficient distilled water. Amplification was performed as follows: initial denaturation at 95°C for 10 min; 35 cycles of denaturation at 94°C for 1 min, annealing at 52°C for 1 min, and extension at 72°C for 1 min; and a final extension at 72°C for 7 min. The PCR products were visualized on 1% agarose gels dyed with Redsafe (120 V, 25 min), purified using a Wizard SV Gel and PCR Clean-UP system (Promega), and sequenced using an ABI 3730xl DNA Analyzer (Applied Biosystems). The sequences of the three regions were aligned using the Clustal algorithm (Thompson, Higgins, & Gibson, 1994). The three plastid fragments were then combined and manually adjusted using Geneious v. 10.2.3 software. Based on the concatenated sequence alignments, the cpDNA haplotypes were determined. The sequences were deposited in the GenBank database (accession numbers MT188515-MT188552).

| Data analysis
Before analyzing the microsatellite data, two individuals identified as homozygous were randomly selected and sequenced by PCR amplification for each locus to account for the presence of null alleles and confirm the accuracy of the markers. The results of the sequencing showed that all markers perfectly amplified their target loci for both samples. We also estimated the null allele frequency using INEst (Inbreeding/Null allele Estimation) software, which calculates the null allele frequency regardless of the effect of inbreeding (Chybicki & Burczyk, 2009). This analysis showed that the null allele frequency of all the loci examined ranged from 1.5% to 6.7%. Therefore, we Population genetics summary statistics for each species were computed using GenAlEx 6.5 (Peakall & Smouse, 2006). These included the number of alleles, mean expected heterozygosity (H E ), mean observed heterozygosity (H O ), number of private alleles (P A ), and fixation index (F IS ). Allele richness (A R ) and genetic differentiation among populations (F ST ) were determined by calculating the overall F IS according to the method of Weir and Cockerham (1984), using FSTAT 1.2 (Goudet, 1995). The statistical significance of F ST was tested using the log-likelihood (G)-based exact test in FSTAT. To test for departures from Hardy-Weinberg equilibrium (HWE) and linkage equilibrium, we conducted exact tests based on a Markov Chain method (1,000 permutations), using GENEPOP 4.0 (Rousset, 2008). The possibility of recent reductions in effective population size was detected using BOTTLENECK 1.2.02 (Cornuet & Luikart, 1996) (10,000 iterations). We utilized two models for evolution-a two-phase model (  the Bayesian Wilcoxon signed-rank test, to evaluate departures from a 1:1 deficiency/excess ratio (Luikart, Allendorf, Cornuet, & Sherwin, 1998).
To infer the population structure of each species, we used a Bayesian clustering approach implemented in STRUCTURE 2.3, as calculated from microsatellite markers (Pritchard, Stephens, & Donnelly, 2000), using 1,000,000 Markov Chain Monte Carlo (MCMC) iterations (100,000 burn-in, with admixture). The simulation used 20 iterations, with K = 1 to K = 8 clusters for B. euphorbioides, and K = 1 to 6 for B. longiradiatum. The optimal number of clusters, K, was found via the K method, using STRUCTURE HARVESTER (Earl & vonHoldt, 2012). CLUMPP v. 1.1.2 (Jakobsson & Rosenberg, 2007) with the Greedy algorithm was used to combine the membership coefficient matrices (Q-matrices) from 1,000 iterations for K = 2 and K = 3, using random input orders.

| Genetic statistical data
Genetic diversity parameters are shown in Table 1

| Genetic diversity and differentiation in B. euphorbioides
Relative euphorbioides, F IS was high, with a mean of 0.192 (Table 1). At the regional level, for B. euphorbioides, the within-popula- Similarly, genetic diversity was higher for the pooled central populations (Table 2). Consistent with these findings, we detected the signatures of recent bottlenecks in the following populations: e-SD (under TPM and mode shift models of evolution), e-SN (mode shift), and e-DN (mode shift); these are all rear-edge populations (

| Population genetic structure
For B. euphorbioides, the STRUCTURE analysis showed that the mean likelihood scores increased, with K values of 1-8, whereas ΔK decreased gradually, with K values of 2-7 ( Figure 3). However, the standard deviation (SD) of L(K) increased considerably from K = 4 ( Figure 3a). Therefore, we chose two as the optimal number of clusters, with three being the next most optimal (Figure 3b)

| Ecological niche modeling
The ENM of B. euphorbioides ( Figure 6) had a high average AUC (0.967), supporting its predictive power. The most important variable was elevation (49.8% contribution), followed by bio_05 (maximum temperature of warmest month; 28.5%) and bio_04 (temperature seasonality; 7.7%). The estimated LGM distribution extended farther north-east than did the LIG distribution, and it expanded into

| Interglacial endurance of Bupleurum euphorbioides
We found that the levels of genetic variability between these congeneric species were conspicuously different. Within-population genetic variation was low or absent in the range-restricted B. euphorbioides, based on both the microsatellite and cpDNA analyses. In contrast, it was high in the widespread B. longiradiatum, which showed low genetic differentiation among populations. Given that these species have almost identical life histories (Sheh & Watson, 2005), these contrasting patterns of genetic diversity suggest that historical factors are the most relevant factors. Therefore, we suggest that these species have different evolutionary histories. Previous findings show that the distribution of B. longiradiatum was fragmented during the last glacial period into two LGM refugia, in the Qinling and Changbai Mountains in China, and expanded during the interglacial period from these refugia (Zhao et al., 2013). Although B. longiradiatum did not undergo extensive latitudinal range shifts, its shift was similar to those of other East Asian temperate plant species that experienced range contractions during the Quaternary glacial periods (Bao et al., 2015;Chen et al., 2012). Our results also reveal that, in South Korea, B. longiradiatum has relatively high genetic diversity with geographic admixing of lineages, implying that these populations were founded by postglacial immigration from multiple glacial refugia (Petit et al., 2003), presumably in China and Japan.
Alternatively, the rear-edge of Baekdudaegan may have acted as a "glacial" refugium for B. longiradiatum, as it has been shown to do for other temperate plants (Chung et al., 2017).
The ENM for B. euphorbioides revealed a maximum range contraction to only the extremely limited southern part of the Baekdudaegan mountain range during the last interglacial (Figure 6c). By compari-

| Rear-edge Baekdudaegan mountain peaks as cryptic southern refugia
The genetic distinction between the central and rear-edge B. euphorbioides populations was revealed in both the microsatellite and cpDNA analyses ( Figure 2). However, genetic structuring within each region seems to have progressed quite differently. There were no shared haplotypes among the central populations, but they had a uniform gene pool according to the STRUCTURE analysis, except for population e-SB, which was distinct according to the microsatellite analysis. In contrast, the rear-edge populations contained one fixed haplotype (E1), although the microsatellite analysis revealed slight genetic admixture.
Likewise, the two analyses produced opposite patterns of genetic differentiation when comparing B. euphorbioides populations within each region. Such nuclear-chloroplast DNA discordance within regions in this species might be explained by the differences in the geographic distances of pollen dispersal (fly-pollinated; Kim & Chang, 2005) and seed dispersal (gravity-dispersed; Bonet & Pausas, 2004;Zając, Zając, & Tokarska-Guzik, 2009). Although flies have a limited flight distance compared to other pollinators (Kim & Chang, 2005), the fly species in these highland regions are probably capable of flying several kilometers (Inouye, Larson, Ssymank, & Kevan, 2015). Given this species' isolation on mountain peaks, the geographic gaps present significant barriers for F I G U R E 6 Potential distributions of Bupleurum euphorbioides during (a) the present, (b) the Last Glacial Maximum, LGM, and (c) the last interglacial, LIG. Distributions predicted by ecological niche modeling; potential distribution during the LGM was averaged from three general circulation models seeds, no matter how close the populations are. For pollen dispersal, however, the distance of ca. 5 km between the central populations may not pose a substantial barrier. Therefore, in the central region, it appears that the genetic traces of chloroplasts are well preserved, whereas the nuclear information seems to have undergone mixing as a result of current gene flow; this is reflected in the positive correlation between genetic differentiation and geographic distance (R 2 = .317, p = .020, Figure 4a). By comparison, rear-edge populations appear to reflect genetic traces caused by historical colonization events well. In general, the smaller the population size, the stronger the random genetic drift is, and this promotes gene fixation. Moreover, the rear-edge contained only one haplotype. These findings contradict the widely accepted hypothesis that genetic diversity decreases as latitude increases (Hewitt, 2000), presenting evidence that at least some glacial southern refugia retain high genetic diversity in pooled regional populations (as proposed by Diekmann & Serrao, 2012). Interestingly, the fixation of a single haplotype at the rear-edge of the B. euphorbioides distribution resembles the genetic footprints of warm-adapted species on scattered islands, and in northern peripheral populations founded by postglacial northward migration from putative southern glacial refugia (Lee, Lee, & Choi, 2013).
Our study is important and novel for several reasons. First, we clearly described the genetic characteristics of rear-edge populations while inferring interglacial refugia in a phylogeographical framework.
This particular case addresses fundamental aspects of species' survival, along with the function of interglacial cryptic southern refugia, and the mechanisms whereby they maintain genetic diversity.
Second, we show that the Baekdudaegan mountain system, an important refugium for many temperate plants during glacial periods, can also play this role in interglacial periods. Low genetic diversity in populations is typically expected to negatively impact species' reproductive fitness and evolutionary potential and, ultimately, lead to extinction (Dyke, 2008;Spielman, Brook, & Frankham, 2004). However, we believe that the genetic impoverishment, and even haplotypic uniformity, of the rear-edge populations of the cold-adapted B. euphorbioides represents strategic selection that will help it to survive the upcoming favorable environment (the next glacial). In light of the fact that all central populations have private alleles, it is highly likely that the central region also represents interglacial refugia. To determine this, future studies should test leading edge samples collected around the northern end of the Baekdudaegan range.

| Implications for conservation
From a long-term conservation genetics perspective, it is especially important for B. euphorbioides that conservation efforts should not be focused on a limited set of populations or habitats. All extant habitats, which represent interglacial refugia, have the evolutionary potential to guarantee the long-term survival of the species. The low genetic diversity of the rear-edge populations may in fact be beneficial to the species' survival. Given that the B. euphorbioides populations we sampled here represent the species' entire South Korean distribution, it should not be difficult to conserve all of these populations in situ. Thus, first, we recommend that all known populations be protected by law to prevent further damage by loss of individuals. An appropriate measure is to control the access roads by which people can enter (e.g., Norihisa & Suzuki, 2006;Willard, Cooper, & Forbes, 2007).
Based on the genetic distinctions, we observed between the central and rear-edge populations, a conservation strategy should be developed to maintain the unique genetic identity of each region. In terms of contemporary genetics, conservation efforts should focus on retaining the genetic uniformity of the rear-edge rather than increasing the genetic diversity. If artificial restoration of habitats is required, it is preferable to select a single donor population after considering its lineage. Ex situ conservation efforts may be effective in the central region (especially on Mt. Seorak) where the haplotypes are not shared among populations. Finally, the fact that pollen transfer could enable gene flow over several kilometers should be considered in restoration attempts.