Predicting Possible Distribution of Tea (Camellia sinensis L.) under Climate Change Scenarios Using MaxEnt Model in China

Climate change has dramatic impacts on the growth and the geographical distribution of tea (Camellia sinensis L.). Assessing the potential distribution of tea will help decision makers to formulate appropriate adaptation measures to use the altered climatic resources and avoid the damage from climate hazards. The objective in this study is to model the current and future distribution of tea species based on the four SSPs scenarios using the MaxEnt model in China. For the modeling procedure, tea growth records in 410 sites and 9 climate variables were used in this paper. The area under the receiver operating characteristic (ROC) curve (AUC) was used to evaluate the performance of the model. The AUC value was over 0.9 in this study, showing the excellent simulation result of the model. In relation to the current distribution, areas of 82.01 × 104 km2 (8.51% of total land area in China), 115.97 × 104 km2 (12.03% of total land area in China), and 67.14 × 104 km2 (6.97% of total land area in China) were recognized as Marginal, Medium, and Optimal climate suitable habitats for tea over China. Compared to the current distribution, most of the Optimal suitability areas in southeast China would be lost in four scenarios. The area of Marginal and Medium suitable habitats would expand in SSP370 and SSP585, especially in 2041–2061 and 2081–2100. The suitable area of tea would expand northwards and westwards, suggesting that additional new suitable habitats could be created for tea production with the future climate change, especially in Shandong, Henan, Guizhou, and Yunnan Provinces. This research would provide vital scientific understanding for policy making on tea production, tea garden site chosen and adopyion of adaptation methods in the future.


Introduction
Global climate change is happening at an unprecedented speed, which is partially accelerated by human activities. The evidence is shown from the IPCC (Intergovernmental Panel on Climate Change) report that temperature of global surface in 2011-2020 was 1.09 • C higher than 1850-1900, with larger increases over land (1.59 • C) than over the ocean (0.88 • C) [1]. In China, the annual average temperature increased at a rate of about 0.15 • C/10 years from 1901 to 2020, and the warming trend was more obvious in 1951-2020 with the rate of 0.26 • C/10 years. The past 20 years have been the warmest period since 1900, which had a profound impact on biodiversity, species distribution, community composition, vegetation pattern, ecosystem structure, function, and service [2].
Tea (Camellia sinensis L.) is widely distributed in tropical and subtropical mountainous areas over 50 countries and regions all over the world [3]. In China, the long history of tea cultivation and suitable environment provide favorable conditions for tea production, the yield and area of tea have been rising continuously and the tea industry has been booming in recent years [4]. As a green cash crop, the successful cultivation of tea is highly dependent on climate conditions, including moderate temperature, rainfall, solar radiation and so on [5][6][7]. Han et al. [8] found that tea yield will decrease by 11% to 35% due to extreme weather events, which may cause the fluctuations of tea prices. Meanwhile, the suitable area can be influenced in some provinces by climate warming. Jin et al. [9] demonstrated that the climate suitability of tea decreased significantly in the 1970s and increased obviously after the 1990s in Zhejiang Province because of the climate change. From 1981 to 2013, the suitability of main tea producing areas in Fujian Province showed a downward trend, and the decline in coastal counties was greater than that in mountainous counties [10]. On the contrary, tea planting suitability increased in the north production region, leading to the northern movement of tea production boundary in China due to the warming [4]. From 1961 to 2019, the limit of the suitable tea plant area in the eastern provinces in China moved obviously to north mountainous areas [11]. At present, the research on tea suitable areas is scattered, and the research contents are mostly based on the previous climate conditions. Future climate change may have a negative impact on the distribution areas of tea in China, but the research on the suitable regional distribution of tea at a national perspective is still relatively lacking.
In the future, the productivity and sustainability of tea gardens could be affected by the rising temperature and decreasing rainfall [12]. It has been found that in Sri Lanka, the monthly rainfall reduction of 100 mm would reduce productivity of 30-80 kg "made" tea/ha/month [13], and the change of rainfall also could reduce the labor demand of the whole tea garden sector in this country by about 1,175,000 man-days per year by 2050, which is likely to have considerable social and welfare implications [14]. However, the existing studies seldom investigate the changes of tea species niche in China. The variation of priority areas suitable for tea industry development under climate change is poorly understood. It is necessary to understand the distribution characteristics change of tea species under the climate change and the impact of climate performance on the potential habitats of tea, which could provide support for tea production planning and management. In modeling suitable distribution, species distribution model (SDMs) is often used to simulate species distribution and its ecological requirements based on the existing distribution data and related environmental factors [15]. Over the past decades, numerous SDMs from a range of methods, such as domain environmental envelope (DOMAIN), Ecological Niche Factor Analysis (ENFA) [16], Genetic Algorithm for Rule-set Production (GARP) [17], Artificial Neural Network (ANN) [18], and Maximum Entropy (Maxent), have been developed. The Maximum Entropy (Maxent) model is a density estimation and species distribution model, and it is a selective method based on the maximum entropy theory [19]. This was developed in 2006 and had played critical roles in species reserves planning, prediction of the potential distribution of invasive species, and spatial distribution of plant and animal species [20][21][22]. The model is robust in prediction capability and stability, and it is not affected by sample size and sample deviations [23]. Therefore, the MaxEnt model was chosen to simulate and analyze the distribution of suitable areas of tea species in the future at a national level in China based on new SSPs climate scenario data instead of RCPs. The results are expected to enhance scientific understanding of variation of priority areas suitable for the development of tea industry under climate change.
Studying the current and potential distribution of tea on a large scale can help to understand the dynamic of overall planting patterns under climate change, which could be used to determine the suitable area for tea cultivation in China. The objective of this study was to use the MaxEnt model to map the current potential distribution of tea (Camellia sinensis L.) by the key environmental variables that highly correlate with tea distribution, predict the potential distribution of tea species under future climate change scenarios, and evaluate the effects of climate change on the suitable habitat area of tea, all to identify the change suitable area under future climate change scenarios.

Species Data of Camellia sinensis
The tea occurrence data in Chinese mainland were collected from databases of Chinese Virtual Herbarium (CVH, http://www.cvh.ac.cn, last accessed on 20 October 2020) and the Global Biodiversity Information Facility (GBIF, http://www.gbif.org, last accessed Agriculture 2021, 11, 1122 3 of 18 on 20 October 2020). These two websites provide the most accurate and reliable location information of tea species. In this study, a total of 3921 occurrences data were obtained, out of which 2030 were collected from CVH and 1891 were collected from GBIF. Firstly, the GPS and Coordinate pick up system of Baidu Map (http://api.map.baidu.com/lbsapi/ getpoint/index.html) were used to carry out geographic registration on the event records with detailed location information. Then, 2854 incomplete records as well as 657 duplicated entries were removed. Finally, 410 accurate existence points remained (Figure 1). The coordinates of latitude and longitude for each record were stored in an Excel database as ".csv" format with the title "species, longitude, latitude" for MaxEnt model building.

Species Data of Camellia sinensis
The tea occurrence data in Chinese mainland were collected from databases of Chinese Virtual Herbarium (CVH, http://www.cvh.ac.cn, last accessed on 20 October 2020) and the Global Biodiversity Information Facility (GBIF, http://www.gbif.org, last accessed on 20 October 2020). These two websites provide the most accurate and reliable location information of tea species. In this study, a total of 3921 occurrences data were obtained, out of which 2030 were collected from CVH and 1891 were collected from GBIF. Firstly, the GPS and Coordinate pick up system of Baidu Map (http://api.map.baidu.com/lbsapi/getpoint/index.html) were used to carry out geographic registration on the event records with detailed location information. Then, 2854 incomplete records as well as 657 duplicated entries were removed. Finally, 410 accurate existence points remained (Figure 1). The coordinates of latitude and longitude for each record were stored in an Excel database as ".csv" format with the title "species, longitude, latitude" for MaxEnt model building.

Bioclimatic Data
Bioclimatic variables are the key factors that determine the distribution of tea varieties [24]. Nineteen bioclimatic variables from WorldClim 2.1 (https://www.worldclim.org/data/cmip6/cmip6_clim2.5m.html, accessed on 20 October 2020) were extracted as environmental data [25]. The new version of 19 bioclimatic variables was released in January 2020 and the data records are from 1970 to 2000 period, which have been widely used to predict the current potential distribution of species in recent research due to the ecological meaning and can reflect the seasonal variation characteristics of temperature and precipitation [26,27]. WorldClim could provide four gridded global climate layers at 10 min, 5 min, 2.5 min, and 30 s in CGS_WGS_1984 projection [28]. Due to the diversity of climatic condition and topography in China in the short term, we chose the spatial

Bioclimatic Data
Bioclimatic variables are the key factors that determine the distribution of tea varieties [24]. Nineteen bioclimatic variables from WorldClim 2.1 (https://www.worldclim. org/data/cmip6/cmip6_clim2.5m.html, accessed on 20 October 2020) were extracted as environmental data [25]. The new version of 19 bioclimatic variables was released in January 2020 and the data records are from 1970 to 2000 period, which have been widely used to predict the current potential distribution of species in recent research due to the ecological meaning and can reflect the seasonal variation characteristics of temperature and precipitation [26,27]. WorldClim could provide four gridded global climate layers at 10 min, 5 min, 2.5 min, and 30 s in CGS_WGS_1984 projection [28]. Due to the diversity of climatic condition and topography in China in the short term, we chose the spatial resolution of 2.5 min (5 km 2 ) for modeling. The datasets downloaded from WorldClim are in "tiff" format and were converted to ASCII files in ArcGIS by "Raster to ASCII", for use in MaxEnt software.

Preprocessing of Bioclimatic Variables
Multicollinearity among variables may lead to over-fitting of the model and misinterpretation of the modeled results [29]. In order to avoid the influence of the cross-correlation of the 19 bioclimatic variables, the contribution rate of each variable was firstly evaluated by the Jackknife test in MaxEnt software version 3.4.1 (https://biodiversityinformatics. amnh.org/open_source/maxent/, accessed on 20 October 2020) with default setting, and then correlation analysis was carried by SPSS software version 22.0 (SPSS, Inc., Chicago, IL, USA) to eliminate the highly correlated variables [30]. If the absolute value of correlation coefficient of two variables was greater than 0.8, only one variable with a greater contribution was selected [31]. At last, the remaining 9 bioclimatic variables (bio2, bio3, bio6, bio7, bio10, bio12, bio14, bio15, and bio18) were selected to simulate the distribution of tea in China under current and future climatic conditions (Table 1, in bold). The change of land use and land cover, human interference, species diffusion and biological interactions were neglected in the modeling. Precipitation of Coldest Quarter mm The bolded variables were used in modeling.

Species Distribution Modeling Procedure
In this study, the MaxEnt software [19] was used to model the current and future distribution of tea in China because it has better prediction capabilities and stability that will not be affected by small sample deviations and fewer samples [23,32,33]. The Species data were inputted in the box of Samples, and the Bioclimatic variables were inputted in the box of Environmental layers. The MaxEnt model was run with 10 replicates and 25% of random test percentage, the convergence threshold of 10 −5 , a maximum number of backgrounds points of 10,000, the maximum iterations as 1000, and a regularization parameter value of 0.5, an auto-feature option, and cloglog output format. Meanwhile, "Create response cures" and "Do Jackknife to measure variable importance" were chosen before running the model. Other values were kept as default [34].
The random test percentage was 25%, which means that 75% of the total database was used as the random sample to train the model, and other 25% of the total database was used to test the model predictions [35]. To avoid over-fitting of the test data, we set the regularization multiplier value as 0.5. The output format was default "cloglog" transform, which has the stronger theoretical justification than the logistic transform [36]. The auto features setting was used in this research based on the number of sample sizes. Phillips and Dudík [37] classified the features for sample sizes: auto features setting for more than 80 records, quadratic and hinge setting for 15-79 record, linear and quadratic for 10-14 record, and linear setting for sample sizes fewer than 10 records [24,38]. The Jackknife test (systematically leaving out each variable) was used to assess the contribution and the importance of each bioclimatic variables to the model [39].

Evluating the Model Performance
The performance of this model was evaluated based on the computed receiver operating characteristic (ROC) curve and the area under the curve (AUC) [19,23]. Taking the false positive rate (the probability of being predicted as positive without the actual distribution of the species) as the horizontal axis and the true positive rate (the probability of being predicted as positive by species) as the vertical axis, the ROC curve was drawn [40]. The area surrounded by the curve and the horizontal axis was AUC (area under curve) value. In general, AUC values vary from 0.5 to 1, which could be divided into five classes: fail (0.5-0.6), poor (0.6-0.7), fair (0.7-0.8), good (0.8-0.9), and excellent (0.9-1) [41]. The closer the AUC value was to 1, the farther away from the random distribution, the greater the correlation between environmental variables and the predicted geographical distribution of species, and the more accurate the performance of the model, while AUC < 0.5 is a contingency difference, which can be regarded as a stochastic forecasting model, and rarely happens [42].

Future Spatial Change Detecation
The same nine bioclimatic variables were employed to analyze the future potential distribution of tea from 2040 to 2100. The global climate model, BCC-CSM2-MR (modeling data from the Beijing Climate Center Climate System Model) was used as future simulated climate at four emissions scenarios driven by different socio-economic assumptions named as Shared Socio-economic Pathways (SSPs) 126, 245, 370, and 585 with a 2.5 min spatial resolution, which were downloaded from WorldClim (http://www.ccafs-climate.org) and both available for three time-periods including 2050s (2041-2060), 2070s (2061-2080) and long-term 2090s (2081-2100) [43,44], which can better simulate the temperature and precipitation variables in Asia [45,46]. SSP126 represents the low end of the forced pathways in the future, updating the RCP2.6, and the radiative forcing in 2100 is 2.6 W/m 2 . SSP245 represents the medium part of the range of future forcing pathways and updates the RCP4.5 pathway. The radiation forcer will be stable at 4.5 W/m 2 in 2100. SSP370 represents the medium to high end of the range of future forcing pathways, which the radiation forcing will stabilize at 7.0 W/m 2 in 2100. SSP585 represents the high end of the range of future pathways, updating the RCP8.5 pathways, and its emissions are high enough to generate a radiative forcing of 8.5 W/m 2 in 2100 [47].
After modeling the species distribution, the spatial changes were calculated for the future compared to the current by in ArcGIS 10 software. The prediction results of current and future were imported into ArcGIS by ASCII to raster" and reclassified as Unsuitable habitat (0-5%), Marginal suitable habitat (5-33%), Medium suitable habitat (33-67%), and Optimal suitable habitat (67-100%) [48,49]. In order to compare the future and current area variation and distributions, Unsuitable habitat, Marginal suitable habitat, Medium suitable habitat, and Optimal suitable habitat was further reclassified as 0, 1, 2, 3 and 4, and the current results layer was subtracted from the future results by the "subtraction" tool in Arctoolbox [50]. Area calculation of the suitable habitat and the subtraction result were carried out in ArcGIS 10.

Model Performance and Variable Contributions
The average AUC value of the test data and training data in this study was 0.902 ( Figure 2), indicating that this model performs excellent in modeling tea distribution in China. It is worth noting that the AUC values tend to be lower for species with wide distribution range [51,52].
Agriculture 2021, 11, 1122 6 of 18 mation by themselves than the other variables did, and more sensitive to Camellia sinensis habitats suitability. Precipitation of Driest Month (bio14), Precipitation of Warmest Quarter (bio18), Mean Diurnal Range (bio2), and Temperature Annual Range (bio7) had moderate gain when used independently. Other variables such as Mean Temperature of Warmest Quarter (bio10), Precipitation Seasonality (bio15), and Isothermality (bio3), had low gains with few information when used in isolation. The importance of the permutation of the predictor variables to the MaxEnt model are shown in Table 2. Among the input environmental variables, Temperature Annual Range (bio7) was most influential and contributed 51% to this model. Precipitation of Warmest Quarter (bio18) was the second most important variable with the permutation importance of 14%, followed by Min Temperature of Coldest Month (bio6) (9.7%) and Mean Temperature of Warmest Quarter (bio10) (8.9%). Precipitation of Driest Month (bio14) and Mean Diurnal Range (bio2) accounted for 4.9% and 4.2%. Annual Precipitation (bio12) and Precipitation Seasonality (bio15) had the smaller influence on the habitat model, which was 3.7% and 2.9%, respectively. Isothermality (bio3) only contributed 0.7%, which was the lowest contribution, as shown in Table 2. Meanwhile, the percent contributions of the environmental variables showed that Precipitation of Driest Month (bio14) had the greatest impact of 63.9% followed by Annual Precipitation (bio12) of 18.9%, while Temperature Annual Range (bio7) contributed 4.5%, and the lowest was Precipitation of Warmest Quarter (bio18), 0.3%. An individual variable with highest training gain in Jackknife test means that it was the most significant bioclimatic variable with more contribution or information to species habitat distributions [53]. The results showed that the Annual Precipitation (bio12) and Min Temperature of Coldest Month (bio6) provided very high gains (>1.2) when used independently (Figure 3), indicating that bio12 and bio6 contained more useful information by themselves than the other variables did, and more sensitive to Camellia sinensis habitats

Variable Response Curves
The quantitative relationship between environmental variables and the logical probability of existence (also known as habitat suitability) is shown as the response curves in Figure 4, which could deepen the understanding of the species niche. Based on the results of the Jackknife test, each environmental factor was selected for modeling to reveal the climate characteristics' impacts on tea suitable distribution area. According to the response curve of Temperature Annual Range (bio7), the probability of tea occurrence was lower than 26 °C (Figure 4a). The suitable range of Precipitation of Warmest Quarter (bio18) was from 500 to 900mm with a peak at 600mm (Figure 4b). The response curve of Min Temperature of Coldest Month (bio6) showed that habitat suitability of tea is associated with areas where the optimal temperature ranged over 0 °C , but it has a fluctuant range at 10 to 16 °C (Figure 4c). The response curve of Mean Temperature of Warmest Quarter (bio10) indicated that 27 °C would be conducive for the development of tea in Warmest Quarter (Figure 4d). The response curve of Precipitation of Driest Month (bio14) demonstrated that there was a wide range of precipitation from 17 to 180 mm to adapt to the tea habitat, and it reaches the peak at 48 mm (Figure 4e). The highest suitability of

Variable Response Curves
The quantitative relationship between environmental variables and the logical probability of existence (also known as habitat suitability) is shown as the response curves in Figure 4, which could deepen the understanding of the species niche. Based on the results of the Jackknife test, each environmental factor was selected for modeling to reveal the climate characteristics' impacts on tea suitable distribution area. According to the response curve of Temperature Annual Range (bio7), the probability of tea occurrence was lower than 26 • C (Figure 4a). The suitable range of Precipitation of Warmest Quarter (bio18) was from 500 to 900mm with a peak at 600mm (Figure 4b). The response curve of Min Temperature of Coldest Month (bio6) showed that habitat suitability of tea is associated with areas where the optimal temperature ranged over 0 • C, but it has a fluctuant range at 10 to 16 • C (Figure 4c). The response curve of Mean Temperature of Warmest Quarter (bio10) indicated that 27 • C would be conducive for the development of tea in Warmest Quarter (Figure 4d). The response curve of Precipitation of Driest Month (bio14) demonstrated that there was a wide range of precipitation from 17 to 180 mm to adapt to the tea habitat, and it reaches the peak at 48 mm (Figure 4e). The highest suitability of Mean Diurnal Range (bio2) was 7 • C, with the decreasing trend beyond 7 • C (Figure 4f). The response curve of Annual Precipitation (bio12) showed that tea prefers precipitation between 1300 and 1800 mm (Figure 4g). Precipitation Seasonality (bio15) is calculated from the ratio of the standard deviation between monthly total precipitation and monthly average total precipitation, indicating that the seasonal value of the highest precipitation of tea was about 70, and habitat suitability decreased with the increase of variation (Figure 4h). Based on the response curves of Isothermality (bio3), the highest suitability for tea occurred in the areas where bio3 was from 25 to 36, and over 45 (Figure 4k). However, the effects of bio15 and bio3 on the model performance are less important than other variables because of the low value of percent contribution and permutation importance (Table 2). the ratio of the standard deviation between monthly total precipitation and monthly av-erage total precipitation, indicating that the seasonal value of the highest precipitation of tea was about 70, and habitat suitability decreased with the increase of variation ( Figure  4h). Based on the response curves of Isothermality (bio3), the highest suitability for tea occurred in the areas where bio3 was from 25 to 36, and over 45 (Figure 4k). However, the effects of bio15 and bio3 on the model performance are less important than other variables because of the low value of percent contribution and permutation importance ( Table 2).

Simulated the Distribution of Suitable Habitats in Current Climate Condition
Based on the current bioclimate variables and tea distribution records, a climatically suitable habitat map was created by ArcGIS 10 ( Figure 5). The result shows that the area from latitude 18° N to 39° N and from 93° E to 122° E was the primary potential suitable region of tea in China. The suitable habitats (including Marginal, Medium, and Optimal suitable habitat) accounted for approximately 26.89% (259.15 × 10 4 km 2 ) of the total land area, primarily located in 20 provinces or autonomous regions: Hainan, Yunnan, Guizhou, Guangdong, Guangxi, Fujian, Zhejiang, Hubei, Hunan, Jiangxi, Anhui, Chongqing, Jiangsu, Taiwan, central and southeast of Sichuan, southeast of Gansu, central and eastern

Simulated the Distribution of Suitable Habitats in Current Climate Condition
Based on the current bioclimate variables and tea distribution records, a climatically suitable habitat map was created by ArcGIS 10 (  (Table S1) of the total land area, including the provinces of Chongqing, Guizhou, Guangxi, north of Guangdong, south of Zhejiang, northwest of Fujian, southwest of Yunnan, southeast of Sichuan, and southwest of Hubei, which was mainly located in the mountainous areas of east and southwest Yun-Gui Plateau, and hills of the southeast coast of China. The Medium suitable habitat accounted for 11.22% (108.15 × 10 4 km 2 ) of the total land area, including Hunan, southern Guangdong, central and eastern Yunnan, southern Anhui, eastern Hubei and Sichuan. Meanwhile, the area of Marginal suitable habitat was lower than that of Medium and Optimal, consisting of 7.42% (71.51 × 10 4 km 2 ) of the total land area, mainly located at southeast of Gansu, central and eastern Shandong, southern Shaanxi, central and southern Henan and some small regions of south Tibet and south Sichuan.
including the provinces of Chongqing, Guizhou, Guangxi, north of Guangdong, south of Zhejiang, northwest of Fujian, southwest of Yunnan, southeast of Sichuan, and southwest of Hubei, which was mainly located in the mountainous areas of east and southwest Yun-Gui Plateau, and hills of the southeast coast of China. The Medium suitable habitat accounted for 11.22% (108.15 × 10 4 km 2 ) of the total land area, including Hunan, southern Guangdong, central and eastern Yunnan, southern Anhui, eastern Hubei and Sichuan. Meanwhile, the area of Marginal suitable habitat was lower than that of Medium and Optimal, consisting of 7.42% (71.51 × 10 4 km 2 ) of the total land area, mainly located at southeast of Gansu, central and eastern Shandong, southern Shaanxi, central and southern Henan and some small regions of south Tibet and south Sichuan.

Changes in the Spatial Distribution of Habitats Suitability in the Future
Using the same classification standard in Figure 5, the future habitats suitability for tea can be divided into four classes ( Figure S1). The distribution of suitable habitats category code is mapped in Figure 6 by comparing the current and future climate conditions under SSP126, 245, 370, and 585 for the years 2050s (2041-2060), 2070s (2061-2080), and 2090s (2081-2100). The category code "0" indicates Unsuitable habitat, "1" indicates Marginal suitable habitat, "2" indicates Medium suitable habitat, and "3" indicates Optimal suitable habitat. The regions with increased suitability in the future will mainly be located in the central and eastern Hubei Province, southern Anhui and Jiangsu Provinces, southern Yunnan Province, and northern Shandong and Henan Provinces. Under the current climatic conditions, these regions are in the Marginal or even Unsuitable habitat. However, the areas with reduced suitability are located at the Optimal area of current climate,

Changes in the Spatial Distribution of Habitats Suitability in the Future
Using the same classification standard in Figure 5, the future habitats suitability for tea can be divided into four classes ( Figure S1). The distribution of suitable habitats category code is mapped in Figure 6 by comparing the current and future climate conditions under SSP126, 245, 370, and 585 for the years 2050s (2041-2060), 2070s (2061-2080), and 2090s (2081-2100). The category code "0" indicates Unsuitable habitat, "1" indicates Marginal suitable habitat, "2" indicates Medium suitable habitat, and "3" indicates Optimal suitable habitat. The regions with increased suitability in the future will mainly be located in the central and eastern Hubei Province, southern Anhui and Jiangsu Provinces, southern Yunnan Province, and northern Shandong and Henan Provinces. Under the current climatic conditions, these regions are in the Marginal or even Unsuitable habitat. However, the areas with reduced suitability are located at the Optimal area of current climate, including the junction of Zhejiang and Fujian Provinces, and the middle of Guangxi and Guangdong Provinces. Two major observations can be drawn as shown in Figures 6 and 7. Firstly, the results showed that the suitability category code can change to sub-level, but hardly change across levels ( Figure 7); that is, suitability can be increased from category code "1 to 2", but hardly from "1 to 3", which means that the impact of climate change on suitability is gradual and periodic. Secondly, the increased area of Marginal and Medium could lead to the increasing of suitable area, especially if the area of "0 to 1" leads to the expansion of the new habitat, which has different trends in different scenarios (Figures 6 and 7).
In SSP126, there is a downward trend on suitability because of the suitability area of increasing category code (57.54 × 10 4 km 2 ) less than decreasing code (75.26 × 10 4 km 2 ) (Table S2), especially the area of category code 3 to 2 in the 2050s and 2070s, corresponding with the suitability decreasing in the coastal areas in southeast China, while the suitability of northern China increases, especially in southwestern Shandong. The area of Marginal and Medium habitat would increase in all study periods, but the area of Optimal habitat is lower than current in in 2070s and 2090s ( Figure S1 and Table S1). Suitable planting area may be stable in SSP245 because the area of increasing of suitability is 77.70 × 10 4 km 2 while with decreasing of suitability is 70.58 × 10 4 km 2 . It is worth noting that the newly occurring tea suitable area (0 to 1) in the 2090s are mainly located over the boundary of Shandong and Henan provinces and the northeast to the central of Hebei province in this scenario, but the Medium area increases to 124.87 × 10 4 km 2 in the 2090s ( Figure S1 and Table S1), which is second only to that of SSP370 of 128.92 × 10 4 km 2 in the 2090s, while the Optimal area increases first and then decreases. In SSP370, the reduced area (84.21 × 10 4 km 2 ) will be more than the increased area (75.95 × 10 4 km 2 ), which is mainly due to the reduction of a large area of category code "3-2", most of which is in the south of China, and even expanded to the southwest edge of Shaanxi in the 2090s ( Figure 6). Meanwhile, Optimal suitability will be increased in the south of Yunnan province in SSP370 with the area of 67.54 × 10 4 km 2 (12.60% of total land area) by comparison with the current climate ( Figure S1 and Table S1). The area in SSP585 will change significantly, and the suitable area of Shandong, Henan, Hubei, and Jiangsu will increase significantly with the increase of "0-1" and "1-2", resulting in a larger increase (79.44 × 10 4 km 2 ) in area as compared to the decrease (44.88 × 10 4 km 2 ) in area ( Figure 6 and Table S2). The southernmost part of the suitable area will reach the southern part of Hunan Province, with the increase of Marginal and decrease of Optimal ( Figure 6 and Table S2).
In summary, compared with the current conditions, climate change in the future may have a negative impact on the optimal suitable habitat, but a positive impact on Medium and Marginal suitable habitat. The optimal suitable habitat will be lost in southeast China, especially in Fujian and Zhejiang province. However, areas of Medium and Marginal suitable habitat will increase, replacing the Optimal suitable area, especially in Guangxi and Yunnan province. At the same time, the new suitable area will be expanded from the middle of Henan and Shandong province to the north. including the junction of Zhejiang and Fujian Provinces, and the middle of Guangxi and Guangdong Provinces. Two major observations can be drawn as shown in Figure 6 and Figure 7. Firstly, the results showed that the suitability category code can change to sublevel, but hardly change across levels ( Figure 7); that is, suitability can be increased from category code "1 to 2", but hardly from "1 to 3", which means that the impact of climate change on suitability is gradual and periodic. Secondly, the increased area of Marginal and Medium could lead to the increasing of suitable area, especially if the area of "0 to 1" leads to the expansion of the new habitat, which has different trends in different scenarios (Figures 6 and 7).  In SSP126, there is a downward trend on suitability because of the suitability area of increasing category code (57.54 × 10 4 km 2 ) less than decreasing code (75.26 × 10 4 km 2 ) (Table S2), especially the area of category code 3 to 2 in the 2050s and 2070s, corresponding with the suitability decreasing in the coastal areas in southeast China, while the suitability of northern China increases, especially in southwestern Shandong. The area of Marginal and Medium habitat would increase in all study periods, but the area of Optimal habitat is lower than current in in 2070s and 2090s ( Figure S1 and Table S1). Suitable planting area may be stable in SSP245 because the area of increasing of suitability is 77.70 × 10 4 km 2 while with decreasing of suitability is 70.58 × 10 4 km 2 . It is worth noting that the newly occurring tea suitable area (0 to 1) in the 2090s are mainly located over the boundary of Shandong and Henan provinces and the northeast to the central of Hebei province in this scenario, but the Medium area increases to 124.87 × 10 4 km 2 in the 2090s ( Figure S1 and

Discussion
Climate change has a crucial factor of changing geographical distribution of many species [54][55][56]. It is shown from the observation facts that climate change has led to the extinction of some species, and many species have escaped from their habitats to avoid the terrible climatic environment [57]. Therefore, predicting the response from species to climate change, especially their potential habitats distribution, can be regarded as a great effort to help the planning and decision-making to improve the sustainability and adaptability of ecosystems [58]. Previous studies have shown that various plants and animal species can transfer suitable areas to high altitude and high latitude regions, so as to avoid the adverse effects of extreme weather and utilize available climate resources [59][60][61][62]. In order to formulate countermeasures to adapt to climate change, it is necessary to predict the impact of climate change on the potential distribution pattern of species [63]. The MaxEnt model has been widely used to predict the potential distribution changes of species under different climatic scenarios because of its simple operation, high fitting accuracy and low requirement on sample numbers [64]. As a famous edible economic species, the potential distribution of tea and the impacts of climate change on it are little known. This study comprehensively analyzed the changes of suitable habitat of tea plants in China under the latest climate scenario, which provided an important reference for formulating sustainable development policies, planning and management strategies.

The Performance of MaxEnt Model
In this study, 410 species location data and 9 environmental variables were analyzed by MaxEnt to predict the distribution of tea in mainland China. Area under the curve (AUC) is an important model quality index to evaluate the model accuracy, because it provides a single measure of the accuracy of the overall model independent of selection threshold [65,66]. Generally, models with an AUC value greater than 0.75 might show good performance for niche model [23,67]. The AUC value of this study under current climate conditions is greater than 0.9, which indicates that the model has excellent fitting ability, and the results can objectively reflect the distribution of tea species. Furthermore, our simulation of the potential distribution of tea is only based on the occurrence data in native regions rather than occurrence data in exotic areas, which indicated that the climatic niche of tea is the realized niche instead of the fundamental niche [68,69]. Previous studies have shown that the realized niche is always smaller than the fundamental niche [70,71]. Our study of tea potential habitat was only based on the impact of climate variables and did not include the hydro-geological characteristics, due to data limitations.

Variable Influence on Tea Suitable Habitats Distribution
MaxEnt is a powerful suitability analysis model based on geographic information, which is mainly used for species prediction based on "existing" data, and environmental variables are very important for the model [23]. Many researchers use this model to predict the influence of all environmental factors or major environmental factors on the geographical distribution of species [53,72]. However, environmental variables are selected according to the different requirements, which could bring about problems such as autocorrelation and multicollinearity among these variables, thus negatively affecting the simulation results [73,74]. This study predicted the current and future climate suitability habitats based on the nine key climate factors with close correlation to tea distribution. The Temperature Annual Range (bio7) was found to be the most influential variable, with the permutation importance of 51% according to the results of the Jackknife test. Previous studies showed that temperature plays an important role in tea production and quality [14,75,76]. Meanwhile, temperature can also affect the distribution of tea by affecting water absorption, photosynthesis, transpiration, respiration, reproduction, and growth. The key Temperature Annual Range was 26 • C, which is higher than the result of Jayasinghe in Sri Lanka [24]. As a kind of rain-fed crop, tea growing needs a large amount of rainfall to maintain its life cycle and ensure the production of fresh leaves. The growth and production of tea plants are threatened by the uncertainty of rainfall brought by climate change [12]. Furthermore, climate change, especially global warming, changes the distribution pattern of precipitation. The precipitation in the dry season is a key indicator of tea survival [77]. In this study, it was found that the Precipitation of Driest Month (bio14) contributed a lot to the MaxEnt model (63.9%), ranging from 17 mm to 180 mm, and was used to evaluate the suitability of tea habitat. Precipitation will be diversified in the future [78], and extreme weather such as rainstorm or drought is expected to occur repeatedly, which will play a key role in the distribution of tea in China in the future.

Current and Future Dynamic Trend of Tea Suitable Habitats Distribution under Climate Change
Previous studies have shown that tea plants in China are distributed from 94 • E to 122 • E longitude, and from 18 • N to 37 • N latitude, but the economic planting areas are concentrated in the region between 94 • E to 121 • E eastwards and 18 • N to 32 • N northwards [78,79]. Our simulated potential distribution areas of tea are mainly located in the southeast of Tibet to the north of Shandong province, ranging from 93 • E to 122 • E and 18 • N to 39 • N, including 20 provinces or autonomous regions. In this study, only the natural geographical distribution of tea was considered, and areas with low temperature in winter, such as western Sichuan and parts of Tibet, will become potential low-level suitable areas.
Because of the uncertainty of future climate, there are obstacles to assessing the impacts of climate change on species distribution [80]. Numerous previous findings proved that global warming would have negative effects on ecosystems and species [81,82]. Our study indicated that the suitable area of tea species will expand in the future by the increasing of Marginal and Medium suitable area, mainly located in the mountainous regions such as southern Yunnan and Anhui Provinces by comparing spatial distribution and calculating area change. This indicated that under the current climate background, the area of the Optimal suitable area has reached the highest value, and there is little potential for future area increase. Other low-level suitable areas will become more conducive to production due to warming, and the potential for future increase is not great. In the period of 2081-2100, the Optimal habitat will decrease at SSP126, but there is an opposite trend between Medium and Marginal. However, in the high emission scenario SSP585, there are negative impacts of climate on Marginal in the 2070s and Medium in the 2050s, respectively. The areas that can be used for expanding production in the future are mainly located in Marginal and Medium areas, while the newly-increased production areas with great development potential are mainly located in the edge areas such as Shandong, Yunnan, and Henan. On the one hand, because of its high economic value, tea tree is different from other species, and it is mainly managed according to the management method of Gramineae crops in actual production. The development of tea to cope with climate change is highly dependent on social and economic input, and the benefits brought by tea production could promote the excavation of potential distribution areas, but this process is meaningful only when the climatic conditions of a certain area meet the requirements of tea production. On the other hand, according to the physiological characteristics of tea tree, as an evergreen plant, it is suitable for growing in a humid and warm environment, with limited resistance to cold and heat [83,84]. Under the climatic conditions of increasing CO 2 concentration and temperature in the future, the temperature in northern Shandong and southern Shanxi will reach the minimum standard for tea growth, which will lead to the northward shift of the suitable regional boundary for its distribution. At the same time, due to the higher elevation in southwest China, the temperature in high altitude area will also rise, which will lead to the distribution boundary moving to southwest China. Moreover, the suitable area of tea will expand northward and west in response to global warming and social and economic development, suggesting that additional new suitable habitats will be created for Chinese tea production with the future climate change.

Other Uncertainty of Tea Distribution in the Future Climate Change
Due to changes of environmental factors caused by climate change, the suitable habitats of tea plants in China have changed, and it is difficult to formulate a suitable long-term tea production management plan for adapting to climate change [38]. The MaxEnt model is an important tool to simulate the suitable habitats distribution of the species, because it can predict suitable distribution areas under future climate change scenarios. It should be noted that, although the simulation results of the MaxEnt model have high accuracy, it does not mean that the suitable area is completely consistent with the actual distribution area [21]. The simulation result of the model is the maximum possible geographical distribution of tea, which is the maximum possible distribution range under ideal conditions. However, due to the limitation of data, the economic value of tea quality has not been considered, and not all simulated suitable areas have better quality and higher economic benefits. Moreover, the distribution of tea is affected not only by climate factors, but also by non-biological factors such as topography and soil. On the other hand, it is also influenced by social and economic development, human intervention, policies, and other human activities. Although there are many assumptions and uncertainties in the tea distribution model, this model is still the key data source for future suitability prediction [85], so as to evaluate scientific adaptation strategies at the level of species, community and ecosystem, and to offset the impact of future warming on tea distribution.

Conclusions
In this study, the present-day occurrence data of tea (Camellia sinensis L.) were used to predict the current and future potential suitable distribution habitats in China by MaxEnt model under four SSPs climate change scenarios. The main conclusions are as follows: (1) AUC value was over 0.9, which reveals that the model performance was excellent to predict the future suitable habitats of tea. The nine dominant climate factors were selected by correlation analysis and Jackknife test, in which bio14 (Precipitation of Driest Month) contributed the most, accounting for 63.9%, and bio7 was the most important to predict tea distribution, with the permutation importance accounting for 51%. (2) Under the current climatic condition, tea distribution areas were mainly located from south of Hainan to the north of Shandong province, including 20 provinces or autonomous regions. The most suitable climate zone was mainly located in the higher elevation zones including the central and eastern Guizhou, Chongqing, Zhejiang, and the coastal areas of Fujian Province. (3) Compared with current distribution, this study revealed that the suitable area would increase in the future climate change scenarios because of the expansion of Medium and Marginal suitable habitats in central and eastern Hubei Province, central Anhui and Jiangsu Provinces, southern Yunnan Province, and even northern Shandong and Henan Provinces. Only considering climate factors, tea production in western provinces such as Chongqing, Guizhou, Yunnan, and other provinces such as Shandong and Henan have great development space, so formulating the adaptative technology of tea production should be considered in these provinces.
The distribution map of tea potential habitat is helpful to strengthen land use management, create new tea gardens, increase farmers' welfare, avoid competing for arable land with food production, and make the future distribution prediction an important backbone for formulating the development plan of the tea industry. Climate factors are considered in the distribution model, but the other factors such as socio-economic development, human intervention, policies, and other human activities have always played a significant role in species distribution, which should be strengthened in future study.