Contrasting genetic patterns between two coexisting Eleutherococcus species in northern China

Abstract Climate oscillations are the key factors to understand the patterns in modern biodiversity. East Asia harbors the most diverse temperate flora, largely because an extensive terrestrial ice cap was absent during repeated Pleistocene glaciation–interglacial cycles. Comparing the demographic histories of species that are codistributed and are close relatives may provide insight into how the process of climate change influences species ranges. In this study, we compared the spatial genetic structure and demographic histories of two coexisting Eleutherococcus species, Eleutherococcus senticosus and E. sessiliflorus. Both species are distributed in northern China, regions that are generally considered to be sensitive to climatic fluctuations. These regions once hosted temperate forest, but this temperate forest was replaced by tundra and taiga forest during the Last Glacial Maximum (LGM), according to pollen records. Using three chloroplast DNA fragments, we assessed the genetic structure of 20 and 9 natural populations of E. senticosus and E. sessiliflorus, respectively. Extremely contrasting genetic patterns were found between the two species; E. sessiliflorus had little genetic variation, whereas E. senticosus had considerably higher levels of genetic variation (15 haplotypes). We speculated that a recent severe bottleneck may have resulted in the extremely low genetic diversity in E. sessiliflorus. In E. senticosus, populations in Northeast China (NEC) harbored all of the haplotypes found in this species and included private haplotypes. The populations in NEC had higher levels of genetic diversity than did those from North China (NC). Therefore, we suggest that both the NC and NEC regions can sustain LGM refugia and that lineage admixture from multiple refugia took place after the LGM elevated the local genetic diversity in NEC. In NEC, multiple genetic hot spots were found in the Changbai Mountains and the Xiaoxing'an Range, which implied that multiple locations in NEC may sustain LGM refugia, even in the Xiaoxing'an Range.


Introduction
Comparative phylogeography of codistributed species seeks to understand general patterns in the evolutionary history of communities (Bermingham and Moritz 1998;Avise 2000). For example, comparative phylogeography studies have provided many insights on the response of temperate forests to climatic oscillation in North America and Europe (Hewitt 1999(Hewitt , 2000(Hewitt , 2001(Hewitt , 2004Soltis et al. 2006). Many phylogeography studies have been taken in the Qinghai-Tibetan Plateau and subtropical region (Qiu et al. 2011;Liu et al. 2012), while few in temperate forest located in northern East Asia.
Using pollen record data, northern East Asia was most likely replaced with tundra and taiga forest during the Last Glacial Maximum (LGM, approximately 21,000-18,000 years ago), whereas temperate forests probably retreated to the South, below 30°N (Harrison et al. 2001;Qiu et al. 2011;Cao et al. 2015). Interestingly, molecular data have revealed that the genetic diversity of several species in northern East Asia is not low and that genetic divergence existed between populations in southern and northern East Asia. For example, a divergence in the chloroplast genetic diversity of Juglans mandshurica has been found in the southern and northern Yanshan Mountains (Bai et al. 2010). Similar patterns have been recovered in Fraxinus mandschurica (Hu et al. 2008), Ostryopsis davidiana (Tian et al. 2009), Quercus mongolica (Zeng et al. 2010(Zeng et al. , 2011, and Acer mono (Guo et al. 2014;Liu et al. 2014), although the locations of the divergence vary from case to case. All of these cases provide evidence that northern cryptic refugia probably existed widely throughout northern East Asia (Qiu et al. 2011). However, whether all species in northern East Asia were capable of surviving in situ during the LGM is unknown. Not enough cases have been taken in northern East Asia to gain a general dynamic picture of the temperate forest there (Qiu et al. 2011;Liu et al. 2012). Moreover, the majority of these studied species have been trees with distant phylogenetic relationships. What patterns should be expected for congeneric species occurring in the same range is unclear.
Codistributed congeneric species usually have similar life histories, mating systems, and overlapped species ranges, which may be expected a similar level of genetic diversity (Leffler et al. 2012). However, many of those codistributed congeneric species have typically diverged in their niche, probably due to competitive exclusion between closely related species (Gause et al. 1934;Cavender-Bares and Pahlich 2009). These species may have selected different altitudes, soil conditions, mountain slopes, or other aspects (e.g., Maliouchenko et al. 2007). Hence, these codistributed close relatives may react differently to climate change. Comparisons of the demographic histories of these codistributed and closely related species may shed light on how climate change influences species ranges. For example, European white oaks have occasionally followed similar colonization routes during glacial cycles. Movements out of the Iberian and Italian Peninsulas have been clearly identified. The species also survived in separate refugia in the eastern Balkans by colonizing the West, which resulted in the reshuffling of haplotype distributions (Petit et al. 2002). As another example, four forest-dependent bat species, endemic to Taiwan, differed in their tolerance to altitude, the lowland species have experienced genetic divergence due to the mountain barrier before the LGM, whereas high-altitude specialists and altitudinal generalists have had greater levels of gene flow (Kuo et al. 2014).
Eleutherococcus senticosus ( Fig. 1A and B) and E. sessiliflorus ( Fig. 1C and D) have similar ranges and biological characteristics. Both species are temperate deciduous shrubs that are primarily distributed in North China (NC) and Northeast China (NEC) within the Sino-Japanese forest subkingdom (Fang et al. 2009) (Fig. 2A). Both species are pollinated by bees or butterflies (Liu et al. 2002a). E. senticosus is a trioecious species with birdingested seed dispersal (Liu et al. 1997), whereas E. sessiliflorus is hermaphroditic and dichogamous (Liu et al. 2002b). Small differences also exist in the habitat selection of the two species. E. senticosus tends to be distributed at higher altitudes and latitudes. The altitude difference is notable in the southern part of the species range. In the southeastern area of Taihang Mountains, E. senticosus is generally found 300 m higher than E. sessiliflorus. This difference is generally consistent with descriptions from the Flora of China (Shang and Lowry 2007), because E. sessiliflorus is distributed at altitudes between 200 and 1000 m in China, whereas E. senticosus can reach 2000 m. This difference in altitude may reflect the differences in the sensitivity to temperature between the two species. Thus, whether these two species would respond differently to climate change is unknown. Because E. sessiliflorus appears to be more sensitive to temperature decreases, the range of E. sessiliflorus may shrink more severely than that of E. senticosus. E. sessiliflorus and E. senticosus are an ideal species pair to study the influence of climate change on the demographic history of species and on the level of genetic diversity. In this study, we used chloroplast DNA (cpDNA) sequence data and ecological niche modeling (ENM) to compare the demographic histories of the two shrub species. Three specific questions were asked: (1) Do the two codistributed species have similar levels of genetic diversity? (2) What are the differences in the population expansion patterns between the two codistributed species?
(3) Where are the LGM refugia for the two species? Did these species occupy a single refugium or multiple refugia?

Population sampling
Leaf samples of 185 E. senticosus individuals from 20 natural populations and 61 E. sessiliflorus individuals from 9 natural populations were collected from 2012 to 2014 ( Fig. 2A; Table 1). These sampled populations covered the ranges of both species. Within each population, a range of 4-15 individuals at least 30 m apart was sampled. Samples were dried and stored in silica gel until DNA was extracted.

Genetic data analyses
Chloroplast DNA sequence data were edited and aligned using Codoncode Aligner 3.6.1 (Codon Code Corporation, http://www.codoncode.com/aligner/) with the Clustal module. Each indel was treated as a single site in all the following analysis. The genetic diversity parameters (Hd, p) (Nei 1987) and (h) (Watterson 1975) of each population and each species were calculated in DnaSP 5.10.01 software (Rozas et al. 2003). A median-joining  network was constructed by NETWORK 4.6.1.1 (Bandelt et al. 1999). The spatial distribution of haplotypes was mapped in ArcMap 9.3 (ESRI, Inc., RedLands City, CA, USA). To account for differences in sample sizes, allele richness after rarefaction (R S ) and the contributions to total gene diversity of each population (C T ) were calculated with the method described by Petit et al. (1998) in Contrib 1.4 (https://www6.bordeaux-aquitaine.inra.fr/biog eco_eng/Scientific-Production/Computer-software/Contrib -Permut/Contrib). Rarefaction size was set to be the same as the smallest population (e.g., 4, as it was in DS in E. senticosus). The rarefaction analyses partitioned the contribution of C T to C S and C D , with C S representing intrapopulation diversity and C D representing interpopulation diversity. Population differentiation (G ST , N ST ) was also calculated by Permut 1.0 (Pons and Petit 1996), with 1000 permutations. In addition, a Mantel test using Gen-AlEx 6.501 (Peakall and Smouse 2006) was performed to detect the relationship between the genetic and geographic distance of population. The genetic distance was calculated in MEGA 5.0 (Tamura et al. 2011). Potential signals of selection in chloroplast sequences were tested by comparing the relative abundance of the nonsynonymous substitution rate (d N ) and the synonymous substitution rate (d S ) between sequence pairs of the matK gene from each haplotype in the two species. We used the Nei-Gojobori model (Nei and Gojobori 1986) in MEGA 5.0 (Tamura et al. 2011). Z-tests were applied to determine whether d N was equal to d S . When d N equaled d S , we inferred that neutral sequence evolution had occurred. The statistic Z = (d N Àd S )/SQRT(Var(d S ) + Var (d N )) was computed, wherein the variance of the difference was calculated using the bootstrap method (500 replicates).
To assess the spatial genetic structure of populations, values of the fixation index F C T were estimated by spatial analysis of molecular variance (SAMOVA) version 2.0 software (Dupanloup et al. 2002). The highest variance among groups (F C T ) was calculated each time with a user-defined group number (K). We obtained a series of F CT values for group numbers (K) from 2 to 8. The optimal number of groups was determined by the highest F CT value.
To infer the phylogenetic relationships between haplotypes in both species, we used Bayesian inference in Beast 1.8.0 (Drummond and Rambaut 2007). A constant-size coalescent tree prior was used in the analysis with a TIM3+G substitution model determined by jModelTest 2.0.1 (Guindon and Gascuel 2003). A Markov Chain Monte Carlo (MCMC) chain length of 10 7 was used and recorded 10 4 trees. Because the two studied species shared no haplotypes, we treated E. sessiliflorus as an outgroup for the haplotype tree of E. senticosus. The convergence of MCMC simulations was determined by the ESS value in Tracer 1.5 (Rambaut and Drummond 2009). Calculations were not stopped until the ESS value was greater than 200, which was considered to be a sufficient mixture between parameters. The tree topology of the Bayesian tree was visualized in FigTree 1.4.0 (http:// tree.bio.ed.ac.uk/software/figtree/).
Mismatch distribution analysis was conducted in Arlequin 3.5 (Rogers and Harpending 1992;Harpending 1994;Excoffier et al. 2005) to test the hypothesis that spatial population expansion had occurred. The sum of squared deviations (SSD) and the raggedness index (Rag) were calculated, and significance was tested by comparing the observed and simulated mismatch distributions. If the spatial expansion model was not rejected, the expansion time (t) was estimated using the relationship s = 2ut = 2klt, where l is the substitution rate per site per year (s/s/y) and k is the length of the sequence under study. An estimator of the substitution rate for angiosperm species (2.0 9 10 À9 s/s/y) was used in this study (Wolfe et al. 1987). We also calculated Tajima's D in Arlequin 3.5. Tajima's D was used to determine whether rare mutations were significantly higher than expected (Tajima 1989), which is an indicator of recent spatial expansion (Kimura 1983). Because no genetic variation was detected in E. sessiliflorus (see Results), the Mantel test, SAMOVA analysis, and mismatch analysis conducted above were applied only to E. senticosus.

Ecological niche modeling
We used maximum entropy modeling (Maxent 3.3.3k) (Phillips et al. 2006) to predict the distribution of E. senticosus and E. sessiliflorus currently and during the LGM. Occurrence records were obtained from our own field sampling sites and online database, including the Chinese Virtual Herbarium (http://www.cvh.org.cn/), Global Biodiversity Information Facility (http://www.gbif.org/), and Tropicos (http://www.tropicos.org/). To mitigate potential effects of spatial sampling biases, we employed the R package spThin (Aiello-Lammens et al. 2015) to rarefy the data to have a minimum distance of at least 10 km among them. After rarefying, 45 occurrence records of E. senticosus and 22 occurrence records of E. sessiliflorus from unique localities remained for modeling. The current and the LGM climate data were comprised of 19 bioclimatic variables from WorldClim (www.worldclim.org) (Hijmans et al. 2005) at a resolution of 30 arc-sec. Of the 19 climate variables, 7 variables were co-correlated below a 0.80 Pearson |r| correlation and were subsequently used for species distribution modeling. These 7 bioclimatic variables were (1) mean diurnal temperature range, (2) temperature seasonality, (3) maximum temperature of the warmest month, (4) minimum temperature of the coldest month, (5) annual precipitation, (6) precipitation of the driest month, and (7) precipitation seasonality.
For LGM prediction, we used climate data from the Model for Interdisciplinary Research on Climate (MIROC), provided by the Paleoclimate Modelling Intercomparison Project (http://pmip2.lsce.ipsl.fr/). Paleocoast lines and the paleoclimate surfaces of the exposed seafloor area (assuming a sea level depression of 130 m) were estimated using seafloor topography data (ETOPO1) produced by the National Geophysical Data Center of the National Oceanic and Atmospheric Administration (Boulder, CO). We validated our models through choosing a cross-validation set at 50 replicates using the mean area under the curve (AUC) of the receiver-operating characteristics (ROC).

Chloroplast genetic diversity
The three sequenced chloroplast fragments of E. senticosus and E. sessiliflorus were aligned with a consensus length of 2666 bp. Among the sampled individuals of E. senticosus, 17 substitutions and 3 indels were recovered and consisted of 15 haplotypes ( Fig. 2A; Table S1 in the Supporting Information). However, only one haplotype was found among all sampled individuals of E. sessiliflorus (Fig. S1). The only haplotype (H16) of E. sessiliflorus was at least 17 mutational steps from those of E. senticosus (Fig. 2B).
In E. senticosus, all 15 of the haplotypes existed in NEC, and only 5 haplotypes (H1, H3, H5, H6, and H9) were shared in NC. Among the NEC populations, private haplotypes occurred in populations BS (H8 and H10), DS (H11), and CH (H14) in the Changbai Mountains and KDQ (H15) in the Xiaoxing'an Range. Most populations harbored chloroplast polymorphisms, except for 3 populations, including SHB in the Yanshan Mountains, PQG in the Taihang Mountains, and QS in the Changbai Mountains, which were fixed with single haplotypes.
High genetic diversity was found in E. senticosus, with Hd = 0.835, p = 0.00086, and h = 0.00110. The average allele richness within populations (R S ) was 2.413, and the total allele richness (R T ) was 3.148 when the rarefaction size was set to 4. Population BS in the Changbai Mountains had the highest levels of haplotype diversity (R S = 2.367 and Hd = 0.889), with nucleotide diversity p = 0.00124 and h = 0.00146. Population CH in the Changbai Mountains had the most abundant haplotypes (7 haplotypes: H1, H2, H3, H7, H9, H12, and H14), with R S = 2.269, Hd = 0.867, p = 0.00100, and h = 0.00081 (Table 1). Populations QS, CH, and WQ (all in the Changbai Mountains) had the highest contributions to the total gene diversity of the species, with C T values of 0.017, 0.013, and 0.011, respectively (Fig. 3). However, the contributions of the different populations to total gene diversity had different origins. For example, the contributions of CH and WQ were primarily due to withinpopulation gene diversity. In contrast, among-population genetic divergence may account for the high value of C T in QS, which had only one haplotype: H7 (Fig. 3).
The Mantel test for isolation by distance (IBD) failed to detect a significant correlation between the genetic and geographic distance with 999 permutations (r = À0.072, P = 0.239). Z-tests of (d N Àd S ) of each haplotype sequence pairs were not significantly different (e.g., for the sequence pair of H16 and H6, P = 0.533, Table S2).

Population genetic structure and spatial distribution of haplotypes
Extremely low population genetic differentiation was found (G ST = 0.278 and N ST = 0.151). The SAMOVA analysis also indicated a low population genetic structure. Among the predefined number of groups (K = 2-8), F C T values were low: the highest value was 0.261 when K = 2 (Fig. S2), with population QS as one group and the remaining populations as the other group.
The low population genetic structure may be due to the spatial distribution of haplotypes. The most common haplotype, H1, composed 17% of the total samples and 70% of the populations (14 of 20 populations). H5 was present in 30% of the individuals and 65% of the populations. Moreover, other haplotypes were also geographically widespread despite the few individuals, such as H9, with 19 individuals found in 8 populations; H3, with 16 individuals found in 7 populations; and H6, with 28 individuals found in 7 populations ( Fig. 2A; Table 1).
In the rooted phylogenetic consensus tree (Fig. 4), there were two supported clades within E. senticosus.

Ecological niche modeling
The Maxent model performed well for both E. senticosus and E. sessiliflorus, with AUC values of 0.971 AE 0.011 and 0.983 AE 0.015 (mean AE SD), respectively. Under the present climate, the distribution of E. senticosus was predicted to occur in the NEC region, the NC region, and the Korean Peninsula with high probability (Fig. 6A). E. sessiliflorus was predicted to occupy nearly all of the same areas as E. senticosus, except for the northernmost area (e.g., the Xiaoxing'an Range) in NEC (Fig. 6B). The predicted distribution was generally consistent with the actual distribution of each species (Fig. 2A).
Under the LGM climate, the probability of each species being distributed in the northern and eastern margins of NEC decreased, whereas the probability of each species being distributed in southern NEC and northern NC as well as in the Qinling Mountains in the south increased  ( Fig. 6C and D). The difference in the predicted distribution between the two species was relatively minor.

Discussion
Although E. senticosus and E. sessiliflorus are two sympatric species of the genus Eleutherococcus, their genetic patterns were quite different. A high genetic diversity was observed for E. senticosus, with 15 haplotypes detected from 2666 bp cpDNA. In contrast, no genetic diversity was detected in E. sessiliflorus (only one haplotype, H16) (Fig. 2B). It might be normal that different species within the same region were detected different counts of chloroplast haplotypes, such as Acer mono (10 haplotypes from 2140 bp sampled fragments) (Guo et al. 2014), Quercus mongolica (9 haplotypes from 3979 bp) (Zeng et al. 2015), Ostryopsis davidiana (3 haplotypes from 1300 bp) (Tian et al. 2009), and Juglans mandshurica (2 haplotypes from 5850 bp) (Bai et al. 2010). However, some species still had obvious low genetic variation even with longer sampled fragment size, such as Juglans mandshurica and E. sessiliflorus. E. sessiliflorus had even lower genetic diversity than Juglans mandshurica because no genetic variation was detected in distant populations despite sampling more than 11,000 bp of 13 chloroplast loci during the pre-experimental stage (see the Methods section). In addition to the same sequencing size sampled in both Eleutherococcus species, all of the sampled locations of E. sessiliflorus, except for DQG and BSh, overlapped with those of E. senticosus. Both species are congeneric shrubs with similar seeds and pollen dispersal systems, generation times, and species ranges. The mating systems were both outcrossing-dominated and self-compatible (Liu et al. 2002b). The aforementioned factors are crucial in affecting the levels of species genetic variation (Cutter et al. 2013). Besides, with the low sample size, it is possible to miss important genetic variation (61 individuals from E. sessiliflorus, compared to 185 from E. senticosus). However, the 61 individuals from 9 populations of E. sessiliflorus are scattered across the species range. The scattered sample strategy should have characterized the important genetic variation of E. sessiliflorus. We propose that two possible scenarios may have generated the extremely low genetic variation in E. sessiliflorus: (1) a severe bottleneck, which resulted in a low extant effective population size (Charlesworth et al. 2003;Leffler et al. 2012); and (2) selective sweeps, which may have eliminated standing genetic variation (Charlesworth et al. 2003;Frankham 2012). Chloroplast is usually considered as neutral, but one or more trans-species selective sweeps of the chloroplasts have been shown to occur in many Salix species, resulting in an unusual shared pattern of chloroplast haplotypes (Percy et al. 2014). However, our study suggested that the second scenario was less likely in E. sessiliflorus, because Z-tests of (d N Àd S ) in matK were not significant and thus do not permit the rejection of the hypothesis that these sequences have evolved neutrally.
A severe bottleneck may have occurred during the Quaternary climatic oscillation, when harsh environments during the glaciation may have eliminated populations in the northern range. For example, a study by Jakob et al. (2007) has indicated that the lack of cpDNA variation in grass species Hordeum gussoneanum, based on both DNA sequencing and microsatellite length variation, might be due to a severe and relatively recent bottleneck that erased all of the chloroplast variation. Although the total absence of genetic variation throughout the range of each species has not been commonly found in publications, extremely low regional genetic variation, especially in the northern portion of the range of a species, is often observed in temperate species (Palm e and Vendramin 2002;Grivet and Petit 2003;Dorken and Barrett 2004;Magni et al. 2005;Bai et al. 2010Bai et al. , 2015Sakaguchi et al. 2012;Guo et al. 2014). The regional and simultaneous absence of genetic variation in the range of northern species provides evidence that many of these northern populations may have experienced a recent severe bottleneck due to the LGM (Brochmann et al. 2003;Hewitt 2004Hewitt , 2011Soltis et al. 2006). As E. sessiliflorus distributes in lower altitude than E. senticosus, it is possible that the former species is cold-intolerant compared to the latter one. Thus, during the glacial period, E. sessiliflorus may have experienced much more extent of contraction than E. senticosus. However, ecological experiments comparing the cold tolerances of the two species should be taken in the future to gain more solid proofs.
In contrast to E. sessiliflorus, E. senticosus has relatively high genetic diversity. However, low genetic differentiation and no phylogeographic structure were recovered throughout the range of species, because no clear spatial genetic group could be divided based on SAMOVA analysis, and the G ST and/or F CT values were considerably lower than those of other temperate species, which could be higher than 0.8 Ye et al. 2015). In addition, no IBD pattern was found. And the signal of spatial expansion was not explicit, as mismatch analysis indicated an old expansion (0.20 Mya, 0.04-0.59), while Tajima's D revealed no signal of expansion. Nonetheless, all the analyses supported that no expansion started after LGM. Collectively, we infer that multiple LGM refugia may have contributed to the high genetic diversity and low genetic differentiation pattern of E. senticosus.
The genetic diversity in the NEC region was generally higher than that in the NC region. Most populations in NC had only 1-2 haplotypes, with a maximum of 4, whereas populations in NEC generally had 3-7 haplotypes. The average rarified allelic richness in NC (R S = 0.973) was lower than that in NEC (R S = 1.707), although this difference was not significant (P = 0.405). The C T value in NC was negative, except for populations SHB and WLS, whereas 6 populations in the Changbai Mountains had positive values (Fig. 3), which indicated that these populations differentially contributed to overall species genetic diversity. In addition, 4 private haplotypes were detected in the Changbai Mountains, but no haplotypes were detected in the NC region. Two scenarios may explain the findings in the Changbai Mountains. First, the Changbai Mountains may have acted as LGM refugia. Long-term survivorship of refugium populations may result in higher genetic diversity and private haplotypes sustained compared to other newly colonized populations (Hewitt 1996(Hewitt , 2000Petit et al. 2003). The Changbai Mountains have been suggested to be LGM refugia for several studied temperate forest species (Bai et al. 2010Guo et al. 2014;Bao et al. 2015;Zeng et al. 2015). Second, the Changbai Mountains have acted as a recent admixture of multiple LGM refugia. Petit et al. (2003) have indicated that the highest genetic diversity of several European forest species occurs at relatively high latitudes, which has been suggested to be a consequence of lineage admixture from 3 independent southern refugia. For E. senticosus, we propose that both of the above-suggested scenarios may have occurred in the Changbai Mountains.
Although NC had low genetic diversity and no private haplotypes, we still could not rule out the possibility of NC as a LGM refugium, due to its special geographic location. NC has been considered an important refugium and source of postglacial expansion to the NEC region for many studied temperate forest species (Hu et al. 2008;Bai et al. 2010Bai et al. , 2015Guo et al. 2014;Zeng et al. 2015). For E. senticosus, if the Changbai Mountains have LGM refugia (discussed above), NC should have a higher potential to sustain LGM populations because NC is located in the South and should therefore have more suitable habitats (Fig. 6C). The NEC region could receive migrants from NC due to northward postglacial expansion, which could elevate the genetic variation in NEC and lead to genetic variation sharing between NEC and NC. Hence, a relatively high genetic diversity in the NEC region could also be a result of lineage admixture from separate refugia including NC.
In our study, some populations in the Xiaoxing'an Range had relatively high levels of genetic diversity, such as population LS. Moreover, the northernmost population KDQ harbored a private haplotype: H15 ( Fig. 2A). Cryptic refugia of E. senticosus probably existed in the Xiaoxing'an Range, as suggested by other codistributed species (Hu et al. 2008;Bao et al. 2015;Zeng et al. 2015). However, the lack of fossil and geological evidence makes this a weak inference (Qiu et al. 2011;Tzedakis et al. 2013;Cao et al. 2015). As the nuclear genome has a higher genetic diversity than chloroplast, it can provide more resolution in future studies to understand the possibility of northern cryptic refugia.
Although the two codistributed Eleutherococcus species had contrasting genetic patterns, our ENM results did not show big differences between them. ENM suggested that both species could sustain populations in NC, south of NEC, and the Korean Peninsula, as well as the Bo Sea, during the LGM (Fig. 6C and D). In recent years, ENM is generally suggested as a standard analysis for phylogeography studies, because the approach provides independent clues for inferring suitable habitats for species during the glacial period. However, an obvious flaw may make the inference unreliable. Firstly, the current niche modeling approach considers only climate variables and does not account for biological interactions, which can also be critical to species survival in a community (Ricklefs 1987;Hooper et al. 2005). Such data are largely based on ecological surveys, which vary from community to community. However, data of this type are not often collected (Martin-Albarracin et al. 2015;Valladares et al. 2015). Thus, species may fail to survive in the predicted, but climatically suitable habitat. Secondly, coarse resolution of the ENM can smooth the differences between microhabitats. Hence, it is hard to find cryptic refugia as well as microhabitats differences between codistributed species (Dellicour et al. 2014(Dellicour et al. , 2015. Although habitats of the two Eleutherococcus species overlap at a regional scale, the small-scale niche differences (e.g., altitude) are ignored by the current simulations due to the low resolution of ancient climate data. Hence, we include ENM results in our study, but considering the low resolution and significant flaws, we do not think it provides reliable supporting information for understanding the contrasting genetic patterns between the two codistributed Eleutherococcus species.

Conclusions and Implications
Although E. senticosus and E. sessiliflorus are closely related and overlap in their ranges, contrasting genetic patterns were detected based on cpDNA sequencing and coherent sampling strategies of these two congeneric species. E. sessiliflorus had no genetic variation across all of the sampled loci, whereas E. senticosus had considerably higher levels of genetic diversity. We infer that a recent severe bottleneck likely occurred in E. sessiliflorus. For E. senticosus, relatively higher genetic diversity was found in the Changbai Mountains than in NC. We suggest that both NC and the Changbai Mountains may have sustained LGM refugia. Lineage admixture from separate LGM refugia may have contributed to the high genetic variation in the Changbai Mountains. There were signals indicating that the Xiaoxing'an Range may have had cryptic LGM refugia. Confirmation of this possibility will require more evidence. Large amounts of nuclear genetic markers will be needed to accurately infer the origin of the contrasting genetic patterns between these codistributed and closely related species.
The Korean Peninsula and the Bo Sea are other sources of potential LGM refugia for E. senticosus, as predicted by ENM. The Bo Sea was converted to land during the LGM due to the lowering of the sea level. Analysis of samples from the Korean Peninsula would be beneficial for inferring demographic history in future studies. Currently, however, we believe that the unsampled status of the Korean Peninsula and the Bo Sea will not diminish the importance of the conclusions that are made in this paper. Both the Korean Peninsula and Bo Sea refugia may contribute to the extant local gene pool of the Changbai Mountains. Some private haplotypes in the Changbai Mountains may have originated from the two above unsampled LGM refugia. However, no IBD pattern was recovered. Populations with high genetic diversity were located at medium to high latitudes (such as WQ, CH, and LS). This location is considerably further north than would have been possible to admix in the region containing the three potential southern refugia (NC, the Bo Sea, and the Korean Peninsula, Fig. 2). Hence, the Changbai Mountains may have sustained separate refugia during the LGM.

Supporting Information
Additional Supporting Information may be found online in the supporting information tab for this article: Figure S1. Geographic distribution of chloroplast haplotype in Eleutherococcus sessiliflorus. Figure S2. SAMOVA of E. senticosus based on cpDNA. Table S1. Designation of the cpDNA haplotypes detected at the 3 chloroplast loci among 20 E. senticosus populations and 9 E. sessiliflorus populations Table S2. Z-tests of (d N -d S ) of each haplotype sequence pairs