Genetic diversity and population dynamic of Ziziphus jujuba var. spinosa (Bunge) Hu ex H. F. Chow in Central China.

Phylogeographic research concerning Central China has been rarely conducted. Population genetic and phylogeography of Ziziphus jujuba var. spinosa (also called sour jujube) were investigated to improve our understanding of plant phylogeographic patterns in Central China. Single-copy nuclear gene markers and complete chloroplast genome data were applied to 328 individuals collected from 21 natural populations of sour jujube in China. Nucleotide variation of sour jujube was relatively high (π = 0.00720, θ w = 0.00925), which resulted from the mating system and complex population dynamics. Analysis of molecular variation analysis revealed that most of the total variation was attributed to variation within populations, and a high level of genetic differentiation among populations was detected (F st = 0.197). Relatively low long-distance dispersal capability and vitality of pollen contributed to high genetic differentiation among populations. Differences in the environmental conditions and long distance among populations further restricted gene flow. Structure clustering analysis uncovered intraspecific divergence between central and marginal populations. Migrate analysis found a high level of gene flow between these two intraspecific groups. Bayesian skyline plot detected population expansion of these two intraspecific groups. Network and phylogeny analysis of chloroplast haplotypes also found intraspecific divergence, and the divergence time was estimated to occur at about 55.86 Ma. Haplotype native to the Loess Plateau was more ancient, and multiple glacial refugia of sour jujube were found to locate at the Loess Plateau, areas adjacent to the Qinling Mountains and Tianmu Mountains. Species distribution model analysis found a typical contraction-expansion model corresponding to the Quaternary climatic oscillations. In the future, the distribution of sour jujube may shift to high-latitude areas. This study provides new insights for phylogeographic research of temperate plant species distributed in Central China and sets a solid foundation for the application of the scientific management strategy of Z. jujuba var. spinosa.

INTRODUCTION

The evolutionary history of plant species has emerged as a complex interaction of biogeography, climate change, and human forces (Feng et al., 2018; Paola et al., 2017). Phylogeography seeks to understand the comprehensive evolutionary history and distribution of organisms and has become a focus of evolutionary biology (Abbott & Comes, 2004; Riddle, 2010). Of particular importance to Chinese phylogeographic research, the Qinghai‐Tibetan Plateau (QTP; Gao et al., 2019; Khan et al., 2018; Qiu et al., 2011), Southwestern China (Du et al., 2017; Hou et al., 2018; Zheng et al., 2017), and Northern China (Bai et al., 2010; Zeng et al., 2016) have attracted more attention resulting from the substantial diversity of plant species, the topographically heterogeneous terrain and the complex climate conditions. These studies shed light on the interactions of geography, ecology, and climate in shaping the distribution and genetic pattern of plant species. However, phylogeographic research concerning Central China (such as the Loess Plateau and adjacent areas) has been rarely conducted (Wen et al., 2020). Some endemic species are narrowly and concentrated distributed in Central China, such as Xanthoceras sorbifolia Bunge (Zhu, 2016) and Elaeagnus mollis Diels (Du et al., 2020). Furthermore, the Quaternary glaciation cycles had little effect on this area, glaciers only formed in scattered areas of the Loess Plateau and lasted for a relatively short period (Zheng et al., 1998). Phylogeographic study of E. mollis, an endangered deciduous shrub species restricted distributed in the Loess Plateau, revealed multiple potential glacial refugia and allopatric divergence in Central China (Du et al., 2020).

Ziziphus jujuba var. spinosa, also known as sour jujube, is a deciduous shrub plant species belonging to Ziziphus, Rhamnaceae. Z. jujuba var. spinosa is widely distributed in the temperate region of China, especially in areas from the Loess Plateau to the Taihang Mountains, which corresponds to its climate and habitat optimum (Zhang et al., 2015; Zhao et al., 2021). It harbors significant ecological values and is typically used for soil and water conservation in Northern China as the high tolerance to drought and salt stress (Wang et al., 2018). The nutritional value of the fruit and the medicinal importance of the seed of sour jujube has kept its economic vitality in China for more than 2000 years (Qu, 1982). It is concluded that Chinese jujube (Z. jujuba Mill) is originated from sour jujube, and the evolutionary path may involve several different patterns (Liu, 1993; Peng, 1991). Genome‐resequencing SNP data have been utilized to investigate the genetic structure of sour jujube with Chinese jujube, and the clear separation of sour jujube and Chinese jujube, and the genome evolution and domestication history of sour jujube were highlighted (Guo et al., 2021; Huang et al., 2016; Shen et al., 2021). As the primitive ancestor of Z. jujuba, sour jujube is considered as a valuable gene pool for the genetic improvement of Chinese jujube (Zhang et al., 2015). Accompanied by the growth of the plantation of Chinese jujube and the overharvesting of the seeds, the distribution range of sour jujube is undergoing severe fragmentation and decrease (Gao et al., 2008). Previous research on sour jujube mainly focused on taxonomic classification, nutritional and medicinal ingredients, and responses to various abiotic stress (Kang et al., 2008; Li et al., 2017; Liu, 1993; Yan et al., 2017). The genetic diversity and genetic structure of sour jujube have not been well characterized, which hinders the comprehensive understanding and utilization of this important genetic resource. Zhang et al. (2015) investigated the genetic diversity and population structure of sour jujube using SSRs and found high level of genetic diversity (H E = 0.659 and H S = 0.674) and moderate differentiation among populations (F st = 0.091, R st = .068). However, the phylogeographic pattern and population demography of sour jujube corresponding to various climate conditions were still yet to be clarified. Phylogeography study of sour jujube may add acknowledgement of the species' evolutionary history in Central China and provide some clues for the potential migration between this area and other adjacent areas.

In the present study, single‐copy nuclear gene markers were utilized to investigate the genetic diversity and population structure of sour jujube distributed across China. Furthermore, chloroplast genome data was used to detect the phylogeographic history of sour jujube. The population dynamic of sour jujube responding to the Quaternary and future climatic changes was also documented. The aims of this study were as follows: (i) determine the genetic diversity and population structure of Z. jujuba var. spinosa, (ii) describe the phylogeographic history and the potential location of glacial refugia, and (iii) investigate the population dynamic of sour jujube in China. The present study will improve our knowledge of plant phylogeography in Central China and set a solid foundation for future management and utilization of this important germplasm resource.

MATERIALS AND METHODS

Population sampling and DNA extraction

In this study, 328 individuals were sampled from 21 sour jujube populations distributed in China. The geographic information of the populations sampled in this study was illustrated in Figure 1 and Table S1. We also sampled one Chinese jujube (Z. jujuba “Hupingzao”) and one Z. mauritiana individual as outgroup in the following analysis. Individuals sampled in the same population were at least 100 m apart. Total DNA was extracted for all sampled individuals from silica gel‐dried leaves using the modified cetyltrimethyl ammonium bromide (CTAB) method.

FIGURE 1

Sampled populations of sour jujube in this study

Locus selection, amplification and sequencing

Single‐copy nuclear gene markers have been widely used in plant population genetic studies, resulting from their cross‐utility, high variable, and easily processed features (Du et al., 2015; Hou et al., 2018; Liu et al., 2016; Wang et al., 2019). In this study, we selected five single‐copy nuclear gene markers developed in our previous research for amplification and sequencing of all the samples (Hu, Du, Wang, & Han, 2021). Detailed information about the locus was listed in Table 1. Polymerase chain reaction (PCR) amplification followed the protocols established in Hu, Du, and Han (2021). Direct sequencing was performed on an ABI 3730XL DNA analyzer (Applied Biosystems) using the same primers that were used for amplification.

TABLE 1

Primers for the single‐copy nuclear gene markers used in this study

Locus	Primer sequences (5′–3′)	Gene annotation	Ta (°C)
SZ2	F:TGGTACAGGATCTACAATTC R:CCTGACTTTCTAATTGCTTC	myb‐like protein X	52
SZ11	F:ATGGCTTTTGCTTGCCTCTC R:GTGCATTCGGGTCATCAATG	UDP‐N‐acetylmuramoyl‐L‐alanyl‐D‐glutamate‐‐2,6‐diaminopimelate ligase MurE homolog	53
SZX2	F:GCTACTCGCTCTGGTTTCCAT R:GAAGAATCCTTGCCGGTTCAG	Trichohyalin	54
SZX12	F:GCCCTTTCGCAAAGCTTTCTT R:CCCAACACTGAGATTACTGGAG	Transcription termination factor MTEF18	54
SZX13	F:AAATGGAAGCGGCCTAGTGA R:TCCAGGAGTTTCCTCAGAGTC	DNA repair endonuclease UVH1	54

Genetic diversity and the neutral test

The assembled contigs of each individual were aligned using Clustal X (Thompson, 1997) and refined manually in Bioedit (Hall, 1999). The number of segregating sites (S), nucleotide diversity parameters, π (Nei, 1987) and Watterson's θ w (Watterson, 1975), and the minimum number of recombination events (R m) were analyzed for the five loci using Dnasp 5.10.0 (Librado & Rozas, 2009). We also calculated Tajima's D (Tajima, 1989), Fu and Li's D* and F* (Li & Fu, 1993) for each locus using Dnasp 5.10.0 to test how well the data conformed to the neutral model of evolution. Furthermore, we performed a standard multilocus Hudson‐Kreitman‐Aquadé (MLHKA) test (Hudson et al., 1987) to test for selection. The sequences of Z. mauritiana were used as an outgroup for the MLHKA test.

Genetic differentiation

The distribution of nucleotide variation among populations at each locus was assessed using AMOVA (analysis of molecular variation) implemented in Arlequin 3.5 (Excoffier et al., 2005) with significance testing evaluated using 1000 permutations of the original data. Genetic variation was hierarchically partitioned into among populations within species and within populations. Pairwise Wright's fixation index, F st (Wright & Maxson, 1968), was used to measure population differentiation Furthermore, to assign the sampled individuals to different genetic clusters (K), genetic structure analysis was performed in Structure 2.3.4 (Falush et al., 2003) with these five loci. Ten independent runs for each possible value of K from 1 to 10 were performed with a burn‐in of 500,000 following 1000,000 MCMC iterations with an admixture model and correlated allele frequencies. The most likely value of K based on the negative natural log‐likelihood of the data (LnP(K)) and ΔK was calculated using Structure Harvester (Campana et al., 2011; Evanno et al., 2005). Distruct 1.1 (Rosenberg, 2004) was used to create and visualize the population bar plots.

Chloroplast genome analysis

To understand the phylogeography history of sour jujube in China, we sequenced and assembled the complete chloroplast genome of 21 individuals (randomly chose one individual from each sampled population). The chloroplast haplotypes were assigned using Dnasp 5.0, and the median‐joining (MJ) haplotype network was generated in Popart 1.7 (Leigh & Bryant, 2015) with the chloroplast genome sequence of Z. mauritiana (NC_037151) as outgroup. The divergence time of the chloroplast haplotypes was detected with BEAST 1.10 (Suchard et al., 2018). Complete chloroplast genome sequences of Z. jujuba (NC_030299), Z. mauritiana, Z. spina‐christi (NC_037152), and Morus indica (NC_008359) were utilized as outgroups. A likelihood ratio test (LRT; Felsenstein, 1988) of the data implemented in MEGA6 (Tamura et al., 2013) rejected the molecular clock hypothesis with p < .01. A GTR nucleotide substitution model determined by jModeltest 3.7 (Posada, 2008) and a Yule speciation process were used. Using complete genome sequences, Liu et al. (2014) revealed that Z. jujuba diverged from other Rosales species at about 79.9 Ma (million years ago, 95% HPD: 55.7–99.4 Ma) and this time with a standard deviation of 5 Ma was used to constrain the maximum root age of the phylogeny. Four independent MCMC replications were performed, and each replication was run for 100,000,000 steps and sampled every 1000 steps. The first 25% of the sampled trees were discarded as burn‐in. The program Tracer 1.7.1 (http://tree.bio.ed.ac.uk/software/tracer) was used to confirm that the runs had reached a stationary distribution with sufficient mixing. Figtree 1.3 (http://tree.bio.ed.ac.uk/software/figtree/) was used to visualize the phylogenetic tree.

Coalescent‐based population demographic history analysis

To reconstruct the demographic history of sour jujube and the two intraspecific groups revealed in Structure analysis (see below) over time, the historical demography was inferred from Bayesian skyline plot analyses (BSP) implemented in BEAST 1.10 (Suchard et al., 2018). This coalescence‐based approach utilized a standard MCMC sampling procedure to evaluate the posterior probability distribution of effective population size backward until the time to the most recent common ancestor of the sampled sequences was reached. All the five loci were used for this analysis. The best‐fit substitution model of the concatenated sequence was determined by jModeltest 2 (Posada, 2008). Independent MCMC analysis was run for 1 × 108 steps, sampling every 100 steps and discarding 25% of sampled trees as burn‐in. For each group, multiple independent analyses were performed with different random seeds to test for convergence, and results of replicate runs were pooled using LogCombiner1.10.4, and skyline plots were visualized with Tracer 1.7.1 (http://tree.bio.ed.ac.uk/software/tracer).

Gene flow analysis

An MCMC maximum likelihood method was utilized to estimate gene flow between two intraspecific groups using Migrate 4.4.3 (Beerli, 2001). Two important population genetic parameters, θ (four times effective population size multiplied by mutation rate per site per generation) and M (immigration rate divided by the mutation rate) were calculated with the implemented F st estimations. Five independent runs by 10 short chains of 5000 steps and 3 long chains of 50,000 steps were run. Genealogies with a sampling increment of 100 and 10,000 burn‐in were recorded.

Landscape genetics

The biogeographic boundaries were calculated by Monmonier's maximum difference algorithm in Barrier 2.2 (Manni et al., 2004) based on the geographic and genetic distance matrix. Permutation and bootstrap tests were conducted with 1000 replicates. Furthermore, to evaluate the effect of geography on genetic differentiation, the correlation between geographic distance (GeoD) and genetic distance (GenD) among populations was calculated by the Mantel test (Mantel, 1967) using Vegan package implemented in R (Simpson et al., 2010).

Impact of environmental factors on genetic structure (isolation by environment)

In order to evaluate the effect of present environmental conditions on the observed pattern of the genetic structure of sour jujube, the correlation between GenD and environmental distance (EnvD) among these sampled populations was calculated by Mantel test using Vegan package implemented in R. Twenty‐three environmental variables of the 21 sampled populations were determined for the current climate layers as used in species distribution modeling analysis (see below). EnvD matrix was computed based on the population score on each bioclimatic variable.

Species distribution modeling

Twenty‐three environmental variables with 30 s spatial resolution were used in species distribution modeling (19 bioclimatic variables and elevation data were directly downloaded from WorldClim dataset [www.worldclim.org], soil variable was downloaded from the National Geomatics Center of China, slope and aspect were generated from elevation variable through surface analysis in Arcgis 10.5 [ESRI, Environmental Systems Research Institute]). These environmental variables were held in a GIS as ESRI grids using China Map Grid coordinates. The multi‐collinearity test was conducted by using the Pearson correlation coefficient (r) to examine the cross‐correlation among the variables and variables with cross‐correlation coefficient value |r| ≥ .8 were randomly discarded one and kept the other to minimize the negative influence on modeling. Finally, 10 environmental variables including Bio3 (iso‐thermality [Bio2/Bio7] × 100), Bio7 (temperature annual range[Bio5‐Bio6]), Bio11 (mean temperature of the coldest quarter), Bio15 (SD of humidity seasonality), Bio18 (precipitation of the warmest quarter), Bio19 (precipitation of the coldest quarter), elevation, slope, aspect, and soil were kept and utilized in species distribution modeling (Table S2).

Species distribution modeling was performed with Maxent 3.4.4 (Phillips et al., 2006; Phillips & Dudík, 2008), which used the maximum entropy method to model species' distribution. Past (including the Last Interglacial, the Last Glacial Maximum, and middle Holocene), present day, and future (including 2050s and 2070s under RCP2.6 [Representative Concentration Pathways], RCP4.5, and RCP8.5) niche models for sour jujube were generated. Distribution data of sour jujube was collected based on the present study and retrieved in National Plant Specimen Research Center (http://www.cvh.ac.cn/) and distribution information listed in Zhang et al. (2015) and Zhao et al. (2021). A total of 519 distribution data were collected, and R package “spThin” was utilized to remove the samples that clustered within 10 km to reduce the sampling deviation impact (Aiello‐Lammens et al., 2015), and delete redundant data. To reduce sampling deviation, only one distribution point was kept for each grid (10 × 10 km), and finally 253 effective distribution records were used in the following analysis (Table S3). In Maxent analysis, 80% of the distribution localities were used to train the model and 20% were randomly selected to test the model. 10 replicate runs were performed, and replicated run type was set to cross‐validate. The random seed was used to ensure consistency in the statistical output between runs. Model performance was evaluated using the area under the receiver operating characteristic curve (AUC). The value of AUC varies between 0.5 and 1, and a higher value means higher model accuracy, although in practice the maximum possible value of AUC is often <1 (Phillips et al., 2006). The Jackknife procedure was used to assess the importance of the variables contributing to the model. Model predictions were visualized in Arcgis 10.5 (ESRI). Based on the definition of species' suitable distribution area in WCRPCMIP6 (World Climate Research Program Coupled Model Intercomparison Project Phase 6, https://esgf‐node.llnl.gov/projects/cmip6/), the value of distribution possibility was reclassified into three hierarchy in Arcgis 10.5: 0–0.33 means unsuitable distribution area; 0.33–0.66 means suitable area and 0.66–1.0 means higher suitable area.

RESULTS

Nucleotide variation and neutrality test

We successfully sequenced 5 single‐copy nuclear gene markers from all the individuals sampled in this study. As shown in Table 2, the length of the single‐copy nuclear gene markers ranged between 190 and 846 bp. The average number of segregating sites and haplotype diversity were 37 and 0.857, respectively. A total number of 237 SNPs (single nucleotide polymorphism) was detected through the five single‐copy nuclear gene markers. Two important nucleotide variation parameters, π and θ w, reached 0.000720 and 0.00925, indicating the high level of nucleotide variation of sour jujube. The average number of recombination events throughout the single‐copy nuclear gene markers was 8, which indicated that the sour jujube genome underwent frequent inter‐chromosome fusions and segmental duplication events (Liu et al., 2014). Results of Tajima's D, Fu, and Li's D* and F* showed significant departure from neutrality at some loci. The result of the MLHKA test showed no significant departure from neutral expectation (χ 2 = .010, p = .9200), indicating that the significant results of Tajima's D, Fu, and Li's D* and F* detected in these loci may result from the presence of recombination within loci or the population structure (Ramos‐Onsins & Rozas, 2002; Zheng & Ge, 2010).

TABLE 2

Result of nucleotide variation of five single‐copy nuclear gene markers and neutrality tests

Locus	L	S	H d	π	θ w	D	D*	F*	R m	N m
SZ2	309	24	0.611	0.00619	0.01164	−1.16927	−4.95722**	−4.07025**	3	1.16
SZ11	846	82	0.940	0.00750	d0.01694	−1.59189	2.18531**	0.28719	17	0.67
SZX2	190	28	0.704	0.00746	0.02347	−1.75806*	−0.27874	−1.13229	3	1.93d
SZX12	648	43	0.923	0.00503	0.01124	−1.50814	−4.08123**	−3.46958**	9	1.09
SZX13	730	60	0.929	0.00765	0.01321	−1.17321	−3.45451**	−2.81122*	8	0.84
Mean	—	37	0.857	0.00720	0.00925	—			8	1.33

p < .05.

p < .01.

Population structure and genetic differentiation

AMOVA analysis was used to investigate the overall distribution of genetic variation (Table 3). Variation among populations ranged from 9.06% at locus SZX2 to 27.96% at locus SZ11 and was significant (p < .001). Most of the total variation was derived from variation within populations (72.04%–90.94%). The average F st was 0.197, indicating the relatively high level of genetic differentiation among populations.

TABLE 3

AMOVA and F st analysis of Ziziphus jujuba var. spinosa populations in China

Locus	SZ2	SZ11	SZX2	SZX12	SZX13	Mean
Among populations (%)	24.37	27.96	9.06	15.84	21.28	19.70
Within populations (%)	75.63	72.04	90.94	84.16	78.72	80.30
F st	0.243	0.279	0.090	0.158	0.212	0.197

Following the method of Evanno et al. (2005), two significant genetic clusters were detected in Structure analysis (Figure 2). According to the estimated ancestry, most sampled individuals showed an average inferred major membership proportion higher than 0.60, which indicated that they can be classified as belonging to one of the two distinct genetic clusters. Among the populations in genetic cluster 1, most were distributed in the area from the Lose Plateau to the Taihang Mountains (83.33%), which was the distribution center of sour jujube in China. Therefore, these populations in cluster 1 were referred to as central populations in the following analysis. Populations in genetic cluster 2 mainly consisted of populations that were distributed far away from the central area (55.56%), and they were referred to as marginal populations. Furthermore, some populations may be classified as hybrid populations of the two genetic clusters, with the proportion of membership lower than 0.6 in both clusters (Zhang et al., 2015). These populations were mainly distributed in the adjacent area of central and marginal populations (Table S4). Migrate 4.4.3 was employed to estimate the historical gene flow between the central and marginal populations. Results showed that the effective population size of the central population was about 1.5 times that of the marginal population. Gene flow between these two groups was asymmetric (Table 4). Mantel test of relationship between GeoD, EnvD and GenD showed significance correlation (GeoD/GenD, r = .2619, p = .02; EnvD/GenD, r = .3873, p = .001). Analyses of landscape genetics using Barrier 2.2 identified four significant barriers among sampled populations (.25 < p < .41; Figure 3a). Most of the barriers appeared between the central and marginal populations, which further proved the intraspecific divergence within sour jujube. BSP analysis revealed that the effective population size of sour jujube and the two intraspecific groups underwent substantial expansion over time (Figure 3b–d).

FIGURE 2

Result of structure analysis (a: result of ΔK, b: genetic composition of Ziziphus jujuba var. spinosa individuals)

TABLE 4

Result of historical gene flow between population clusters

	θ	M
		CP → MP	MP → CP
Central population (CP)	0.01497 (0.00940–0.02160)	1696.7 (1432–1978.7)
Marginal population (MP)	0.00943 (0.00113–0.02)		499.3 (64–301.3)

FIGURE 3

Result of barrier analysis and geographic location of the genetic barrier was indicated by thick red lines (a) Bayesian skyline plots for the three groups, showing effective population size as a function of time. The upper and lower limits of light blue trend represented the 95% confidence intervals of HPD analysis. (b) Total sour jujube populations; (c) central populations; (d) marginal populations

We successfully sequenced and assembled the complete chloroplast genome of 21 sour jujube individuals randomly selected from each sampled population. The length of the complete chloroplast genome varied between 159,399 and 161,279 bp, and GC content ranged from 36.51% to 37.30%. Detailed information about the sour jujube chloroplast genome was listed in Table 5. A total of 21 chloroplast haplotypes (H1–H21) were designated and the outgroup Z. mauritiana contained a private haplotype (H22). The MJ network of the chloroplast haplotypes showed that the chloroplast haplotype native to population YCRC was ancestral, which was also supported by BEAST analysis (Figure 4). Both network and BEAST analysis revealed two intraspecific lineages within sour jujube, which was consistent with the result of single‐copy nuclear gene marker Structure analysis. The divergence time of these two intraspecific lineages was dated back to about 55.86 Ma (95%HPD: 31.82–79.37 Ma). Furthermore, in BEAST analysis, chloroplast genome of Chinese jujube occupied a terminal position in the phylogenetic tree, which indicated that Chinese jujube may originate or artificially domesticate from sour jujube. The exact domestication process may need further research with more data.

TABLE 5

Detailed information about the chloroplast genome of 21 sour jujube individuals

Individual	Length (bp)	GC content (%)	LSC (bp)/GC(%)	SSC (bp)/GC(%)	IR (bp)/GC(%)
HBCD19	160,876	36.92	88,526/34.7	19,392/30.9	26,479/42.6
HBQHD5	160,730	36.69	88,690/34.5	19,364/30.8	26,338/42.7
HBXT5	160,691	36.86	88,411/34.7	19,358/30.9	26,461/42.7
CZCZ1	160,750	36.86	88,378/34.7	19,378/30.9	26,497/42.6
DTYG9	160,795	36.84	88,481/34.7	19,356/30.8	26,479/42.6
JZTG11	160,764	36.69	88,737/34.5	19,353/30.8	26,337/42.7
LLXY7	159,754	36.89	87,701/34.8	19,107/30.8	26,470/42.9
YCRC14	161,259	36.74	88,945/34.5	19,356/30.9	26,479/42.6
SDLS10	161,211	36.75	88,897/34.6	19,356/30.9	26,478/42.6
LYMY19	160,705	36.70	88,674/34.5	19,359/32.4	26,336/43.0
SDTA10	161,278	36.77	88,964/34.6	19,356/30.9	26,479/42.6
TJWQ2	161,100	36.98	88,862/34.6	19,366/30.9	26,436/42.6
NXYC4	161,062	36.80	88,746/34.6	19,358/30.9	26,479/42.6
SXHZ9	161,279	36.76	88,888/34.6	19,361/30.9	26,515/42.6
SXWN9	160,918	36.79	88,640/34.6	19,356/30.9	26,461/42.7
HNNY2	159,994	36.99	87,751/34.9	19,259/31.0	26,492/42.6
HNSMX5	160,347	36.93	87,995/34.8	19,368/30.9	26,492/42.6
HNLY9	159,399	37.30	87,244/35.5	19,127/31.1	26,514/42.6
HNZZ14	161,100	37.01	88,862/34.6	19,361/30.9	26,436/42.6
AHLA1	160,668	36.51	88,573/34.2	19,339/30.5	26,378/42.7
HBXY14	161,080	36.85	88,751/34.6	19,361/30.9	26,484/42.7

Abbreviations: LSC, long single‐copy region; SSC, short single‐copy region; IR, inverted repeat sequence; LGM, Last glaciation maximum.

FIGURE 4

Network and BEAST analysis of sour jujube chloroplast genome sequences (a: network analysis; b: BEAST result. Number next to nodes indicated divergence time, and dark bar indicated 95% HPD)

Distribution model

Using 10 environmental variables, the population distribution dynamic of sour jujube in different periods in China was modeled with Maxent (Figure 5). The result showed that the average value of the AUC of training data in different periods was 0.943 while that of the test data were 0.932, meaning higher accuracy of the distribution model (Table S5). The contribution rate of Bio11 was the highest (35.4%), and soil was the least (0.1%) to the model in the Jackknife procedure (Table S6). Because the average contribution rate of Bio11, Bio18, Bio7, and elevation to the Maxent model was higher than 10% and the accumulated contribution rate of these four variables was >80%, these four variables were classified as dominant variables affecting the distribution of sour jujube in different periods.

FIGURE 5

Distribution dynamic of Ziziphus jujuba var. spinosa in different periods (a: present; b: last glacial maximum; c: middle holocene; d: last interglacial; e: 2050s, RCP2.6; f: 2050s, RCP4.5; g: 2050s, RCP8.5; h: 2070s, RCP2.6; i: 2070s, RCP4.5; j: 2070s, RCP8.5)

The present potential distribution of sour jujube predicted in Maxent was highly consistent with the actual species' distribution (Wu, 1982). During the LGM, the range of sour jujube contracted considerably and showed a southward range shift. The proportion of the higher suitable area decreased from 2.08% at present to 0.80%, and the species was probably restricted to areas between the south of the Loess Plateau and the north of the Qinling Mountains, and areas adjacent to the Tianmu Mountains in East China. During the MH, the total suitable area of sour jujube increased from 9.79% in LGM to 11.45% and the higher suitable area mainly located in Shandong, Hebei, Shanxi, and Shaanxi provinces. In LIG, with the improvement of the climate, the area of total suitable area of sour jujube increased to 17.84% and the higher suitable area mainly located in Central China (from the Loess Plateau to Shandong Peninsula) and around Bohai gulf. In the future (2050s and 2070s), the area of the total suitable area of sour jujube decreased but that of the higher suitable area increased. The distribution range of sour jujube clearly showed a northward shift to high‐latitude areas in future environmental conditions (Figure 5 and Table S6).

DISCUSSIONS

Nucleotide variation of sour jujube in China

With the development of sequencing and analytical technology, genome data, such as genome‐resequencing and RAD‐seq (Restriction site‐associated DNA sequencing) data containing tens of thousands of SNPs have been widely applied in population genetic research (Lange et al., 2021; Mu et al., 2020). Using genome‐resequencing SNP data, the nucleotide variation of sour jujube was found to be higher than that of Chinese jujube (Guo et al., 2021; Huang et al., 2016). Though the number of SNPs in the present study (237 SNPs in five single‐copy nuclear gene markers) was small, but the basic population genetic information of sour jujube was revealed to some extent and added to our acknowledgement of this important plant species. Future genome data may apply in the genetic‐related study of sour jujube, such as the domestication history of sour jujube and the underlying mechanism for metabolite differences between sour jujube and Chinese jujube.

Although nucleotide variation of sour jujube at the locus level varied substantially, sour jujube distributed in China showed a relatively high level of nucleotide variation (π = 0.00720, θ w = 0.00925). Nucleotide variation of sour jujube was higher than that of other species distributed in China, such as Populus davidiana (π = 0.00440, θ  = 0.00750; Du et al., 2015), Tamarix austromongolica (π = 0.002 59; Wen et al., 2020) and Pugionium cornutum (π = 0.00532; Wang et al., 2013), and even comparable to some annual plant species such as Hordeum vulgare ssp. spontaneum (θ w = 0.0109), Zea mays (π = 0.0086–0.0133, θ w = 0.0094–0.015; Wright & Gaut, 2005), Arabidopsis thaliana (π = 0.0080–0.0165, θ w = 0.0049; Ramos‐Onsins et al., 2010). This was consistent with research on the genetic diversity of sour jujube using other genetic markers (i.e., SSR and RAMP; Zhang et al., 2015; Zhang et al., 2014). Furthermore, genome‐resequencing data also showed that the nucleotide diversity of sour jujube was higher than that of some perennial crops, such as peach (Li et al., 2019), but lower than that of date palm (Hazzouri et al., 2015), apple (Duan et al., 2017), and pear (Wu et al., 2018). The variation between this study and other research was partly attributed to the different parts of the genome analyzed.

A series of factors may contribute to the nucleotide variation of plant species, including representative sampling, natural selection, mating system, and demographic history (Wright & Gaut, 2005). Representative sampling can be ruled out to a great extent because the sampling area of this study covered the main distribution range of sour jujube in China, especially the areas from the Loess Plateau to the Taihang Mountains. Natural selection can also be excluded because the test statistics indicated no significant deviation from neutrality. Sour jujube was obligated to outcross and inbreeding seldom occurred in natural populations (Zhang, 2013), which contributed to a relatively high recombination rate and effective population size, and ultimately to a high level of nucleotide variation (Charlesworth, 2003). Considering how living organisms responded to global climate fluctuations, specially the glacial–interglacial cycles in the Quaternary (Qiu et al., 2011), we primarily speculated that sour jujube may undergo a process of population contraction and expansion during the evolutionary history. But the BSP analysis revealed that the effective population size of sour jujube underwent a continuous increase. This paradox can be clarified that the Quaternary glaciation had little effect on areas in Central China, especially the Lose Plateau (Zheng et al., 1998). In other central distribution areas of sour jujube, such as the Taihang Mountains, the population of plant species often migrated to high‐altitude areas to prevent severe bottlenecks from occurring, and maintained a higher effective population size (Hewitt, 2000). Consequently, the mating system and complex population dynamics have acted to retain a high level of genetic diversity in sour jujube distributed in China.

Genetic structure of sour jujube in China

In the present study, AMOVA analysis revealed that most of the total variation was attributed to variation within populations and a high level of genetic differentiation was detected (F st = 0.197, Table 3). This indicated that gene flow among populations was restricted to a relatively low level. Zhang et al. (2013) investigated the genetic structure of 3 sour jujube populations along the Yellow River and found a high level of gene flow and a low level of genetic differentiation among these populations. They speculated that seed dispersal depending on the water flow of the Yellow River may accelerate gene flow among these populations. They further investigated the genetic structure of 34 sour jujube populations and also found moderate differentiation (F st = 0.091; Zhang et al., 2015). The populations sampled in Zhang et al. (2015) were mainly distributed along the reach of some rivers, such as the Yellow River, the Jing River and the Luo River, seed dispersal by the river over long distance was the most probable explanation for gene exchange. However, there was no obvious connection among these populations sampled in the present study. Furthermore, the Mantel test revealed a significant correlation between geographic, environmental distance, and genetic distance. Differences in the environmental conditions and long distance resulted in a relatively high level of genetic differentiation among the sampled populations in this study. The long‐distance dispersal capability and the vitality of sour jujube pollen were relatively low (Shao et al., 2020). In general, it appears that forces such as isolation and adaptation that tend to increase genetic differentiation have been much stronger than homogenizing forces such as gene flow (Savolainen et al., 2007). Our results revealed a central, namely the Loess Plateau‐Taihang Mountains populations, genetic group surrounded by marginal populations that exhibited deep phylogeographic divergence as revealed in Network and phylogenetic analysis. The Migrate analysis on group‐level revealed a relatively high level of gene flow between these two intraspecific groups (Table 4). The existence of putatively hybrid populations may increase the estimation of gene flow. Furthermore, megafauna capable of long‐distance dispersal of seed cannot be completely ruled out, which played a certain role in shaping the genetic structure of some plant species (Feng et al., 2018). It is not yet possible to untangle natural gene flow over evolutionary time and recent human‐mediated gene flow to determine what traits characterized sour jujube population structure. Our analysis increased knowledge of the genetic diversity and population structure of sour jujube in China, revealed the potential geographic location of different gene pools, and could help to rationalize and prioritize reservoirs of genetic diversity.

Population demography of sour jujube in China

Zhao et al. (2021) compared 48 parameter combinations in R package ENMeval and Biomod2 to predict the distribution demography of sour jujube and found the Maxent model was optimal. In the present study, Maxent analysis yielded relatively accurate results that could be utilized to interpret the population demography of sour jujube corresponding to various environmental conditions in different periods. It was found that Bio11 contributed the most to the Maxent modeling, which reflected that low temperature was the main factor limiting the distribution of sour jujube. Sour jujube was not sensitive to soil conditions and its tolerance to barren soil was high (Wang et al., 2018; Zhao, 2016), which resulted in the lowest contribution of soil to Maxent modeling. The effects of Quaternary glaciation‐interglaciation cycles on the distribution pattern of plant organisms in the Northern Hemisphere have been substantially documented (Hewitt, 2004; Qiu et al., 2011). It has been proposed that temperate forests retreated southward to approximate 30°N during the Quaternary glacial periods based on paleovegetation data from East Asia and the populations located in northern areas must have recolonized from southern glacial refugia (Qian & Ricklefs, 2001). The paradigmatic contraction–expansion model of the latitudinal shift of temperate plant species responding to Quaternary climate oscillations has been drawn based on numerous phylogeographic surveys (Gonzales et al., 2008; Lafontaine et al., 2010). During the LGM, it was found that the distribution of sour jujube in China illustrated a typically southward retreat to about 30°N, and the higher suitable area mainly located in areas between the southward of the Loess Plateau and the northward of the Qinling Mountains, and areas adjacent to the Tianmu Mountains in East China. As mentioned above, glaciers only developed in southern scattered areas of the Loess Plateau and lasted for a relatively short time (Zheng et al., 1998). The mountainous areas were affected little by the Quaternary glacier and the climatic condition was relatively stable in the high altitudinal area. Therefore, appropriate climate conditions in these regions can serve as the highly suitable areas for sour jujube to survive the glaciation. During MH and LIG, the distribution of sour jujube showed a northward recolonization and expansion resulting from the improvement of climatic conditions (Yu et al., 2000; Zhao & Piperno, 2000). In the future (2050s and 2070s under 3 RCP conditions), the distribution of sour jujube in China showed northward migration to high latitudinal areas. Compared with the present distribution, the area of the higher suitable area increased but that of the total suitable area decreased. This was consistent with the research that under warming climate conditions, the distribution area of plant species could reduce and shift to high latitudinal areas (Ji et al., 2020; Walther et al., 2005). By simulating and predicting the distribution range of sour jujube in different periods, the historical causes of its formation and its future population dynamic responding to climate changes were deduced. Awareness of the possible effect of past climatic changes on the current population may provide insight into this species' future range dynamics in the light of climatic changes and be useful for germplasm management strategies (Paola et al., 2017; Ren et al., 2017).

Phylogeographic history of sour jujube in China

The effects of Quaternary climatic oscillations on the phylogeographic structure of species in the mid‐ to high‐latitude regions of Europe and North America (Emerson et al., 2010; Ren et al., 2017), high‐altitude areas (French et al., 2019; Gao et al., 2019; Khan et al., 2018), and Northern China (Hou et al., 2018; Zeng et al., 2016) have already been substantially documented. However, few studies have focused on the phylogeographic history of species native to temperate China. Therefore, the study presented here provided an opportunity to uncover the detailed population demographic history of a widely distributed plant species in temperate China and to better understand the processes playing roles in shaping the distribution pattern.

The fossil record was the primary clue to detect the distribution history of plant species. However, the fossil record for Ziziphus was not as complete or reliable as that of other species in Rhamnaceae (Burge & Manchester, 2008). Fossil records of Paliurus, the closest family to Ziziphus within Rhamnaceae based on molecular phylogenetic and morphological analysis, were numerously explored in North America, Europe, and Asia (Han et al., 2016; Islam & Simmons, 2006). The earliest identical fruit fossil record of Paliurus in Eastern Asia dated back to about the middle to late Eocene (Han et al., 2016). Although it cannot precisely determine the divergence time of these two closely related genera, the appearance of the diagnostic characteristics of modern Paliurus implied that species of Ziziphus were relatively old. Based on genomic sequences of species from Rosales, Liu et al. (2014) found that Z. jujuba diverged from other Rosales species at about 79.9 Ma. It was revealed that the two intraspecific groups of sour jujube diverged at about 55.86 Ma. The initial desertification in the Asian interior is thought to be one of the most prominent climate changes in the Northern Hemisphere during the Cenozoic era (An et al., 2001; Guo et al., 2002; Ruddiman & Kutzbach, 1989). Regional tectonic changes and ongoing global climate oscillations were probable causes of the intraspecific divergence revealed in sour jujube. Detailed information about this may need further investigation with more data, such as genomic SNPs.

The uplift of mountain regions may provide suitable habitats for plant species to survive during the Quaternary glaciation. It has been widely accepted that the Zhongtiao Mountains locating in areas between the Loess Plateau and Qinling Mountains in Central China harbor relatively rich vegetation diversity and serve as main glacial refugia for temperate plant species in the Quaternary glaciation, such as Acer mono Thunb. ex Murray (Liu, 2013) and X. sorbifolium (Zhu, 2016), resulting from appropriate temperature and adequate precipitation (Wang et al., 2004). In Maxent analysis, it was revealed that the higher suitable distribution area of sour jujube in China located at the Loess Plateau in different periods. During field sampling, it was found that sour jujube was distributed continuously in this area. Furthermore, nucleotide variation of central populations distributed in the Loess Plateau was much higher than that of the marginal populations (Table S8). Network analysis found that the chloroplast haplotype located in the Loess Plateau was more ancient. Combining these findings, we can speculate that a potential glacial refuge for sour jujube may locate in areas between the Loess Plateau and the Qinling Mountains, probably in areas adjacent to the Zhongtiao Mountains. Another higher suitable distribution area of sour jujube in LGM is located at areas adjacent to the Tianmu Mountains, which was also thought to serve as the main glacial refuge for temperate plant species in China, such as Ostryopsis davidiana Decne. (Tian et al., 2009), Juglans mandshurica Maxim. (Bai et al., 2018) and Ginkgo biloba L. (Zhao et al., 2019). However, in the present we did not collect samples from this area, maybe because the environmental conditions were not suitable for sour jujube, such as high moisture. These mountainous areas were important for the glacial survival of sour jujube and for preserving most of the extant nucleotide variation. This provides additional evidence for the in situ survival of plant species in Central China during the Quaternary glaciation and local expansion during interglacial or postglacial, which was similar to the scenarios revealed in other regions (Wang, Abbott, et al., 2010; Wang, Ikeda, et al., 2010). The most important refugia for sour jujube appear to have been located at the Loess Plateau in Central China and areas adjacent to the Tianmu Mountains in East China, which provides new insights for phylogeographic research of temperate plant species distributed in Central China. Furthermore, samples from other sour jujube distribution areas, such as Inner Mongolia Autonomous Region and Liaoning Province were not collected in the present study, but based on the above analysis, we may speculate that the populations distributed in these areas may also migrate or colonize from the potential refugia or central populations revealed in the present study. Future research may include samples from these areas to explicitly test our speculation.

CONCLUSION

This study provided new insights into the genetic diversity, population structure, and demography of Z. jujuba var. spinosa. The genetic diversity of sour jujube was relatively high. Results of AMOVA showed that most of the total variation was attributed to variation within populations and high genetic differentiation among populations was detected. Geographic distance and environmental difference contributed to the genetic differentiation among populations. Structure analysis successfully uncovered two intraspecific genetic groups, central population and marginal population corresponding to their geographic locations. Population demographic history of sour jujube illustrated a contraction‐expansion model responding to the Quaternary climate oscillations. Under future environmental conditions, the distribution habitat of sour jujube may shift to high altitudinal areas. Two potential glacial refugia were uncovered. We are confident that the information provided in the present study will be very helpful to the sour jujube management, conservation, and breeding activities.

AUTHOR CONTRIBUTIONS

Shuhui Du: Formal analysis (equal); writing – original draft (lead). Xiaoyan Hu: Formal analysis (equal); writing – original draft (supporting). Wendong Yu: Investigation (equal). Zhaoshan Wang: Investigation (equal).

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.

Appendix S1.

Click here for additional data file.