Genomic Data Reveals Profound Genetic Structure and Multiple Glacial Refugia in Lonicera oblata (Caprifoliaceae), a Threatened Montane Shrub Endemic to North China.

Characterizing genetic diversity and structure and identifying conservation units are both crucial for the conservation and management of threatened species. The development of high-throughput sequencing technology provides exciting opportunities for conservation genetics. Here, we employed the powerful SuperGBS method to identify 33, 758 high-quality single-nucleotide polymorphisms (SNP) from 134 individuals of a critically endangered montane shrub endemic to North China, Lonicera oblata. A low level of genetic diversity and a high degree of genetic differentiation among populations were observed based on the SNP data. Both principal component and phylogenetic analyses detected seven clusters, which correspond exactly to the seven geographic populations. Under the optimal K = 7, Admixture suggested the combination of the two small and geographically neighboring populations in the Taihang Mountains, Dongling Mountains, and Lijiazhuang, while the division of the big population of Jiankou Great Wall in the Yan Mountains into two clusters. High population genetic diversity and a large number of private alleles were detected in the four large populations, while low diversity and non-private alleles were observed for the remaining three small populations, implying the importance of these large populations as conservation units in priority. Demographic history inference suggested two drastic contractions of population size events that occurred after the Middle Pleistocene Transition and the Last Glacial Maximum, respectively. Combining our previous ecological niche modeling results with the present genomic data, there was a possible presence of glacial refugia in the Taihang and Yan Mountains, North China. This study provides valuable data for the conservation and management of L. oblata and broadens the understanding of the high biodiversity in the Taihang and Yan Mountains.

Introduction

The Earth’s biodiversity consists of approximately 9 million types of organisms (e.g., plants, animals, protists, and fungi), while biodiversity losses have exerted a profound impact on the ecology and society as a whole (Chapin et al., 2000; Cardinale et al., 2012). A total of 10% of tree species (> 8,000) on the earth are threatened with extinction (Cavender et al., 2015). Endemic species are important members of biodiversity hotspots and are of great value to biodiversity conservation (Myers et al., 2000). They are often characterized by narrow geographic ranges and specialized niche requirements, and they tend to have a small population size (Işik and Kani, 2011). Due to the sensitivity to internal factors (e.g., genetic bottleneck, inbreeding depression, and genetic drift) and/or external factors (e.g., human disturbance and environmental randomness), species with small populations are more vulnerable, which may lead to reduced fitness of certain individuals (Charlesworth and Willis, 2009), or even extinction. Moreover, low genetic variation and gene flow are frequently detected in rare plants compared to common plants (Cole, 2003). Understanding the genetic structure of endangered woody species can provide valuable information for their conservation and management.

The Sino-Japanese floristic region holds one of the oldest floras in the North Hemisphere with high species richness (Chen et al., 2018b; Lu et al., 2018). Climatic oscillations and geological events have greatly influenced the genetic pattern and distributional range of many plant species, particularly during the Quaternary glacial–interglacial cycles. Compared to North America and Europe, the fauna and flora in Asia were reported to be less affected by glaciation events. Northern China was less affected by massive ice sheets (Qiu et al., 2011), and an arid belt was suggested to have occurred in this region (Guo et al., 2008). Composing the main part of the Sino-Japanese floristic region, the northern China flora exhibits unique genetic patterns of forestry species. Southern warm-temperate forests mixed with northern cool-temperate forests occur in this region, indicating complicated genetic diversity and genetic differentiation patterns (Ye et al., 2017). It is commonly reported based on limited molecular markers that during the northward expansion of species, there would be a gradual decline in population genetic diversity [e.g., Acer mono Maxim. (Guo et al., 2014; Liu et al., 2014) and Quercus mongolia Fisch. ex Ledeb. (Zeng et al., 2015)]. In addition, single/multiple refugia in northeast China and high genetic diversity in northern populations are also observed [e.g., Juglans mandshurica Maxim. (Bai et al., 2010), Eleutherococcus senticosus (Rupr. and Maxim.) Maxim. (Wang et al., 2016), and Schisandra sinensis C. Bailley (Ye et al., 2019)]. North China is an important geographic and floristic region of northern China; it contains both numerous rugged mountains and the North China Basin and provides a north–south migration corridor for wildlife. Compared to the substantial progress made by previous literature focusing on the species with a broad geographic range, the genetic patterns of woody plants endemic to North China remain to be clearly understood, especially in the era of population genomics.

Two great mountains are located in North China, namely, the NNE trending Taihang Mountains (THMs) and the EW trending Yan Mountains (YMs), which harbor the highest seed plant diversity in northern China, as well as a high endemism rate (Wang et al., 1995). The THMs lie between the Ordos-Shanxi plateau on the west and the North China Basin on the east, with a north latitude from 34° to 40° (Wang and Li, 2008). The THMs are also a transitional region located between the second staircase (altitude 1,000–2,000 m) and the third staircase (altitude < 500 m) in the three-step topography of China. The intense uplift during the Late Pliocene to Pleistocene (Wu et al., 1999), along with heterogeneous local geographic and environmental changes (Guo et al., 2008; Wang et al., 2018), provided a unique habitat for several species endemic to this region, such as Opisthopappus taihangensis (Ling) Shih and Taihangia rupestris Yuet Li. The YMs lie in the north of the North China Basin with an east longitude from 115° to 119°. The west end of the YMs is connected to the northeast end of THMs. Most mountains in the THMs and YMs are extremely steep, particularly in the THMs. These two great mountains meet in the northwest of Beijing, and the valley between Changping District on the westside and Yanqing District on the northside are considered as the division. Uncovering the genetic composition of species endemic to the THMs and YMs may help us understand the genetic pattern and preservation of biodiversity in North China and can also allow us to further explore the role of these mountains in the survival and dispersal of endemic species.

High-throughput sequencing (HTS) technology can rapidly detect thousands of single-nucleotide polymorphisms (SNPs) in a cost- and time-efficient manner on a genome-wide scale. Compared to Sanger sequencing method, which generates limited loci, the large data generated from HTS can provide robust phylogeny, reducing the potential incomplete lineage sorting, and enable us to address a wide range of evolutionary questions. Furthermore, a small sample size per population (e.g., a minimum of two individuals) can also generate big data in the population genomic study based on the HTS approach (Nazareno et al., 2017). Thus, it is an effective tool for the assessment of population genetic diversity compared to traditional genetic markers (Glover et al., 2010; Li et al., 2020). Restriction site-associated DNA (RAD) is one of the most effective genotyping-by-sequencing (GBS) methods that allow for extensive SNP discovery on the genome level (Davey et al., 2011; Sonah et al., 2013; Puritz et al., 2014). It has been widely applied in population genetics for numerous model and non-model species (e.g., Wang et al., 2013; Luo et al., 2019; Feng et al., 2020; Xiong et al., 2020; Boukteb et al., 2021; Qiao et al., 2021). In particular, the UGbS-Flex (also denoted as SuperGBS) is a powerful tool for population genetic research with the ability to generate long reads (300–700 bp) and maximize SNP callings irrespective of the species ploidy level, breeding system, and reference genome availability (Qi et al., 2018; Cheng et al., 2020).

Lonicera oblata K. S. Hao ex P. S. Hsu and H. J. Wang (Caprifoliaceae) is a critically endangered montane shrub endemic to North China (Zhu et al., 2019; Wu et al., 2021), and it is ranked second in the list of national key protected wild plants in China.[1] Results from our 13-year field investigation identified seven highly fragmented populations that typically grow in the THMs and YMs (Figure 1). The majority of individuals occurred in habitats of stony and steep cliffs in open forests with an elevation of ca. 1,100 m. Although the red berries of L. oblata are frequently observed in the wild, no seedlings are reported. The habitats of several populations are exposed to human activities (e.g., tourism, logging, and grazing), and some individuals are injured directly by cutoff. The highly threatened condition of L. oblata may be attributed to the small number and size of populations scattered in such a vast geographic range, unique but extremely fragmented habitats, strong human disturbances, and sensitivity to climate change (Wu et al., 2021). The highly isolated distributional pattern and mountain-top preferred habitat of L. oblata is similar to species in the sky islands in southwest China (He and Jiang, 2014), which may suffer both morphological and genetic variation because of local geographic vicariance and population isolation. With the exception of its endangerment status, little is known about the biological and genetic characteristics of L. oblata. Thus, research on the genetic diversity and population structure of L. oblata at the genome level is imperative for the adequate determination of conservation strategies for policy-makers.

FIGURE 1

Geographic distributions of sampled Lonicera oblata populations. The two red circles on the northside represent populations from Yan Mountains, and the five blue circles represent populations from Taihang Mountains. The circle size corresponds to the population sample size.

In this study, we employed the novel approach of RAD sequencing, SuperGBS (Qi et al., 2018), on the threatened montane shrub L. oblata endemic to the two great mountains, the THMs and the YMs, in North China to (i) assess the population structure and genetic diversity across the full range of the species; (ii) infer the genetic signatures of potential refugia during the Quaternary climate oscillations in the THMs and YMs; and (iii) identify conservation units (CU) and provide suggestions for the conservation of this species. Our results provide both a theoretical basis and a valuable genetic database for the conservation planning of L. oblata and contribute to improving our understanding of the preservation and evolution of biodiversity in North China.

Materials and Methods

Sample Collection and DNA Extraction

A total of 134 accessions from all seven natural populations covering the entire geographic range of L. oblata were sampled. Based on the number of individuals, the seven populations can be classified into four large populations [Heduling (HDL, N = 24), Wutai Mountains (WTS, N = 34), Jiming Mountains (JMS, N = 20), and Jiankou Great Wall (JK, N = 36)] located at the species marginal distribution region and three small populations [Bijia Mountains (BJS, N = 6), Dongling Mountains (DLS, N = 6), and Lijiazhuang (LJZ, N = 8)] located at the central region (Figure 1). The only individual found in Beijing Songshan National Nature Reserve, which is near the JK population in Beijing, was not included in this study. Among them, JK and JMS are distributed in the YMs, while the remaining populations are located in the THMs (Figure 1 and Supplementary Table 1). A minimum distance of 30 m from one individual to another was set during the sampling. Leaves were dried in silica gel in the field and stored at −80°C until further use. Voucher specimens were preserved in the herbarium at Beijing Forestry University. Genomic DNA was extracted from the above leaf tissue using the DNAsecure Plant Kit (Tiangen Biotech, Beijing, China) following the manufacturer’s protocol. DNA quantity and quality were determined by Agarose Gel and NanoDrop.

Preparation of Libraries and Sequencing

Genotyping-by-sequencing was performed following the method described by Qi et al. (2018) using the enzyme combination of PstI/MspI. Briefly, 250 ng of DNA from each sample was double-digested with PstI and MspI, followed by the barcode adapter ligation at the PstI site and a common Y-adapter at the MspI site. Unligated adapters were removed by the recovery system of the improved magnetic bead. Recovered fragments with a length of 300–700 bp were PCR-amplified and tested for density using Qubit2.0 to ensure density values greater than 5 ng/μl. The libraries were subsequently sequenced using the Illumina HiSeq X Ten, PE150 platform at OE Biotech Co., Ltd., Qingdao, China.

Single-Nucleotide Polymorphisms Calling

Raw reads were split by barcode using Stacks version 2.4 (Catchen et al., 2013) with the option “process_radtags –renz-1 -r -s 0 –retain_header –adapter_mm 1 –adapter_1 –adapter_2 –b.” The quality of passed reads was checked and filtered using Fastp version 0.20.0 (Chen et al., 2018a). Restriction sites, base quality values < 20, and the last 5 bp of the raw reads were likely to contain errors and were thus removed to obtain clean reads with the option “–n_base_limit 5 –cut_window_size 4 –cut_mean_quality 20 –length_required 75 –qualified_quality_phred 15.” As the whole genome information of L. oblata was not yet reported at the time of this study, a de novo assemble of the GBS reference genome was generated in Stacks version 2.4 with the option “ustacks -M 3, cstacks -n 3, sstacks, tsv2bam, gstacks, populations” following Qi et al. (2018). Clean reads were aligned against the obtained GBS reference genome using Bowtie 2 version 2.3.4.1 with default parameters (Langmead and Salzberg, 2012), and subsequently genotyped and screened using GATK version 3.8–1 (McKenna et al., 2010) for SNP and INDEL predictions. Finally, the obtained SNPs were filtered using VCFtools version 0.1.16 (Danecek et al., 2011) with option “–maf 0.05 –max-missing 0.8 –minDP 4 –min-alleles 2 –max-alleles 2.”

Genetic Diversity, Population Structure, and Phylogenetic Analysis

The number of private alleles, expected and observed heterozygosity (H and H), polymorphism information content (PIC), nucleotide diversity (π), the effective number of alleles (N), and Wright’s F statistics [fixation index (F) and inbreeding coefficient (F)] were calculated using the module “population” in Stacks version 2.4, VCFtools version 0.1.16, and genepop version 1.1.4 (Rousset, 2008). Reynold’s genetic distance (DR) among populations was determined from F as follows:

Evolutionary clusters were identified using Admixture version 1.3.0 (Alexander et al., 2009) with the default parameters. The predefined genetic clusters (K) ranged between 2 and 10, each of which was repeated 10 times. The optimal value of K was determined using cross-validation (CV) error, which has the minimum error in Admixture. Principal component analysis (PCA) was performed by GCTA version 1.26.0 (Yang et al., 2011) to explore the genetic structure of the species. A maximum likelihood-based phylogenetic tree was performed using IQ-TREE 2 (Minh et al., 2020) to clarify the genetic relationships among populations with the coalescent SNP dataset, and node support values (ML) were calculated with 1,000 ultrafast bootstrap replicates.

Demographic History Inference

The demographic history of L. oblata was inferred using Stairway Plot 2 (Liu and Fu, 2020). The mutation rate was set as 7.7e-9 per site per year following Pu et al. (2020). The generation time was set as seven, based on our observations from the cultivated individual since its seed stage in 2016. The tallest individual is about 0.5 m in height, and it may bloom in the spring of 2022. Though the RAD sequencing did not cover the full genome, the obtained SNP are most likely a part of the latter, which may reflect the trend of its historical demography. Two representative populations, namely, the southmost large population HDL and the small LJZ located at central were employed for demographic history inference. All samples were treated as one group, a one-dimensional site frequency spectrum (1D-SFS) was constructed using ANGSD software (Korneliussen et al., 2014), and bootstrap iterations of 200 were implemented. The result was visualized in R.

Results

Single-Nucleotide Polymorphisms Discovery

A total of 625, 652, 430 raw reads (91.93 Gb) were generated, resulting in 594, 664, 390 clean reads (86.82 Gb) and an average of 0.65 Gb per sample (Supplementary Table 2). The average sequencing depth of all samples was 47.33 × , with a coverage range of 87.82–96.86%. A total of 33,758 SNPs were obtained and employed for downstream analysis.

Population Genetic Diversity and Differentiation

Table 1 reports the descriptive statistics of genetic diversity (H, H, PIC, π, and N) for the seven populations of L. oblata. PIC values in all seven populations were very low (<0.19), with a value of 0.2364 at the species level. Relatively higher values of genetic diversity were observed among the four large populations (JK and JMS from YMs, WTS, and HDL from THMs) at the marginal regions compared to the three small populations located at the center (DLS, LJZ, and BJS). The northwest population WTS exhibited the highest values of H (0.2291), H (0.2054), PIC (0.1853), π (0.2328), and N (1.3826) (Table 1), while the lowest values of H (0.1477), PIC (0.1172), π (0.1625), and N (1.2563) were determined for BJS, one of the smallest populations. Furthermore, a large number of private alleles were detected in the four large populations (N = 1581 in JK, N = 256 in JMS, N = 724 in WTS, and N = 1,489 in HDL), while none were detected in the three small populations.

TABLE 1

The statistics of the number of samples (N), expected heterozygosity (H), observed heterozygosity (H), polymorphism information content (PIC), nucleotide diversity (π), efficient allelic number (N), and the number of private alleles among populations (Pops).

Mountains	Pops	N	H E	H O	PIC	π	N E	Private alleles
Yan Mountains	JK	36	0.2015	0.1787	0.1618	0.2046	1.3421	1,581
	JMS	20	0.1928	0.1828	0.1550	0.1982	1.3270	256
Taihang Mountains	DLS	6	0.1541	0.1664	0.1232	0.1694	1.2645	0
	LJZ	8	0.1888	0.1857	0.1524	0.2027	1.3165	0
	BJS	6	0.1477	0.1932	0.1172	0.1625	1.2563	0
	WTS	34	0.2291	0.2054	0.1853	0.2328	1.3826	724
	HDL	24	0.2044	0.1794	0.1642	0.2091	1.3469	1,489
ALL		134	0.2863	0.1863	0.2364	0.2874	1.4554	4,050

A high level of genetic differentiation (F = 0.3245) was observed at the species level, along with a low inbreeding coefficient (F = 0.0986). The highest genetic differentiation was detected between the most northeast population JK in the YMs and the most southwest population HDL in the THMs (F = 0.3963), while the value was minimized between the two closest (geographically) neighbored populations located at the center of its distributional range, LJZ and DLS (F = 0.1859, Table 2). The two YM populations, JK and JMS, exhibited large genetic variations (F = 0.3099). The southernmost population HDL differed greatly from the other THM populations (F of 0.2626–0.3805). Notably, WTS, one of the four largest populations on the northwest side, exhibited less variation with its neighbors, namely, BJS, DLS, LJZ, and HDL, with F values varying from 0.2320 to 0.2723 (Table 2). The higher the F between two populations, the longer the DR (Table 2), which may also suggest a lower gene flow.

TABLE 2

Genetic differentiation coefficient (F) and Reynolds genetic distance (DR) among populations.

Pops	JK	JMS	LJZ	DLS	BJS	WTS	HDL
JK	–	0.3709	0.3683	0.4385	0.4974	0.3895	0.5047
JMS	0.3099	–	0.2497	0.3245	0.4460	0.3455	0.4861
LJZ	0.3081	0.2210	–	0.2057	0.3692	0.2657	0.4214
DLS	0.3550	0.2771	0.1859	–	0.4763	0.3179	0.4788
BJS	0.3919	0.3598	0.3087	0.3789	–	0.2640	0.4740
WTS	0.3226	0.2921	0.2333	0.2723	0.2320	–	0.3046
HDL	0.3963	0.3850	0.3439	0.3805	0.3775	0.2626	–

The lower triangle presents interpopulation F

Population Structure and Phylogenetic Analysis

The genetic structure of L. oblata was investigated by Admixture (Figure 2A). The CV analysis suggested an optimal K value of 7 (Supplementary Figure 1). However, the seven suggested clusters did not fully correspond to their seven natural populations. Admixture suggested that LJZ and DLS formed one cluster, and it identified JK to have a fine-scale genetic structure and divided it into two clusters (Figure 2A). The first principal component (PC1) of the PCA, which explained 15.79% of all genetic variance, was able to differentiate the seven geographic groups (Figure 2B). The four large natural populations were clearly identified, and the three small populations were plotted close to each other (Figure 2B). The SNP phylogeny revealed a strong phylogeographic pattern (Figure 2C), consisting of the YM and THM clades. Seven distinct lineages were detected, corresponding exactly to the seven natural populations with full clade supports (ML = 100) (Figure 2C and Supplementary Figure 2).

FIGURE 2

Genetic structure and phylogenetic relationships of the seven Lonicera oblata populations. (A) Admixture proportions of genetic clusters for each individual of the seven populations. The scenario of K = 7 is the best value according to cross-validation analysis. (B) Principal component analysis plot for the 134 L. oblata individuals based on the first two principal components. (C) A maximum-likelihood tree based on 33,758 single-nucleotide polymorphisms (SNPs) of the nuclear genome, with seven fully supported lineages (ML = 100), which exactly correspond to species’ natural populations. The two pink circles on the northside represent populations from Yan Mountains, and the five blue circles represent populations from Taihang Mountains.

Demographic History of Lonicera oblata

The variations of effective population size dating from one million years ago were inferred using Stairway Plot 2 (Figure 3). The two populations of L. oblata, namely, HDL and LJZ, both underwent drastic decline after the Middle Pleistocene Transition (ca. 1.2–0.7 million years ago, Clark et al., 2006). Although expanded later, the population size contracted dramatically again after the Last Glacial Maximum (LGM, ca. 22 ka years ago). A recent population expansion following the LGM was suggested for HDL, while the population size of LJZ was decreasing.

FIGURE 3

Demographic history of two populations of Lonicera oblata, namely, Heduling (HDL) and Lijiazhuang (LJZ) inferred by Stairway Plot 2. The x-axis indicates the time before the present, and the y-axis represents the historical effective population size.

Discussion

Genetic Composition and Variation

The genome-wide data generated by HTS are relatively easier to obtain than traditional markers, and the high number of informative loci provide us not only the solid basis for accurately characterizing the population structure (Nazareno et al., 2017) but also the exciting opportunity to address a wide range of factors for conservation genetics (Médail and Baumel, 2018; Li et al., 2020). We investigated the population genetic diversity of L. oblata, a critically endangered montane shrub endemic to North China, based on nuclear genomic data from the SuperGBS approach. The results reveal a low genetic diversity and strong genetic structure. PIC is an important index for the evaluation of genetic diversity and can be classified into three scales, such as 0 < PIC < 0.25, 0.25 < PIC < 0.5, and PIC > 0.5. The higher the PIC value, the higher the polymorphism information in the population (Botstein et al., 1980). In this study, low PIC (< 0.25) values were detected at the population and species level, suggesting a low degree of genetic diversity within L. oblata. H values determined here were markedly lower than other Lonicera species (e.g., H = 0.2863 in this study vs. H = 0.78 in L. maackii (Rupr.) Maxim.) (Barriball et al., 2015), as well as other endangered plants species [e.g., H = 0.3482 in Tetraena mongolica Maxim. (Cheng et al., 2020), H = 0.364 in Firmiana danxiaensis H. H. Hsue and H. S. Kiu (Chen et al., 2014; Table 1)]. Though only eight individuals were sampled, the LJZ population possesses a relatively high level of genetic diversity among the three small populations (JMS, DLS, and LJZ, Table 1). There are lots of mountains with altitudes of ca. 1,200 m in LJZ, whose ridges are connected. The vertical slopes between the top and the hillside in most mountains may decrease human interference and provide the opportunity for frequent genetic exchange within the population. Note that the genetic diversity of the central populations (BJS, DLS, and LJZ) was averagely lower than that of the four marginal populations (JK, JMS, WTS, and HDL) in L. oblata (Table 1), so the scenario of declining genetic diversity during the northward expansion of the species is not supported in this study. In contrast to the low genetic diversity, we detected pronounced genetic structure within L. oblata (F = 0.3245), implying a long history of independent evolution and less gene flow among populations.

The population genetic composition is generally a function of both the species-specific biological characteristics and the corresponding ecological factors. The low genetic diversity and high genetic differentiation detected for L. oblata may be attributed to the following factors. (i) The extremely low numbers both at the population and individual level. The only seven uneven populations are highly fragmented and scattered at the top of mountains in a broad geographic range, and the limestone-preferred habitat of the species has been disturbed by long-term anthropogenic activities throughout history. The strong genetic differentiation among populations of L. oblata implies that the phenomenon of sky islands, which is significant in southwest China, may also occur in the THMs and the YMs in North China. (ii) The combined side effect of its reproductive system and fruit dispersal characteristics. A mixed mating system of biased outcross and partial self-fertilization is observed in L. oblata, and insects are necessary for pollination (Wu et al., 2022), while pollen limitation occurs during its flowering period. The harsh climate conditions during its flowering period (strong wind, sudden drop in temperature, snow, and/or dusty wind in the early spring) greatly interfere with the activities of effective pollinators and the pollination success of L. oblata. The highly fragmented populations and persistent strong winds at the top of hills and cliffs further inhibit pollination. The red fleshy fruits of L. oblata are attractive to birds, thus providing the possibility of long-distance dispersal of seeds and genetic exchange among populations. However, only a small-bodied bird Zosterops was recorded to feed on the fruits of L. oblata during our field investigation in Beijing in July. Zosterops live in a distance of ca. < 2 km in thick forests during the breeding season (May to July). Extending the distributional range of L. oblata over the contribution of birds like Zosterops is difficult. (iii) Seedlings barely survive in the field, which further decreases the genetic exchange between populations. Seedlings were difficult to find in the community during our field investigations, and an obvious inverted triangle population pyramid was observed (unpublished data). Hence, multiple negative effects may contribute to the strong genetic structure of L. oblata, whose destiny may be further challenged by future climate changes (Wu et al., 2021).

Population Structure and Glacial Refugia

The genetic structure of a population provides further information on the evolution of a species. In this study, a significant phylogeographic structure was identified from the structure, principal component, and phylogeny analyses, suggesting the division of seven groups (Figure 2). However, the seven clusters identified in Admixture did not correspond to those in the latter two analyses. In particular, the Admixture analysis placed DLS and LJZ into a single cluster, while the eastmost large population JK was divided into two clusters (Figure 2A). When K = 5, samples of JK did not differentiate intrapopulations, and the three small populations, namely, BJS, DLS, and LJZ, compose one cluster (Supplementary Figure 3). Such a scenario suggests an Admixture identity of these three small populations distributed in the species’ central region. The presence of two clusters in JK under the scenario K = 7 is difficult to explain. The six samples from one of the JK clusters (JK20–JK25) were randomly selected in the field just like the other samples. There were no distinct characteristics (e.g., elevation, slope direction, and canopy density) identified for the six samples compared to the remaining samples. JK is strongly exposed to tourism activity compared to the other six populations. The Jiankou Great Wall, a grand section of the Great Wall, attracts thousands of tourists each year. Many L. oblata individuals grow right on the Great Wall and are threatened by tourists during their climb. Some deaths of individuals due to the direct damage from tourists were witnessed during our fieldwork. Future work will include fine-scale sampling and landscape genomic study in order to understand the fine-scale genetic structure within JK.

The PCA and phylogenetic analysis revealed the seven identified groups to fully correspond to the natural populations of the species. Clade support values do not generally reach high levels in the study of population genetic diversity, and the sample(s) of one population may sometimes be resolved in the lineage of another population, such as cushion willow (Chen et al., 2019), grape (Liang et al., 2019), Prunus (Liang et al., 2019), and Q. aquifolioides Rehder and E. H. Wilson (Du et al., 2017). In this study, we obtained a well-resolved phylogeny of L. oblata at the population level. The phylogenetic tree consisted of two key clades, namely, the THMs and YMs, and the seven natural populations were all resolved as monophyly with full clade supports on the molecular tree (ML = 100) (Figure 2C). Both the strong phylogeographic pattern and the fully supported clade values at the population level suggest the seven natural populations to be evolutionary significant units following the definition of Moritz (1994). Strong phylogeographic characteristics combined with great genetic variation among populations (Table 2) also imply these populations be defined as management units following Avise (2000). Among them, the four large populations (JK and JMS in the YMs and WTS and HDL in the THMs), which have plenty of private alleles, are of greater importance as primary conservation units. The high level of genetic differentiation and well-resolved phylogeny in L. oblata may also suggest that long-term vicariance exists among these highly fragmented populations, as well as a low gene flow.

The NNE trending THMs and EW trending YMs act as dispersal corridors not only between the northeast cool-temperate and southern warm-temperate forests but also among fragmented populations of endangered species in China. These two great mountains in North China may also provide potential refugia for species during the global climate changes of the Quaternary (Zeng et al., 2011; Hou et al., 2020; Lin et al., 2021). Our results demonstrate many of the extant populations of L. oblata are genetically isolated. A high level of genetic diversity and plentiful private alleles were detected in the four large populations located in the marginal region of the L. oblata distribution range (N = 1,581 in JK, N = 256 in JMS, N = 724 in WTS, and N = 1,489 in HDL, Table 1). In contrast, low genetic diversity and no private alleles were detected in the populations located in the central region (LJZ, DLS, and BJS). Our ecological niche modeling analysis (ENM) of L. oblata suggests a wide historical distribution range, covering its current range since the LGM, followed by a dramatic decline to the present (Wu et al., 2021). Affected by the intensive uplift of both the Qinghai-Tibetan Plateau during the Late Cenozoic (Li et al., 2001) and THMs and YMs during the Late Pliocene to Pleistocene (Wu et al., 1999), which intensified monsoon while enhanced aridity in the Asian interior and repeated glacial events during the Quaternary, historical climate and geographic characteristics of northern China varied greatly. The demographic history of L. oblata may also have suffered from these factors, particularly the LGM, and its population size decreased sharply (Figure 3). This might also explain that no private allele was detected in the three small current populations, namely, BJS, DLS, and LJZ. Taken together, the genomic data and ENM suggest the existence of LGM refugia for L. oblata in the THMs and YMs, and the hypothesis of in situ survival is supported. The two greatest mountains in North China, namely the THMs and the YMs, may have acted as “Noah’s Ark” for numerous plant lineages during the LGM, providing vital refugia for the postglacial preservation of biodiversity in North China.

Threats and Conservation Suggestions for Lonicera oblata

The endangerment of species may be a result of both internal and/or external factors. Our field investigation suggests the possible worsening of the L. oblata endangerment situation. Although it occupies a relatively large distribution range in North China, only seven populations and a total of ca. 1,000 individuals of the species were recorded. Its limestone-specific habitat may greatly restrict its expansion, and the limited populations are highly fragmented and affected by human activities (e.g., agriculture, tourism, and logging). Furthermore, pollen limitation and harsh climate conditions during its flowering period restrain its survival and reproduction (Wu et al., 2022). Our results reveal a low genetic diversity and high genetic differentiation, which may suggest a reduction in the species fitness. Moreover, potentially suitable regions for the expansion of the species are limited (Wu et al., 2021). The dry Chinese Loess Plateau to the west of the THMs, the cold and dry desert in the Yin Mountains to the north of the YMs, and the cold temperature in Northeast China all constrain the future expansion of L. oblata severely. Due to the great conservation achievements of forestry and ecology across China in recent decades, the majority of previously bare mountaintops are now densely occupied by woody plants. This may also act as a challenge to the expansion of L. oblata, which favors an open habitat. Thus, this critically endangered species endemic to North China is currently facing the plight of “nowhere to go” (Nogués-Bravo et al., 2007; Loarie et al., 2009).

The exploration of population genetic patterns can provide vital information for the conservation and management of threatened species. Considering the aforementioned threatening factors, several conservation suggestions are proposed here: (i) forceful in situ conservation should be provided at both the population and the individual level; (ii) JK, JMS, HDL, and WTS populations should be given priority among the seven conservation units, including the construction of mini-reserves for the former three that are not distributed in national parks or nature reserves, and forceful conservation actions should be performed; (iii) the collection and long-term preservation of seeds from multi-populations should be accomplished, with the timely determination of seed germination and artificial propagation techniques; and (iv) as the intersection regions between the north of the THMs and west of the YMs are predicted to provide stable climate conditions in the future (Wu et al., 2021), this area should be given the priority to experiments of reintroduction and ex situ conservation.

Conclusion

In this study, we performed a genome-wide population genetic investigation on L. oblata, a highly threatened montane shrub endemic to North China, whose natural distribution range lies in the contact zone between southern warm-temperate and northern cold-temperate forests. Based on the SuperGBS method, we determined a low level of genetic diversity, a high degree of genetic differentiation, and a strong phylogeographic structure. Sharp population size contraction was suggested, which occurred after the middle Pleistocene Transition and the LGM, respectively. The results from the current genomic data combined with our previous ENM analysis suggest the existence of LGM refugia in the THMs and YMs. Affected by multiple factors relating to geography, climate, and biological and ecological characteristics, the future of L. oblata is challenged by the “nowhere to go” scenario. Four key conservation units were identified for L. oblata, and several conservation suggestions were provided. Further work characterizing the whole genome of L. oblata and a fine-scale landscape genomic study should be performed especially on the JK population with the aim to understand how geographic and ecological factors shape its genetic structure and evolutionary history. Our population genomic study provides valuable information for the conservation and management of L. oblata and contributes to the further understanding of the role of the THMs and YMs in preserving the biodiversity in North China.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/genbank/, PRJNA787048.

Author Contributions

X-YM conceptualized, supervised the study, analyzed the data, wrote, and revised the draft. Y-MW, X-LS, LT, F-WL, X-FX, and YN contributed to the fieldwork and collected the materials. Y-MW contributed to the species distribution maps. All authors contributed to the article and approved the submitted version.

Conflict of Interest

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

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.