Effects of Climate Change on the Distribution of Akebia quinata

Akebia quinata, also known as chocolate vine, is a creeping woody vine which is used as Chinese herbal medicine, and found widely distributed in East Asia. At present, its wild resources are being constantly destroyed. This study aims to provide a theoretical basis for the resource protection of this plant species by analyzing the possible changes in its geographic distribution pattern and its response to climate factors. It is the first time maximum entropy modeling (MaxEnt) and ArcGIS software have been used to predict the distribution of A. quinata in the past, the present, and the future (four greenhouse gas emission scenarios, namely, SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5). Through the prediction results, the impact of climate change on the distribution of A. quinata and the response of A. quinata to climate factors were analyzed. The results showed that the most significant climatic factor affecting the distribution pattern of A. quinata was the annual precipitation. At present, the suitable distribution regions of A. quinata are mainly in the temperate zone, and a few suitable distribution regions are in the tropical zone. The medium and high suitable regions are mainly located in East Asia, accounting for 51.1 and 81.7% of the worldwide medium and high suitable regions, respectively. The migration of the geometric center of the distribution regions of A. quinata in East Asia is mainly affected by the change of distribution regions in China, and the average migration rate of the geometric center in each climate scenario is positively correlated with the level of greenhouse gas emission scenario.


INTRODUCTION
Akebia quinata, also known as chocolate vine, is a creeping woody vine widely distributed in East Asia (Wang et al., 2021). As a valuable Chinese herbal medicine, its fruit and stem have analgesic, diuretic, and anti-inflammatory effects (Park et al., 2018;Song et al., 2018). Over the years, the wild resources of this species have been destroyed due to increased cutting, and this species is facing more and more severe survival challenges. The protection of A. quinata is very urgent, and some reasonable and efficient protection schemes need to be implemented. This requires us to understand the suitable growth conditions of A. quinata and understand its geographical distribution and the impact of climate change upon it. At present, this species has only been reported in the fields of pharmacological activity (Sung et al., 2015;Lee et al., 2017), chemical composition (Jiang et al., 2006;Mimaki et al., 2007), and genome research (Li et al., 2016), but there are still gaps in knowledge around its geographical distribution and climate response.
The geographic distribution data of species are the basic information for studying their life process. In recent years, due to human activities and climate change, the habitat of many species has been destroyed and changed, and the survival of species has been threatened, or some even driven to extinction. Accurate simulation and prediction of species' distribution area is the key to their conservation. To achieve this goal, many factors need to be considered, especially climatic factors. On the one hand, climatic factors play a decisive role in the life process of species, and the accelerated prosperity or extinction of species can even depend on it (Lenoir et al., 2008;Acevedo et al., 2020). On the other hand, climate change also has a profound impact on the distribution pattern of species (Moraitis et al., 2019;Wang et al., 2019;Wilson et al., 2019) and will lead to a series of changes in the survival rate, dominance, community structure, and replacement rate of species (Yi et al., 2018). Therefore, it is more important to make clear the distribution of species and predict the change of their distribution in the future. Accurate distribution data can not only be used to predict the various possibilities of species distribution change but are also critical to assess the potential impact of changing ecosystems.
Climate is one of the main determinants delimiting the geographical distribution of plant species on large scale (Ferrarini et al., 2019). There is a considerable amount of research declaring that climate change leads to expansion or retraction in plant species' ranges (Thuiller et al., 2005;Ferrarini et al., 2018). To assess the vulnerability of plant species under a rapidly changing climate, we can use species distribution modeling (SDM) to predict species climate niches and project their potential future range shifts (Huntley et al., 1995;Pearson and Dawson, 2003;Thuiller et al., 2005;Alavi et al., 2019).
Maximum entropy modeling (MaxEnt) is a very powerful presence/pseudo-absence algorithm. Many authors have suggested that it is one of the most efficient approaches for predicting the potential distributions of species (Elith et al., 2006(Elith et al., , 2011Phillips et al., 2006a). The model can simulate and predict the potential geographical distribution of species by using the current distribution data and various environmental parameters (Phillips et al., 2006a;Phillips and Dudík, 2008). It has the advantages of small sample size, fast running speed, and stable operation (Phillips et al., 2006a;Estes et al., 2013;Li J. et al., 2020). Even in the case of insufficient species distribution information, it also has good accuracy and can test the accuracy of prediction results (Saatchi et al., 2008;Yi et al., 2017). Therefore, it is widely used in many aspects of species distribution analysis (Yang et al., 2013;Qin et al., 2017;Zhang et al., 2019).
In this study, the MaxEnt model was used to simulate and predict the distribution of A. quinata in different climatic scenarios. The purpose of this study was (1) to find the key climatic factors that restrict the distribution of A. quinata; (2) to predict the distribution pattern of A. quinata under different climate scenarios; (3) to evaluate the impact of climate change on the distribution pattern of A. quinata; and (4) to predict the concentrated distribution region of A. quinata, analyze the migration of its concentrated distribution region, and explore the migration reasons of its concentrated distribution region.

Location Data Sources of A. quinata
The geographic data on A. quinata distribution were collected from two sources. Data were collected within China using a GPS device (GARMIN GPSMAP 63SC, Kansas City, KS, USA) from the field survey in 2017-2019, covering Northwest China, Central China, South China, and East China. Data from other parts of the world were collected from the website of the Global Biodiversity Information Facility (GBIF, https://www.gbif.org). Based on the above 1,024 effective records (Figure 1), the prediction model was established. The actual distribution of A. quinata was analyzed using ArcGIS (version 10.2, ESRI, Redlands, CA, USA) software. The geographic distribution map of richness ×1 was drawn using a geographic information system (GIS).
We used two methods to filter the data downloaded using GIBF. First, we deleted the data with missing longitude and latitude information and fuzzy geographic location information and, second, we deleted the data that did not conform to the geographic coordinate system (World Geodetic System 1984) used in this study.

Climate Data Sources
It is reported that 19 bioclimatic variables (Table 1) are the most typical and important variables to form the potential species distribution model (Molloy et al., 2013;Li J. et al., 2020). In this study, the current , past, and future climate scenarios were downloaded from WorldClim Version 2.1 (this version was released in January 2020) (www.worldclim.org) (Fick and Hijmans, 2017). These data include 19 climate variables. All variables were cast to USA Contiguous Albers Equal Area Conic (NAD 1983) and resampled using nearest-neighbor to a 30 arcsecond resolution using ArcGIS. Using ArcGIS conversion tools, the environment variables were converted to ASCII format.
This study used the scenarios of the Last Interglacial, the Last Glacial Maximum, and the Mid-Holocene to predict the species distribution in the past. These three scenarios were provided by the Fourth Phase of Community Climate System Model (CCSM4) in the Fifth Phase of Coupled Model Intercomparison Project (CMIP5).

Prediction Using the MaxEnt Model
The occurrence data of A. quinata and climate data were input into the MaxEnt (MaxEnt 3.4.1) model (Phillips et al., 2006b). We randomly selected 25% of the point data as the test set and 75% of the point data as the training set. The model was run 10 times to evaluate the average results. The area under curve (AUC) of the receiver operating characteristic (ROC) curve was selected to evaluate the performance of the model. The ROC curve analysis is a method to verify the performance of the MaxEnt model. Its principle is to judge the prediction accuracy of the model by calculating the area enclosed by the curve and abscissa, that is, the AUC value. Generally, the model has five accuracy levels. When the AUC value is between 0.50 and 0.60, the prediction fails. The prediction accuracy between 0.60 and 0.70 is poor, the prediction effect between 0.70 and 0.80 is average, the prediction effect between 0.80 and 0.90 is good, and the prediction accuracy between 0.90 and 1.00 is excellent . The response curve of climate variables generated by the model reflects the relationship between the value of climate variables and the existence probability of A. quinata. A jackknife test and statistical table of contribution rate of climate variables were used to screen the climate variables with high importance. Too many climate variables will reduce the computational efficiency of the MaxEnt model in predicting the distribution of species on a large regional scale, and the climate variables with obvious collinearity will affect the prediction accuracy Sillero and Barbosa, 2020). Therefore, after running the MaxEnt model with 19 climate variables, this study screened climate variables and compiled the model again for prediction. The screening process was as follows: (1) The Pearson correlation coefficient in SPSS (Statistical Product and Service Solutions, version 26.0, Armonk, NY, USA) software was used to calculate the correlation between climate variables. (2) We removed all variables (variables whose percent contribution is <1%) whose percent contribution in the model prediction is lower than the contribution threshold setting. Next, among the variables with high correlation (the absolute value of the correlation coefficient is greater than or equal to 0.8), the variable with the highest contribution rate was retained, and other variables were removed.

Suitability Division of the Distribution Regions
By using the To Raster tool in ArcGIS, the American Standard Code for Information Interchange (ASCII) file exported using the MaxEnt model converted the raster into a grid layer, and the suitable distribution regions of species were obtained. In the predicted suitable distribution regions, the existence probability of species was between 0 and 1. According to the actual distribution and field survey results, using the Reclassify tool of ArcGIS and artificial classification method, the prediction results were divided into four grades: high suitability (>0.66), medium suitability (0.33-0.66), low suitability (0.15-0.33), and no suitability (<0.15). This study used the world climate data to run the MaxEnt model. The map data and results of specific regions in this study were extracted from the global prediction results.

Geometric Center Analysis of the Distribution Regions
Raster Calculator tool in ArcGIS was used to separate the grid layer of species distribution, and the suitable distribution regions were separated. Next, the Raster Domain tool was used to convert the grid layer of the suitable distribution regions into a face, and then the Mean Center tool was used to find out the geometric center of the suitable distribution regions.

Model Accuracy Evaluation
In this study, the average AUC of the MaxEnt model is 0.956 (Figure 2A), which indicates that the prediction accuracy is excellent, and the results can be used.

Important Climate Variables
Supplementary Table S1 shows the contribution rates of climate variables in the MaxEnt model. Supplementary Table S2 shows the correlation test between 19 climate variables. According to the screening principle, eight climate variables (bio1, bio2, bio6, bio9, bio12, bio14, bio15, and bio18) are retained for recompiling the MaxEnt model operation. Figure 2B shows the Jackknife test (using AUC on test data). When using a single variable, the climate variable with the highest gain is bio12 (annual precipitation), and its gain value is >0.93. In addition, bio1 (annual mean temperature) and bio6 (minimum temperature of the coldest month) are the two most gain variables after bio12. They are three important climatic variables that restrict the geographic distribution of A. quinata. Figures 3A-C shows the response curves between the above three climate variables and the probability of the existence of A. quinata.
The probability of the existence of A. quinata is close to 0 when the annual precipitation (the most significant variable, bio12) is <570 mm, then increases rapidly and reaches the maximum when bio12 is 2,300 mm. The probability of existence decreases when bio12 continues to increase, but when bio12 is more than 2,700 mm, the probability of existence does not continue to decrease and remains in the low suitability range (0.15-0.33). There is no clear upper limit but there exists a clear lower limit of annual precipitation for the suitable distribution regions of A. quinata. According to the division of suitability, to meet the minimum existence probability (>0.15) of A. quinata, at least bio12 should be >1,000 mm. To achieve medium suitable survival conditions for A. quinata (probability of existence >0.33), bio12 should be in the range of 1,160-2,700 mm. The performance of the temperature factor in the suitable distribution regions of A. quinata is different from precipitation, and the temperature factor has clear upper and lower limits. When the temperature of bio1 (annual mean temperature) is <5.0 • C, the probability of existence is close to 0. Due to the continuous increase of bio1, the probability of existence increases rapidly. When the temperature of bio1 increases to 15.5 • C, the probability of existence reaches the maximum and is close to 0.63, then decreases rapidly and finally decreases to 0. The response curve of bio6 (minimum temperature of coldest month) is similar to bio1, except that the temperature threshold of existence probability is different. The two thresholds of existence probability approaching 0 are −14.0 and 19.0 • C, respectively, and −3.0 and 0.0 • C is the temperature range with the highest probability of existence (probability of existence is close to 0.60).

Potential Distribution Regions of A. quinata
According to the prediction results of the MaxEnt model, the prediction of the current and the future distribution is relatively clear, but the distribution of A. quinata is not found in the past three periods (Supplementary Figure S1).

Current Potential Distribution
According to the current prediction results (Figure 4A), most of the suitable distribution regions of A. quinata were mainly in the temperate zone, and a small number of suitable distribution regions were in the tropical zone. The total suitable area was 592.87 × 10 4 km 2 , the low suitable area was 447.35 × 10 4 km 2 , the medium suitable area was 143.99 × 10 4 km 2 , and the highly suitable area was 1.53 × 10 4 km 2 . In East Asia, the occurrence data of A. quinata were the most intensive. The suitable distribution regions were also mainly located in East Asia (Figure 4B), and the areas of medium and high suitability were 73.55 × 10 4 and 1.25 × 10 4 km 2 , respectively, accounting for 51.1 and 81.7% of the global medium and high suitable area, respectively. East Asia is the region with the most concentrated distribution of A. quinata, and so is a region worthy of attention in this study.

Supplementary Table S3
shows the suitable area of A. quinata in the world under different climate scenarios. Figure 5A shows the suitable area change of A. quinata in the form of a broken line diagram. Worldwide, regardless of the transition from the current scenario to any scenario in 2021, the suitable area of A. quinata shows an increasing trend. In the SSP1-2.6 scenario, the suitable area shows a downward trend after 2021 and an upward trend from 2041. In the SSP2-4.5 scenario, the suitable area shows an upward trend before 2041 and begins to decline after 2041. In the SSP3-7.0 and SSP5-8.5 scenarios, the suitable area increased significantly from 2041 to 2061 and then remained stable. In the first two scenarios, the suitable area remains relatively stable, while in the latter two scenarios, the suitable area shows a significantly increasing trend, especially in the SSP5-8.5 scenario. Supplementary Figures S2-S5 shows the prediction picture of the suitable distribution regions of A. quinata in the world. A remarkable phenomenon is that from 2061 to 2081, with the upgrading of the greenhouse gas scenario, the suitable area for A. quinata in Europe will expand significantly, and there is a trend to expand to the northeast. With the upgrading of climate scenario, the suitable area of A. quinata in North America increases and tends to expand to the north, while the suitable area in South America shows a downward trend. In the rest of the regions except East Asia, the change in the suitable area of A. quinata is not obvious.
In East Asia, the suitable area of A. quinata shows a continuous decreasing trend under the SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios, especially in SSP5-8.5 (Figure 5B). In the SSP1-2.6 scenario, the suitable area shows a downward trend before 2041, an upward trend in 2041 and 2061, and then continues to decline. A. quinata is mainly distributed in China, Korea, and Japan in East Asia. Figures 6A-H shows the suitable distribution regions of A. quinata in East Asia under the SSP1-2.6 and SSP5-8.5 scenarios. It is obvious that under the SSP5-8.5 scenario, the suitable distribution region of A. quinata in China shows a trend of continuous fragmentation from 2041 to 2081. Under the SSP1-2.6 scenario, the suitable distribution region in China can remain in a relatively stable state. From 2061 to 2081, under the SSP2-4.5 and SSP3-7.0 scenarios, the suitable distribution region in China also showed a fragmentation trend (Supplementary Figure S6). In Japan and South Korea, the suitable distribution regions of A. quinata remain stable in all scenarios. It can be seen from the line chart of suitable area that the suitable area of A. quinata in East Asia is closely related to that of China (Figures 5B,C). Supplementary Tables S4, S5 show the specific suitable area values in East Asia and China, respectively.

Geometric Center of Suitable Distribution Regions and Its Migration
The overall change of the suitable area can be expressed by the shift of the geometric center of the suitable distribution regions of A. quinata in East Asia. Based on the predicted potential distribution, the geometric centers of the distribution regions under different climate scenarios are obtained (Figure 7). Figure 7A shows the geometric centers under the current climate scenario and the 16 different future scenarios. It can be seen intuitively that these geometric centers generally show a trend of migration to the northeast. The migration distance of the geometric center is different in four different climate situations. In the SSP1-2.6 scenario (Figure 7B), the migration distance of the geometric center is relatively small, in the SSP2-4.5 and SSP3-7.0 scenarios (Figures 7C,D), the migration distance of the geometric center is relatively medium, and in the SSP5-8.5 scenario (Figure 7E), the migration distance of geometric center is relatively large.
To express them more clearly, the migration of the geometric center are quantified. Considering the geometric center under the current climate scenario as the origin, the migration rate of the geometric center farthest from the origin is expressed as 1. The ratio of the distance between other geometric centers and the origin to the farthest distance is the migration rate of these geometric centers ( Table 2). The climate scenario with the largest migration rate is SSP5-8.5 in 2061 and 2081, and the migration rate is 1. The climate scenario with the lowest migration rate is SSP3-7.0 in 2021, and the migration rate is 0.17. The average migration rates of the four climate scenarios are 0.29 (SSP1-2.6), 0.49 (SSP2-4.5), 0.55 (SSP3-7.0), and 0.69 (SSP5-8.5), respectively.
In the SSP1-2.6 scenario, the migration rate of the geometric center increases slightly from 2021 to 2041, decreases in 2061 (the geometric center makes a return motion), and increases again in 2081, which is in a relatively stable dynamic equilibrium (the migration rate fluctuates between 0.22 and 0.37). It can be speculated that the distribution area of A. quinata may continue to maintain a relatively stable state in this scenario. In other scenarios, the geometric center does not make a return motion. And, with the continuous upgrading of the scene, the average migration rate of the geometric center also increases.
The decrease of the suitable area of A. quinata in East Asia is mainly affected by the change of the suitable area in China, while the suitable area of Japan and Korea has almost no change. Therefore, the weight of the suitable area of Japan and Korea in the total suitable area of East Asia has increased. The change of the suitable area weight causes the geometric center of the distribution area to move to the northeast, and the average mobility of geometric centers under different climate scenarios is positively correlated with the level of greenhouse gas emission scenarios.

DISCUSSION
The change of plant distribution patterns is different under climate warming. The related research on the prediction of Cunninghamia lanceolata distribution shows that an increase in greenhouse gas emissions may lead to the decrease of the suitable area of C. lanceolata . In the study of two species of peony (Zhang et al., 2018), the suitable areas of Paeonia delavayi and Paeonia rockii will increase under the low concentration greenhouse gas emission scenario (RCP2.6), but the suitable area of P. rockii will increase and the suitable area of P. delavayi will decrease under the high concentration greenhouse gas emission scenario (RCP8.5). According to the related research of Coptis herbs, the suitable areas of Coptis chinensis and Coptis teeta will decrease, and the suitable area of Coptis deltoidea will increase in the future RCP8.5 scenario (Li J. et al., 2020). The prediction in this study shows that the suitable area of A. quinata in different regions of the world changes differently. In East Asia, when transitioning from the current scenario to three greenhouse gas emission scenarios (SSP2-4.5, SSP3-7.0, and SSP5-8.5), the suitable area of A. quinata will be significantly reduced. Compared with other scenarios, in the low concentration greenhouse gas emission scenario (SSP1-2.6), the suitable area change of A. quinata in East Asia is more conservative.
Global climate change will not only cause temperature changes in different regions but also change the distribution pattern of precipitation, resulting in changes in the distribution of A. quinata. Generally, plants can adapt to the fluctuation of climate factors within a certain threshold range, but when the change of climate factors approaches or even exceeds the threshold range, it will lead to the migration of their distribution (Camille and Gary, 2003). Plants need enough water to grow, but drought will limit their growth. When the precipitation in the driest month increases, it helps to prolong the growing season of the plants and promote their migration to more suitable habitats (Vaganov et al., 1999). In addition, extreme temperatures also significantly affect the growth of plants. If the minimum temperature of the coldest month drops, it will undoubtedly aggravate freezing and chilling injuries and cause plant death (Camille and Gary, 2003). The increase of maximum temperature in the warmest month may destroy the water balance in plants and hinder their metabolic function (Lemmens et al., 2006). The changes in these climate factors are directly reflected in the increase or decrease of suitable area. If the climate change is too large, it will cause more serious changes, that is, habitat fragmentation.

CONCLUSION
It is of great significance to predict the distribution pattern of A. quinata in different climatic conditions and analyze the response relationship between A. quinata and climatic factors for its protection and research. The results show that the concentrated distribution region of A. quinata is in East Asia. And, bio1 (annual mean temperature), bio6 (minimum temperature of the coldest month), and bio12 (annual precipitation) are the main climatic factors affecting the distribution pattern of A. quinata. In East Asia, when transitioningfrom the current scenario to three greenhouse gas emission scenarios (SSP2-4.5, SSP3-7.0, and SSP5-8.5), the suitable area of A. quinata will be significantly reduced. Compared with other scenarios, under the low concentration greenhouse gas emission scenario (SSP1-2.6), the change in suitable area of A. quinata in East Asia is more conservative. The geometric center of the distribution area of A. quinata in East Asia will move to the Northeast under the climate warming, which is mainly due to the decrease of the distribution area of A. quinata in China. And, the average migration rate of the geometric center under each climate scenario is positively correlated with the level of greenhouse gas emission scenario.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
WW conceived and X-QX designed the study. J-MZ and WW processed the data, performed the analyses and analyzed the results, and wrote the manuscript. M-LS, Z-JL, X-YP, SS, and BL edited the manuscript. All authors read and approved the final version of the manuscript.