Predicting potential distribution of Ziziphus spinosa (Bunge) H.H. Hu ex F.H. Chen in China under climate change scenarios

Abstract Ziziphus spinosa (Bunge) H.H. Hu ex F.H. Chen is a woody plant species of the family Rhamnaceae (order Rhamnales) that possesses high nutritional and medicinal value. Predicting the effects of climate change on the distribution of Z. spinosa is of great significance for the investigation, protection, and exploitation of this germplasm resource. For this study, optimized maximum entropy models were employed to predict the distribution patterns and changes of its present (1970–2000) and future (2050s, 2070s, and 2090s) potential suitable regions in China under multiple climate scenarios (SSP1‐2.6, SSP2‐4.5, SSP3‐7.0 & SSP5‐8.5). The results revealed that the total area of the present potential suitable region for Z. spinosa is 162.60 × 104 km2, which accounts for 16.94% of China's territory. Within this area, the regions having low, medium, and high suitability were 80.14 × 104 km2, 81.50 × 104 km2, and 0.96 × 104 km2, respectively, with the high suitability regions being distributed primarily in Shanxi, Hebei, and Beijing Provinces. Except for SSP‐1‐2.6‐2070s, SSP‐5‐8.5‐2070s, and SSP‐5‐8.5‐2090s, the suitable areas for Z. spinosa in the future increased to different degrees. Meanwhile, considering the distribution of Z. spinosa during different periods and under different climate scenarios, our study predicted that the low impact areas of Z. spinosa were mainly restricted to Shanxi, Shaanxi, Ningxia, Gansu, Liaoning, Inner Mongolia, and Jilin Provinces. The results of core distributional shifts showed that, except for SSP1‐2.6, the center of the potential suitable region of Z. spinosa exhibited a trend of gradually shifting to the northwest.


| INTRODUC TI ON
Climate change is considered to be a key factor in altering the geographical distribution of species in the 21st century (Record et al., 2013;Santos-Hernández et al., 2021;Waltari et al., 2007). According to the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC), the average global surface temperature is expected to rise by 0.3-4.8°C by the end of the 21st century due to the continuous increases in greenhouse gas emissions (Braunisch et al., 2013;Cahill et al., 2013;Carroll et al., 2010). In response to this warming trend, many studies have suggested that species would change their currently suitable habitats in response to changes in environmental conditions, particularly as species distribution increases in elevation and migrates to northern latitudes (Heikkinen et al., 2006;Kujala et al., 2013;Qiu et al., 2011;Schweiger et al., 2008). In addition to changes in spatial habitats, climate change is modifying sensitive ecological responses, including flowering periods and the duration of growing seasons (Hampe et al., 2013;Wang et al., 2013).
With the emergence of novel computational statistics technologies and the development of the Global Information System (GIS), direct correlations between environmental factors (e.g., climate, topography, meteorological data, species data) have become possible, which is extensively used in ecological applications (Ye et al., 2020).
Ecological niche models (ENMs), also known as species distribution models (SDMs) (Brown, 2014;Brown et al., 2017), are employed to estimate the relationships between species presence and environmental factors through the extrapolation of multiple algorithms at multiple temporal and spatial scales (Conolly et al., 2012), which can be used to predict the potential distribution of species (Conolly et al., 2012). Over the last few decades, ENMs have played an important role in predicting the potential geographic distribution of species in the context of climate change, and have been broadly used in the domains of biology, ecology, biogeography, evolutionary biology, and species conservation (Araújo & Guisan, 2006). Among various ENM/SDM methodologies, maximum entropy (MAXENT) modeling has exhibited a better predictive ability, to become one of the most widely used models at present (Phillips et al., 2017;Radosavljevic & Anderson, 2014;Zeng et al., 2016). To date, the MAXENT model has been used to predict trends in the potential habitats of many plant species (particularly endangered species), such as Mimusops laurifolia (Forssk.) Friis (Hall et al., 2010) and Semiliquidambar cathayensis H.T.
Chang (Ye et al., 2020).  . Ziziphus spinosa possesses high nutritional, economic, and medicinal value, as its pulp is rich in sugars, acids, proteins, and vitamins, has a long flowering period, and can be a source of nectar. Moreover, Ziziphi Spinosae Semen (i.e., dry mature seeds of Z. spinosa) has such functions as nourishing the heart and liver, and treating insomnia Song et al., 2020). The past decade has witnessed an increasingly high market demand and price for Z. spinosa due to limited yields and supplies; thus, it has become urgent to promote its artificial planting and development.
To date, previous investigations of Z. spinosa have focused primarily on cultivation technologies and the pharmacological effects of Ziziphi Spinosae Semen Song et al., 2020). The present study employed an optimized maximum entropy model to predict and analyze the distribution areas of Z. spinosa under both present  and future (2050s, 2070s and 2090s) climate scenarios (SSP1-2.6, SSP2-4.5, SSP3.70 & SSP5-8.5).

| Collection and screening of sample data
Over the last 3 years (2019-2021), our research group conducted extensive field surveys in Shaanxi, Shanxi, Gansu, Hebei, Henan, and Shandong Provinces, and collected a total of 106 occurrence points.
In addition, 321 occurrence points were obtained from the previ-  (Figure 1).

| Variable environment screening and data processing
Bioclimatic variables are the major determinants of the ENMs of species and are frequently used for the development of plant ENMs (Tang et al., 2018). The set of 19 bioclimatic variables used in the present study were downloaded from the WorldClim website (https:// www.world clim.org), which involved a recent  and three future periods (2050s, 2070s, and 2090s). Considering the impacts of climate change scenarios on the accuracy of model development (Santos-Hernández et al., 2021), we selected four shared socioeconomic pathways (SSPs; SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) for three general circulation models (GCMs; BCC-CSM2-MR, CNRM-CM6-1, and MIROC-ES2L) in future climate data (Xu et al., 2020). Consequently, a total of 37 sets of bioclimatic data were included in this study, with one recent, and 36 future sets. Secondly, pairwise Pearson correlation coefficients (R) among the 19 bioclimatic factors were calculated using ENMTools v1.4.4 , and any pair of factors with |R| ≥ 0.8 were considered to be significantly correlated . Finally, for each pair of significantly correlated variables, only the ones with higher contributions to the model were retained (Cahill et al., 2013;Václavík & Meentemeyer, 2009Wisz et al., 2013;Zeng et al., 2016).

| Model establishment, optimization, and evaluation
Maxent v3.4.1 (Phillips et al., 2017) was used to construct the maximum entropy model for this study. Considering that the selection of general circulation models (GCMs) would lead to the uncertainty of the prediction results, we carried out arithmetic average processing on the prediction results of the three GCMs for the future periods.
To ensure the probability of Z. spinosa distribution being close to the normal distribution, we selected 70% of the data for model training and the remaining data for model testing (Phillips & Dudík, 2008). The other major parameters were set as follows: <Maximum Iterations: 5000; Replicated Run Type: Crossvalidate; No. of Replicates: 10>.
In this study, the R package <kuenm> (Cobos et al., 2019)  Thus, a total of 600 parameter combinations were multiplied by the FC and RM. On the basis of optimal model determination, the model (OR_AICc) with a statistically significant omission rate that was lower than the threshold value (0.05), and a delta AICc value of less than 2 was selected Ye et al., 2018).

| Classification of suitable region and reliability test of model
The suitability of species habitats is typically represented by the value range 0-1, where the higher the value is, the more suitable a certain area is for the species to grow. The selection of thresholds has an important impact on the prediction of suitable regions of different grades, which affects the calculation of different suitable areas (Arenas-Castro et al., 2020;Hu et al., 2020). Tang et al. (2018) proposed that a maximum test sensitivity plus specificity (MTSPS) threshold was superior to other threshold options in the grade division of suitable regions.
Thus, MTSPS was employed as the threshold value for this study, and those areas with suitability values lower than MTSPS were considered unsuitable for the growth of the species. The suitability range between the MTSPS and 1 was subdivided into three equal parts, which corresponded to the low, moderate, and high suitability regions, respectively (Li et al., 2019;Wang et al., 2017). The area sizes of different suitable regions as well as their changes in different future periods were calculated by using DIVA-GIS v7.5 (https://diva-gis.org).
Following model construction, the area under the receiver operating characteristic curve (AUC) was used to evaluate the accuracy F I G U R E 1 The occurrence data (406 points) of Ziziphus spinosa in China of the predictive model (Lobo et al., 2008). The mean AUC value was in the range of (0, 1), where AUC > 0.9 indicated that the model results were excellent and accurate (Warren & Seifert, 2010). Further, we considered the difference between the training AUC and the test AUC, where the smaller the absolute value of the difference, the higher the reliability of the model (Warren & Seifert, 2011).

| Analysis of low impact areas
Low impact areas refer to those where species are relatively less affected by climate change, which can be projected by superposing the binary prediction maps of suitable regions in different periods and taking the completely overlapping parts (Pan et al., 2020). In DIVA-GIS v7.5, the distribution maps of the potential suitable regions of different periods were overlaid to reclass the spatial units with suitability values greater than the MTSPS threshold as the suitable regions (Zhao, Zhang, et al., 2021). Those spatial units with suitability values lower than the MTSPS threshold were reclassed as unsuitable regions, which established the unsuitable and suitable matrices of Z. spinosa. Subsequently, the completely overlapping parts in the overlay layers were selected. After processing, the overlaid layers of different periods were imported into DIVA-GIS v7.5, and the potential low impact areas of Z. spinosa were visualized. This study predicted the low impact areas of four different shared socioeconomic pathways (SSP1-2.6, SSP2-4.5, SSP3-7.0 & SSP5-8.5) in current and future periods (2050s, 2070s, and 2090s).

| Analysis of spatial pattern changes
Spatial pattern changes refer to those of potential suitable regions of species across different periods, which could be obtained by superposing binary prediction maps of suitable regions during different periods (Santos-Hernández et al., 2021;Wu et al., 2021;Zurell et al., 2020). For this study, we created a prediction chart of the spatial pattern changes of four different SSPs (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) during current and future periods (2050s, 2070s, and 2090s), resulting in a total of 12 pattern change predictions. This was used to analyze the change rule of potential suitable regions of Z. spinosa in the recent period and three different future periods under various SSPs. In DIVA-GIS v7.5, the distribution maps of potential suitable regions of different periods were superposed to establish both the unsuitability and suitability matrices of Z. spinosa.
Based on the matrix table, the changes in spatial patterns of the suitable distribution regions under current and future climate change scenarios were further analyzed.

| Core distributional shifts
SDMToolbox V2.4 toolkit (Brown, 2014;Brown et al., 2017) of ArcGIS v10.2 was employed to calculate the variation trend of different suitable regions for Z. spinosa, and the central points of different regions were compared. We considered the suitable region of Z. spinosa as a whole, simplified it to a vector particle, and used the change of the centroid position to reflect the size and direction of the suitable region of Z. spinosa. Finally, the SDMToolbox toolkit was used to track the centroid of Z. spinosa, to investigate the distribution of the centroid during different periods and under various climate conditions, and to evaluate the migration distance of suitable regions via longitude and latitude coordinates (Smith et al., 2019;Zurell et al., 2020).

| Future potential suitable regions
The distribution and changes in potential suitability regions for Z.
spinosa during the three future periods differed under various climate scenarios; however, there were some similarities in the trends of changes (Figures 3-5; Tables 2 and 3). Under different climate scenarios, the area of low suitability regions increased significantly compared with its current value, while that of the moderate suitability regions decreased significantly. Conversely, the changing trends of total and highly suitable regions were not completely consistent under different climate scenarios (Tables 2 and 3).
Under the SSP1-2.6 scenario, the total area of the potential suitable region for Z. spinosa changed slightly (98.86%-104.62% of the current corresponding value), and showed only a small contraction in the 2070s. The area of the high suitability region increased to varying degrees. Notably, the area of high suitability region in the 2090s was 2.89 × 10 4 km 2 , with an increase of 201.26% compared with the current area.
Under the SSP2-4.5 scenario, the total area of the suitable region for Z. spinosa showed a slight increase (102.03%-107.46% of the current corresponding value). The area of the highly suitable region initially increased, contracted, and then finally increased. The values during the three future periods were 1.55 × 10 4 km 2 (2050s), 0.94 × 10 4 km 2 (2070s), and 1.85 × 10 4 km 2 (2090s), accounting Under the SSP5-8.5 scenario, the total area of the suitable region of Z. spinosa initially increased and then contracted, and the area of the suitable region was only 139.53 × 10 4 km 2 in the 2090s. In the 2090s, the area of the moderate suitability region was 31.16 × 10 4 km 2 (only 38.23% of the current corresponding value).
Except for the 2050s, the area of the high suitability region showed a gradual contraction trend. In the 2090s, the high suitability region was only 0.12 × 10 4 km 2 in area (12.82% of the current corresponding value).

| Low impact area
The prediction of the low impact area differed under various climatic scenarios ( Figure 6, Figure S1;

| Shift in the distribution center of the suitable region
The potential suitable region for Z. spinosa showed a trend of gradu-  . Previous studies revealed that the quality control of species occurrence points, screening of environment variables, selection of GCMs and SSPs, threshold selection, and the optimization of RM and FC all had significant effects on the predictive results of the model Yang et al., 2021;Ye et al., 2020;Zeng et al., 2016). For this study, these parameters were systematically optimized so as to ensure predictive accuracy to the maximum extent.

| Effects of environmental variables on species distribution
The geographical distribution of plants is restricted mainly by climatic variables, where hydrothermal conditions play a leading role in their distribution patterns . Precipitation is likely to increase or decrease as the climate changes and will affect soil moisture, which can cause plants to fail to reproduce, grow, and survive . Our study revealed that the main bioclimatic variables that affected the potential distribution of Z. spinosa were the mean temperature of the coldest quarter and the precipitation of the warmest quarter (Table 1). The mean temperature of coldest quarter for Z. spinosa is −9.1-5.8°C; it is much more likely to affect survival or cause a mismatch between phenology and season. The suitable range of the precipitation of the warmest quarter for Z. spinosa is 93-675 mm, which implies that Z. spinosa may prefer relatively dry environments. This was evidenced by our field surveys: Z. spinosa prefers arid and mild soil, and generally grows on sunny or semi-sunny slopes, with little requirement for water or fertilizers.
Furthermore, it is more suitable for temperate monsoon and temperate continental climates; thus, it appears that Z. spinosa is more suitable for planting in Northern China.
Under the SSP2-4.5 and SSP3-7.0 scenarios, the total area of the potential suitable regions for Z. spinosa increased by varying degrees. This indicated that within a certain range, the rising temperatures and increased rainfall caused by the greenhouse effect were more conducive for the growth of Z. spinosa. However, when the effects of global warming exceed the optimal tolerance range of species, the potential distribution region would also contract to varying degrees. For example, under the SSP5-8.5 scenario, the potential suitable region of Z. spinosa decreased significantly in the 2070s and 2090s.

| Changes in spatial patterns of potential suitable regions
In general, over time, the potential suitable regions for Z. spinosa gradually migrated to high-latitude areas (Figures 4 and 7). For example, under the SSP1-2.6 scenario, the suitable region of Z. spinosa increased, except for the 2070s; however, the distribution area changed. Jiangsu, Guizhou, Hubei, Anhui, southern Henan, and adjacent regions were close to being no longer suitable for its growth. Compared with previous studies , we also predicted the distribution of potential suitable areas in 2090s and the shift in the distribution center of the suitable region in different ages.
We found that during the period of SSP5-8.5-2090s, the potential suitable area of Z. spinosa would shrink by 14.19%; this indicated that with the intensification of the greenhouse effect, when the average temperature line exceeded the tolerance limit of Z. spinosa, the distribution area of Z. spinosa would gradually decrease and move to higher latitude. Therefore, by using the prediction results of the model, it is possible to take reasonable protection measures to Z. spinosa as soon as possible, which would alleviate the impact of climate change. In addition, the total area of potential suitable areas of Z. spinosa tends to increase, and the expansion area is mainly located in middle and high latitudes, while the decrease area is mainly located in low latitudes, which is consistent with the view of Zhao, Zhang, et al. (2021).

| Development and protection of germplasm resources
Ziziphus spinosa is a woody plant with high economic value that can tolerate cold and dry environments, and is suitable for planting in northern China Song et al., 2020). Furthermore, particularly in the Loess Plateau area, it can also be employed as a shelter forest species. Our study predicted that Shanxi, central and northern Shaanxi,

| Study limitations
For this study, we considered only the effects of bioclimatic variables on species distribution, which was also practically influenced by a variety of biological factors (e.g., interspecific competition, predation, and disease) and abiotic factors (e.g., soil, topography, and anthropogenic activities) . Moreover, our study assumed that species would have a sufficient dispersal capacity to migrate to any climatically suitable area under climate change.
However, we did not consider factors such as the species migration rate, as well as geographical and ecological isolation, all of which can lead to potential discrepancies between predicted and actual distributions. Elucidating the influences of all these factors requires a more comprehensive niche modeling approach, which has yet to be done in future studies (Wilting et al., 2010).
ENMs/SDMs make important assumptions about the relationship between species distributions and their environment that may limit their ability to accurately predict future species distributions. In particular, SDMs in theory assume stable fundamental niches, but in practice, they assume stable realized niches. The assumption of a fixed

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad at https://doi.org/10.5061/dryad.gqnk9 8snf.