Climate Change May Imperil Tea Production in the Four Major Tea Producers According to Climate Prediction Models

The threat of accelerating climate change on species distribution now and in the future is a topic of increasing research interest. However, little work has been undertaken to assess how shifting climates will affect the suitability of tea cultivation. Therefore, we used MaxEnt modelling to project the impact of current and future climatic scenarios on the potential distribution of tea across the four tea-producing countries of China, India, Kenya and Sri Lanka. Projections were made for the years 2050 and 2070 with three Representative Concentration Pathways (RCPs) using seven bioclimatic predictors under three global circulation models (GCMs). The current and future habitat suitability for tea predicted by the models produced a high accuracy rate, with high areas under the receiver operating characteristic curve (AUCs) for all tested RCPs under the three GCMs for the four countries. The mean true skill statistic (TSS) values for tea in Sri Lanka, Kenya, India and China were 0.80, 0.91, 0.91, and 0.74, respectively. The kappa values (k) of the current and future models for all four countries ranged from 0.40 to 0.75, which indicates that the overall performance of the model was good. The precipitation seasonality and annual precipitation were found to be the most influential variables in Sri Lanka and India, respectively, while annual mean temperature was the most effective contributor for determining the suitability of habitat for tea in Kenya and China. An important proviso is that some existing tea-growing areas will face reduced suitability for future tea cultivation suggesting that by 2050 there will be a drastic reduction in the optimal suitability by averages of 26.2%, 14%, and 4.7% in Kenya, Sri Lanka and China, respectively. The optimal suitability will be reduced by 15.1%, 28.6% and 2.6% in Kenya, Sri Lanka and China, respectively, by 2070. India displays an advantage in projected future climates as it gains optimal suitability areas of 15% by 2050 and 25% by 2070.


Introduction
Tea, the most widely consumed non-alcoholic beverage in the world, is produced from the tender young shoots of Camellia sinensis (L.) O. Kuntze. Global tea production continues to surge dramatically, reaching 5.8 million metric tonnes (t) in 2018 while the export quantity came to about 1.76 million t per year [1]. Major types of tea are all gaining increasing popularity across the world, thanks to their health benefits, great taste, and ample availability in the global market [2]. Driven by the increasing demand for tea, world tea consumption has risen by 4.5% to 5.5 million t over the decade to 2016 [1]. The tea industry is important in respect of its contribution to the gross national product (GNP), So far no one has studied how climate suitability will shift among major tea-producing countries in order to visualise a comprehensive picture of challenges and opportunities in the future tea market. While some modelling of tea projection has been undertaken, there has been no comparison between major tea producing countries using a similar set of global circulation models (GCMs), RCPs as well as time periods. Therefore, attention needs to be paid to conducting research to model current and future habitat suitability for tea using similar GCMs, RCPs and time periods for comparing how the major tea producing countries will fare in the future.
The primary objective of this study was to investigate the potential current and future tea distribution and climatic suitability on the basis of three GCMs for the years 2050 and 2070 under RCP scenarios of 2.6, 6.0, and 8.5. We highlight relative declines and changes in climate suitability for the world's major tea-producing countries of China, India, Kenya and Sri Lanka. We identify the main bioclimatic variables influencing tea cultivation and study the comparative shifts in areas suitable to growing tea. Quantifying climate characteristics of each tea producing country will enable more effective decision-making and policy formulating in order to sustain and enhance future tea production in the competitive international market. Further, it is worth looking to see if there is a migration of tea growing to newer and more suitable areas in the four selected countries if the traditional regions become less suitable.

Study Area and Its Ecological Significance
This study focused on the top four tea-producing countries, namely Sri Lanka, Kenya, India and China. The comprehensive characteristics of each country selected in this study are given in the Supplementary Materials Table S1.

MaxEnt Description
MaxEnt stands for maximum entropy modelling which has been widely used in niche modelling since it only needs presence-only species records [26]. The MaxEnt model was considered to be the best in both predictive efficiency and model stability compared to other species distribution models (SDMs) [27]. Also, MaxEnt predicts the habitat suitability of the species by analysing the occurrence data of the target species and environmental variables [28,29].

MaxEnt Modelling
The MaxEnt software (MaxEnt version 3.3.3, http://www.cs.princeton.edu/~schapire/MaxEnt/, accessed on 20 February 2020) algorithm was run with the default convergence threshold (10 −5 ) and the maximum number of iterations (5000). The cloglog output format was used and models were each run 10 times to measure internal model variability with cross-validation run type [26]. Maximum training sensitivity plus a specificity threshold was used with 50,000 background points to run the models which allowed the models to have adequate time for convergence. The 'fade-by-clamping' feature was used to avoid dubious projections by the model which removes heavily clustered pixels from the final predictions [30]. In addition, we evaluated the effect of different combinations of feature types and varying degrees of regularisation in the MaxEnt program. The auto-feature configuration performed best with a regularisation multiplier of two to calibrate our final model [27].
We produced bias files for habitat suitability modelling for each country to employ the target-group background method recommended by Phillips, Dudík [31] to avoid sampling bias. The information on the dependence of the predicted suitability on a specific variable as well as the range under which the variable reached its optimal suitability when provided by response curves [27].

Estimating Species Occurrence and Climatic Predictors
The presence-only species collection techniques were used in this study. The occurrence data of tea were obtained from three sources: published literature [23,[32][33][34][35], the Global Biodiversity Information Facility (GBIF, http://www.gbif.org, accessed on 24 February 2020) and Google Earth (http://ditu.google.cn/, accessed on 26 February 2020). A total of 4139 geo-referenced occurrence data of tea were collected within its known native distribution over four countries. Occurrence data of tea for each country are given in the Supplementary table (Table S2). However, for tea, this published literature and GBIF sources yielded very few records of the known native distribution of tea, especially from China and Sri Lanka. Hence, we included data from Google Earth Pro (Version 7.3.1, Google Inc.). We used expert knowledge to confirm the reliability of extracted occurrence points of each county since precise location information would ensure improved geographic accuracy and generate a better model. R software version 2.51 was used with an additional package of 'spThin' to eliminate spatial clusters of concurrence points. Occurrence points were thinned at 5 km with 100 iterations for each study area. Finally, a total of 536 occurrence data of tea were included in this study (Table S2) and the occurrence points were converted into a comma-separated values (CSV) format.

Bioclimatic Variables
For modeling purposes, the set of bioclimatic variables available on the WorldClim website has been widely used [23,36,37]. Generally, the 'current' climate data sourced from the WorldClim use the years 1950-2000 to calculate climatic averages. This is a common method to work with current climate condition (interpolations historical observed data, representative of 1950-2000 climates) as it is easy and straightforward. That does not, however, account for recent (post-2000) climate change. Current climate: for current climatic conditions , 19 'bioclimatic' variables were downloaded from the WorldClim database (http://www.worldclim.org/, accessed on 12 February 2020) at a spatial resolution of 30 arc-seconds (~1 km 2 ). Based on the previous study conducted by Jayasinghe and Kumar [23], the growth conditions of tea as well as the operability and feasibility of the data, seven bioclimatic variables [annual mean temperature (BIO1), mean diurnal range (BIO2), temperature seasonality (BIO4), annual precipitation, (BIO12) precipitation seasonality (BIO15), precipitation of driest (BIO14) and wettest month (BIO16)] were chosen for this study using the Pearson correlation coefficient test (|r| < 0.70). A strong correlation between the bioclimatic variables in the model may lead to misinterpretations due to the high level of collinearity among variables. Highly correlated variables (r > 0.80) can be excluded from the model based on the Pearson correlation coefficients and the less correlated variables can be taken (r < 0.8) for the MaxEnt analysis [38].
Future Climate: Three GCMs were used with 1 km 2 spatial resolution for the years of 2050 and 2070 to predict the future habitat suitability. The GCMS used were CCSM4 (Community Climate System Model, version 4), HadGEM2-ES (Hadley Centre Global Environmental Model 2-Earth System) and MIROCH-H (Model for Interdisciplinary Research on Climate) [39]. It has been reported that MICROC5 has a better simulation of the climatic parameters than its past versions [40] and it has been used for many other studies in the south Asian and Himalayan regions due to its improved capabilities in capturing features well [23,41,42]. As Kenya is one of the study areas in the present study, HadGEM2-ES was used, as previous studies indicated that the model simulations were close to the multimodal ensemble mean over Africa [43]. The CCSM4 has improved El Niño-Southern Oscillation (ENSO) variability with more reasonable frequency distribution [44] and it is anticipated that the rainfall and temperature are better modelled by CCSM4 in the Asian region [45]. The models were run using three RCPs of RCP2.6, RCP6.0 and RCP8.5 in order to investigate how suitable habitats of tea will change under several GCMs, viz., MIROC5, CCSM4 and HadGEM2-ES by 2050 and 2070, respectively, compared to the current climate distribution in all four tea-producing countries. The RCP2.6 and 6.0 have lower (2.6-3.0 W/m 2 ) and medium (6.0 W/m 2 ) radiative forcing levels, respectively, whereas RCP8.5 is considered as the high emission scenario (business as usual) that is based on a radiative forcing level of +8.5 W/m 2 .

Evaluation of Model Accuracy
Models were calibrated using 75% of the occurrence points as training data, and the remaining 25% were used for model validation. The area under the receiver operating curve (AUC), Cohen's kappa, and true skill statistic (TSS) were used to evaluate the model performance [27,46]. AUC is a measure of the overall performance of a model [47]. The AUC ranges from 0.5 (random) to 1.0 (perfect discrimination), and values >0.8 indicate satisfactory performance [48]. TSS and Cohen's kappa are the threshold-dependent measures of model performances that evaluate their classification accuracy after selecting a threshold value [47]. We used a suitability threshold of 0.5, which has been widely used in species distribution models [49]. The TSS and Cohen's kappa values >0.8 were suggested as excellent, 0.4-0.8 useful, and <0.4 poor model performance [50]. The confusion matrix, formulas and measures of accuracy indicators are given in the Supplementary tables (Table S3).
Furthermore, we used three methods to estimate the contributions of climate indices to simulated models [51]: (1) percentage contribution, (2) permutation importance and (3) the jack-knife test. We explored the contribution of environmental variables by: (a) assessing their permutation importance, and (b) with jack-knife tests, which indicate the change in model fit when subsequently withholding each predictor and refitting models [49].

Reclassification and Change Detection
The suitability maps, current and future, produced by the model ranged from 0 (least) to 1 (most). To define habitat versus non-habitat for tea, we used the "maximum training sensitivity plus specificity" threshold that provided very accurate predictions. For further analysis, we classified tea habitat suitability into four classes using the "natural/jenks" option in ArcGIS, i.e., 'optimal' (>0.6), 'medium (0.4-0.6), 'marginal' (0.2-0.4), and 'unsuitable' (<0.2). The area calculation for the identified classes was carried out in the ArcGIS setting.
Then, the change detection method was used to check for spatial changes in suitability areas from current to future scenarios. The suitability values for the future climate were subtracted from the current climate to see where the changes were high and low. This gave an indication of how the climate will have a spatial impact on the tea suitability of the four countries. Tea suitability is a term used in the present study to highlight those areas that have suitable habitats, in terms of climatic conditions, for successful tea cultivation in the current and/or future scenarios.

Habitat Suitability Model Validation, Parameter Sensitivity
Potentially suitable habitats with optimal suitability thresholds for tea were predicted in the four countries where populations of tea are already known to exist. The average test AUCs for the replicate runs were 0.92 (standard deviation (SD) = 0.009), 0.98 (SD = 0.003), 0.97 (SD = 0.005), and 0.90 (SD = 0.013) for Sri Lanka, Kenya, India and China, respectively for current climate. AUCs for species distribution maps for all four countries showed a high accuracy, indicating the reliability of model predictions ( Table 1). The mean TSS values for tea in Sri Lanka, Kenya, India and China were 0.80, 0.91, 0.91, and 0.74, respectively, which indicates very good model performances. The kappa values (k) of the current and future models for all four countries were above 0.40 (0.4 < k < 0.75), showing that the overall model performance was good (Table 1).  Table 2 gives permutation importance of the predictor variables to the MaxEnt model of each country. In Sri Lanka, current and future tea habitats are highly influenced by precipitation seasonality (BIO15) and annual mean temperature (BIO1), accounting for almost 80% of the model prediction. For Kenya, annual mean temperature (BIO1) is the significant climatic predictor for current (M = 47%) and future (50%) tea cultivation. For India, annual precipitation (BIO12) and annual mean temperature (BIO1) are the most influential variables, contributing more than 50% to the tea distribution models under current and the future climate. For China, annual mean temperature (BIO1) (M = 46%) and precipitation of the wettest quarter (BIO16) (M = 31%) are the two most important bioclimatic variables influencing the distribution of tea under current as well as future climatic conditions.

Key Environmental Factors Influencing Habitat Distribution
The jack-knife analysis of variable importance also resulted in annual mean temperature (BIO1) and precipitation seasonality (BIO15) being the two most important predictors of tea distribution in Sri Lanka while annual precipitation (BIO12) was the most significant contributor to the habitat suitability distribution of tea in Kenya and India. It also showed that annual mean temperature (BIO1) and mean diurnal range (BIO2) were the main variables for the impact on the distribution of suitable habitat for tea in China. The results of jack-knife analysis is given in the supplementary material ( Figure S1).

Importance of Bioclimatic Variables in the Habitat Suitability of Camellia Sinensis
Response curves reflect the quantitative relationship between the MaxEnt predicted probability of tea habitat suitability and bioclimatic variables (Figure 1). The x-axis shows how each bioclimatic variable changes as the cloglog prediction changes between 0 and 1. All bioclimatic variables were set to their average value over the set of presence occurrences. According to the response curves of annual mean temperature (BIO1) (Figure 1a), high probabilities of habitat suitability for tea were skewed sharply towards a specific range (10-25 • C) in both India and China. In Sri Lanka, high probabilities of presence were predicted when BIO1 ranged from 13 to 28 • C, while temperatures around 20 • C are ideal (Figure 1b). The highest habitat suitability for tea shows a symmetric distribution with the BIO1 ranging from 15-20 • C in Kenya. Figure 1b shows that habitat suitability of tea in Sri Lanka, China, and India is associated with areas where the mean diurnal temperature (BIO2) values ranged from 4.5 to 10 • C while the diurnal temperature of Kenya exhibited a range from 13 to 16 • C. The response curve for BIO4 (Figure 1c) showed that the habitat suitability of tea in Sri Lanka and India drastically decreases with a high-temperature variation that is above 40 • C and 48 • C, respectively. The response curve of annual precipitation (BIO12) (Figure 1d) showed that tea preferred high precipitation ranges from 2000 to 5000 mm in Sri Lanka. Its habitat suitability was strongly constrained by BIO12 at 2500 mm in India and 1700 mm in China. The corresponding response curves for precipitation of the driest month (BIO14) show a clear downward trend (Figure 1e), representing positive relationships between increasing values of BIO14 and predicted habitat suitability for tea in all tested tea producing countries. The habitat suitability of tea of four countries changed by varying degrees with increasing precipitation of the wettest quarter (BIO16) (Figure 1g). and precipitation seasonality (BIO15) being the two most important predictors of tea distribution in Sri Lanka while annual precipitation (BIO12) was the most significant contributor to the habitat suitability distribution of tea in Kenya and India. It also showed that annual mean temperature (BIO1) and mean diurnal range (BIO2) were the main variables for the impact on the distribution of suitable habitat for tea in China. The results of jack-knife analysis is given in the supplementary material ( Figure S1).

Importance of Bioclimatic Variables in the Habitat Suitability of Camellia Sinensis
Response curves reflect the quantitative relationship between the MaxEnt predicted probability of tea habitat suitability and bioclimatic variables (Figure 1). The x-axis shows how each bioclimatic variable changes as the cloglog prediction changes between 0 and 1. All bioclimatic variables were set to their average value over the set of presence occurrences. According to the response curves of annual mean temperature (BIO1) (Figure 1a), high probabilities of habitat suitability for tea were skewed sharply towards a specific range (10-25 °C) in both India and China. In Sri Lanka, high probabilities of presence were predicted when BIO1 ranged from 13 to 28 °C, while temperatures around 20 °C are ideal (Figure 1b). The highest habitat suitability for tea shows a symmetric distribution with the BIO1 ranging from 15-20 °C in Kenya. Figure 1b shows that habitat suitability of tea in Sri Lanka, China, and India is associated with areas where the mean diurnal temperature (BIO2) values ranged from 4.5 to 10 °C while the diurnal temperature of Kenya exhibited a range from 13 to 16 °C. The response curve for BIO4 (Figure 1c) showed that the habitat suitability of tea in Sri Lanka and India drastically decreases with a high-temperature variation that is above 40 °C and 48 °C, respectively. The response curve of annual precipitation (BIO12) (Figure 1d) showed that tea preferred high precipitation ranges from 2000 to 5000 mm in Sri Lanka. Its habitat suitability was strongly constrained by BIO12 at 2500 mm in India and 1700 mm in China. The corresponding response curves for precipitation of the driest month (BIO14) show a clear downward trend (Figure  1e), representing positive relationships between increasing values of BIO14 and predicted habitat suitability for tea in all tested tea producing countries. The habitat suitability of tea of four countries changed by varying degrees with increasing precipitation of the wettest quarter (BIO16) (Figure 1g).

Current Habitat Suitability Distribution
The potential distribution maps show that the suitable habitats of tea will change in the future, compared to the current climate (Figure 2a, Figure 3a, Figure 4a, Figure 5a). In Sri Lanka, the predicted medium and optimal suitability areas for tea under the current climate conditions are 5769 km 2 (8.8%) and 6090 km 2 (9.3%), respectively ( Figure 2a) and (Table 3a). The medium and optimal current habitat for tea in Kenya represent 11,078 km 2 (4.5%) and 10,333 km 2 (4.2%), respectively ( Figure 3a) and (Table 3b). In India, an area of 66,725 km 2 (2%) has optimal climate suitability while 131,268 km 2 (4%) has medium suitability for tea distribution under the current climate ( Figure 4a) and (Table 3c). In China, the simulation output shows that 1,804,999 km 2 (19.3%) and 1,003,219 km 2 (10.7%) of land have optimal and medium suitability for tea under the current climate, respectively ( Figure 5a) and (Table 3d). The potential distribution maps show that the suitable habitats of tea will change in the future, compared to the current climate (Figures 2a-5a). In Sri Lanka, the predicted medium and optimal suitability areas for tea under the current climate conditions are 5769 km 2 (8.8%) and 6090 km 2 (9.3%), respectively ( Figure 2a) and (Table 3a). The medium and optimal current habitat for tea in Kenya represent 11,078 km 2 (4.5%) and 10,333 km 2 (4.2%), respectively ( Figure 3a) and (Table 3b). In India, an area of 66,725 km 2 (2%) has optimal climate suitability while 131,268 km 2 (4%) has medium suitability for tea distribution under the current climate ( Figure 4a) and (Table 3c). In China, the simulation output shows that 1,804,999 km 2 (19.3%) and 1,003,219 km 2 (10.7%) of land have optimal and medium suitability for tea under the current climate, respectively ( Figure 5a) and (Table 3d).

Future Habitat Suitability Distribution and Change Detection
The projected habitat suitability areas for tea for future climate scenarios for both the years 2050 and 2070 (Figures 2a-5a) (Table 3) vary compared with the current distribution in all four tea-producing countries. The details of suitable habitats for tea in Sri Lanka, Kenya, India and China under current and future climate conditions, predicted using the three GCMs separately, are given in supplementary tables (Tables S4-S7).
There will be a decrease in habitat suitability for tea under a future climate compared to the current climate in Sri Lanka (Table 3a). Compared to the mid-and high elevation areas, most of the

Future Habitat Suitability Distribution and Change Detection
The projected habitat suitability areas for tea for future climate scenarios for both the years 2050 and 2070 (Figure 2a, Figure 3a, Figure 4a, Figure 5a) (Table 3) vary compared with the current distribution in all four tea-producing countries. The details of suitable habitats for tea in Sri Lanka, Kenya, India and China under current and future climate conditions, predicted using the three GCMs separately, are given in supplementary tables (Tables S4-S7).
There will be a decrease in habitat suitability for tea under a future climate compared to the current climate in Sri Lanka (Table 3a). Compared to the mid-and high elevation areas, most of the optimal and medium suitability areas at low elevation (e.g., Galle and Matara) will be lost at a higher rate for all tested RCPs (Figure 2a). In the year 2050, most current optimally suitable areas for tea in the central hills will increase. The change detection maps for Sri Lanka (Figure 2b), clearly show that southern parts of the current tea-growing areas will decrease in climate suitability in the future. By 2050, under average GCMs for RCP2.6, the southern part of Ratnapura district will gain climate suitability.
Compared to current habitat suitability in Kenya, the average optimal and medium suitability will be reduced in the future. In particular, the optimal class around Mt Elgon and Mbeere will be totally absent under all RCPs by 2050 and 2070. Gucha and Kiambu areas will also see drastic reduction in optimal climate suitability for tea under RCP2.6 by 2050 and it will be completely absent under RCP6.0 and 8.5 by 2050, as well as under all RCPs by 2070. It is anticipated that there will be an increase in the suitability areas in Buret, Kericho and southern parts of Nandi in Kenya by 2050 and 2070, while suitability will decrease in some parts of Meru, Gucha and southern parts of Mt Elgon (Figure 3b).
In India, the optimal areas will expand in Sikkim, Manipur, upper Assam and Karnataka regions under projected future climates. However, overall medium suitability for tea will be significantly lost in all GCMs under RCP2.6 by 2070, especially in Meghalaya and Tripura regions. The suitable regions will increase in the future in Karnataka, some parts in Assam, Manipur, Uttarakhand and Himachal Pradesh while they decrease in some parts of Assam and northwest districts of Himachal Pradesh.
By the year 2070 northern and eastern of Assam will gain suitability under RCP6.0 and 8.5, respectively ( Figure 4b).
By 2050 and 2070, optimal and medium suitability areas in China show large variations compared to the current climate. It is clear that optimal and medium suitability areas for tea increase in Yunnan region in the future. The range of optimally suitable habitat would decrease in south Guangxi, Jiangsu and Shandong by 2070. The area of optimally suitable habitat that may be lost was 76,398 km 2 by 2050 and 55,481 km 2 by 2070. By 2050 under the average of three GCMs for RCP2.6, it is projected that the climate will become unsuitable for the currently optimal tea cultivated areas in parts of Hunan, Shandong, Yunnan and Sichuan ( Figure 5b) and (Table 3d). Also, current medium and marginal areas in Shandong and Hubei will see a reduction in suitability under RCP6.0 by 2050. By 2070, suitability will be drastically reduced under the average of tested GCMs for RCP8.5, especially in some areas of Jiangsu, Anhui, Hunan, Guangxi, southern parts of Chongqing and the middle part of Jiangxi and Hubei.

Percentage of Suitability Area Loss/Gain for Camellia Sinensis from the Current to the Future Climate Scenarios in the Four Major Tea-Producing Countries
In Sri Lanka, a high percentage of medium suitability land will be lost under RCP2.6, indicating an overall loss of 21% by the year 2050. Tea will experience loss in areas of optimal suitability at higher percentages under RCP8.5 by 2070 (Table 4a). The greater percentage of loss (22%) in optimal suitability for tea is recorded under HadGEM2-ES under both future climates (Table S8). The average loss of percentage of optimal climate suitability will be 14% and 15% by the years 2050 and 2070, respectively.
Proportionately, medium and optimal habitat suitability categories in Kenya will decline by a high percentages under future climate projections. The reduction of average loss in the medium suitability for tea will be 39% and 38%, respectively, by 2050 and 2070, while the loss in habitat area of optimally suitable land will be 26.2% and 27%, respectively, by 2050 and 2070 (Table 4b). A higher percentage of loss in optimally suitable areas will be seen under HadGEM2-ES and MIROC5 compared to CCSM4 by both years of 2050 and 2070 (Table S9).
In India, projected optimal climate suitability will increase by 12.4% and 21.9% by 2050 and 2070, respectively, under all RCPs (Table 4c). The highest gain in the optimal areas (40.5%) in 2070 will record in RCP6.0 under HadGEM2-ES (Table S10). However, medium suitability areas will decrease by 24.2% by 2050 and 21.6% by 2070 (Table 4c). Note: A summary of average area (km 2 ) of optimal suitability with the percentage of loss (−) or gain (+) is given in the last column compared to current climate.

of 21
For China, the results are explicit that future climate in both 2050 and 2070 can affect the optimal suitability of tea habitats, mainly under RCP6.0 in 2050 and RCP8.5 in 2070. The average percentage gain of medium climate suitability areas will be 5.7% and 7.2%, respectively, by the year 2050 (Table 4d). Climate change may precipitate significant loss in areas of optimal and medium suitability in China at higher percentages under HadGEM2-ES, representing an average percentage loss of 25% and 26% by 2050 and 2070, respectively. Moreover, an average of 5.6% of medium suitability area will be gained under all RCPs under MIROC5 (Table S11).
As per the summary of percentage of loss (−) or gain (+) of optimal suitability areas (%) in the four tea-producing countries in 2050 and 2070 (Table 4), it is clearly seen that Sri Lanka, Kenya and China will lose more optimal suitability areas in the future whereas India will gain an advantage as it gains optimal suitability areas, compared to the current climate.

Discussion
Climate change poses enormous challenges for the agriculture sector across the globe [52,53]. Tea is a relatively delicate plant and is especially sensitive to changes in climate [6]. This study shows that MaxEnt models fitted with presence-only data can provide strong estimates of habitat suitability for tea species in four different countries, with the high AUC, TSS and kappa values (Table 1). Our analysis further shows that the impacts of climate change on habitat suitability of tea would be very variable for the four major tea-producing countries. Overall, the impact of future climate in Sri Lanka, Kenya and China is predicted to be negative, while India would gain suitable habitats for tea in the future. The findings highlight that by 2050, a potential future decline of optimal suitability areas of 26.2%, 14%, and 4.7% in Kenya, Sri Lanka and China, respectively. By 2070, the optimal habitat suitability is also projected to decline by 28.6%, 15.1% and 2.6% in Kenya, Sri Lanka and China, respectively. India will benefit from future climate considering the optimal suitability area is likely to expand by 15% and 25% by the years 2050 and 2070, respectively (Tables 3 and 4). All the four countries would maintain some suitable habitats, although within each country this would vary a lot.
This study projects negative consequences of climate change on tea cultivation in Sri Lanka, as previously suggested by Jayasinghe and Kumar [23]. The potential distribution of tea in the current climate was in line with the existing distribution of tea cultivation in Sri Lanka (Figure 2). The highest percentage of areas having medium and optimally suitable climatic conditions will decrease under HadGEM2-ES and RCP2.6 by 2050 ( Figure 2). Since the radiative forcing of RCP2.6 peaks at 3 W/m 2 in 2050 and returns to 2.6 W/m 2 by 2100, the effect of RCP2.6 will be less in 2070 than in 2050 [54]. Under RCP2.6 using MIROC5 in 2050, a significant expansion occurs in the areas of optimal suitability distribution in Sri Lanka ( Figure 2). This may due to the fluctuations in temperature and precipitation seasonality in the future climate under different GCMs (Table S12); tea habitats may shift from marginal to unsuitable, or optimal to medium, or vice versa. The most noticeable effect of climate change is the decrease in those areas that display "optimal" suitability, with approximately 13.7% and 15% of the decline in the field from current to future scenarios. In business as usual scenario (RCP8.5), the tea-producing area further concentrates in the higher elevated areas. As temperature is highly related to elevation, tea habitats having optimal climatic suitability may shift towards the high elevation areas in the future (Figure 2), whereas at lower altitudes the areas will suffer from a warming climate [55] that may decrease their climate suitability for tea.
Our model predicts that tea is unlikely to be established in the areas of Kenya where it had not previously existed. The tea-growing areas in Kenya will decline by 26.2% and 26.7% by 2050 and 2070, respectively ( Figure 3) and (Table 3), which may have severe repercussions for the sustainability of tea production. These results are quite similar to the previous study undertaken by Han,Li [56] as they indicated that suitability of tea-growing areas is expected to decline by 22.5% by the year 2075. Our results indicated that the impacts were likely to be more severe on the west side of the Rift Valley in Kenya, a major tea-growing area in the country, as more optimal, medium and marginal suitability habitats convert to unsuitable habitats by 2050, which is comparable to the previous study [57].
These outcomes can be caused by the predicted variations in temperature and rainfall in Kenya under future climate (Table S12). A rise in mean air temperatures beyond the optimal temperature for shoot growth and development (i.e., 23.5 • C) is likely to occur. The uneven distribution of rainfall and prolonged drought could lead to undesirable conditions for tea plants [58]. The result is that the distribution of the suitability of tea growing in the existing tea-growing regions of Kenya will decrease mainly due to the above climate shifts. In Kenya, future optimal habitat suitability will be reserved in high elevations (2000 to 2300 m), especially the higher altitudes around Mount Kenya (Figure 3), which is similar to the previous study conducted by Leshamta [34].
For India, our results revealed that the optimal tea growing areas will expand in Sikkim, Manipur, upper Assam and Karnataka regions under projected future climate (Figure 4). In the current climate, the south bank region which covers parts of Upper Assam and Cachar are optimal suitable regions for tea whereas the north bank region is comparatively less suitable, which is similar to the results obtained from the Indian Tea Association [59]. Our study also predicted some areas in Assam, such as Sikkim, Meghalaya and Mizoram will more suitable by 2050 and 2070. By contrast, the climate change forecasting maps developed by the Ethical Tea Partnership organisation [21] and the Tea Research Association (TRA) demonstrated that much of Assam will become far less suitable for tea farming. One of the reasons for this dissimilarity could be the uneven distribution of occurrence data, or that we considered different GCMs and RCPs and different bioclimatic variables. A clear increase in temperature in the future would allow areas in high and medium elevations (i.e., Assam Manipur, Darjilling, Dooars, Teri and Arunachal Pradesh) to reach the optimal temperature range for tea cultivation, and would allow for expansion of tea-growing areas. Drought and low rainfall in India have been observed since the 1960s and it has had a huge impact on tea cultivation in past [60]. However, rainfall is projected to increase over all the states with very low change in the western part and highest change in the north-eastern and southern regions under RCP8.5 [61]. In the future, areas with unfavourable temperatures and rainfall conditions for tea growth could shift to desirable climate conditions compared to the current scenario. This is the reason why the future habitat distribution of tea in India is not negatively affected compared to other countries.
The modelled current and future optimal suitability habitats are concentrated around existing tea-growing areas in China (Fujian, Yunnan, Zheijiang, Sichuan, Hunan, Hubei and Anhui) ( Figure 5). By 2050 and 2070, optimal and medium suitability areas in China show large variations under different GCMs. The reason for this result may be due to the harsh weather conditions in China over past years [10,62]. The results of this study are consistent with previous studies showing that the habitat suitability of the tea plant species is predicted to improve in some regions under climate change conditions [63,64]. It is very clear that optimal and medium suitability areas for tea increase in Yunnan regions under all tested RCPs in the future. In China, the average temperature for the near-term forecast (2020-2050) is similar for both scenarios of RCP6.0 and RCP8.5. Nevertheless, for the long-term forecast (2060-2090), the average temperature for the RCP8.5 scenario is higher than for the RCP6.0 scenario. The temperature has increased in the near term by a median of 1.1 • C and 1.3 • C for RCP scenarios 6.0 and 8.5, respectively, and projected to increase in the long term by a median of 1.7 • C and 2.9 • C for the RCP scenarios 6.0 and 8.5, respectively [65]. This may be the reason for decreasing optimally suitable habitat in south Guangxi, Jiangsu and Shandong in 2070 under the scenario of RCP 8.5. The main tea-production area could gradually shift from south to north, depending on future climate projections in China [33].
According to the relative importance of bioclimatic factors, the most influential variable was precipitation seasonality (BIO15) (Figure 1) and (Table 2) for determining the habitat suitability for tea in Sri Lanka, which is comparable to the present study [23]. Annual mean temperature (BIO1) was the key indicator for determining the suitability of tea habitat in Kenya and China ( Figure 1) and (Table 2). Similar to the findings of previous studies, the highest probability of tea occurrence in the four countries was related to areas having a mean annual temperature around 20 • C [16,58]. Annual precipitation (BIO12) emerged as the strongest range predictor to determine current and future suitable habitats for tea in India and Kenya as precipitation is critical and plays an important role in sustaining tea cultivation [16].
Our study provides an indication for each country how their tea cultivation could be affected by future climate, with some countries more severely affected than others. Here, we have only considered climate suitability in the four major tea-producing countries in the world. Future research could expand this study in other major tea-producing countries. It is also important to include other non-climatic factors in assessing tea habitat suitability, as the actual areas available in the future are not only affected by climate factors; however, such data are currently not readily available for many countries or not available at suitable resolutions for modelling exercises. Inclusion of such data is likely to increase model accuracies. It should also be noted that land use was not included in this analysis. While there may be areas where climate is projected to be suitable, those lands may not actually be available for tea cultivation if they are currently being used for other higher-value cultivation or if they are centres of residential or commercial development.

Conclusions
We successfully modeled the current and future habitat suitability of tea in Sri Lanka, Kenya, India and China, which allowed us to identify potential distribution and spatial shifts of tea. In this context, the four tea-producing countries mapped for current and future 'suitability habitat' provided results that could be helpful in estimating the geographically based impacts of climate change. The predictions suggest that, by the years 2050 and 2070, climate would have an adverse impact on the habitat suitability of tea in Kenya, Sri Lanka and China, while India will gain an advantage. The current and future maps of habitat suitability and identification of suitable bioclimatic variables provides a significant guide for the tea sector across the globe to implement adequate adaptation strategies to circumvent the effects of climate change. Here, our results are based in part on the assumption that the potential area with suitable climatic conditions for tea also comprises other land uses. Thus, the resulting maps should be combined with the knowledge of local authorities of each country to identify potential new suitable habitats in order to expand tea cultivation. Our findings are a harbinger of what could happen in the four major tea-producing countries in general. Further analysis is required to identify land-use changes and decide the effective area of suitable lands using ensemble' approaches with other GCMs and inclusion of non-climatic factors.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/10/1536/s1, Figure S1: MaxEnt jackknife tests of the relative importance of bioclimatic variables with their regularized training gain for tea species in the four tea producing countries, Table S1: The comprehensive characteristics of each country selected in this study, Table S2: Occurrence data of tea for each country, Table S3: Confusion matrix, measures and formulas of predictive accuracy indicators for the developed model, Tables S4-S7: Areas (km 2 ) of suitable habitats for Camellia sinensis in Sri Lanka (S4), Kenya (S5), India (S6) and China (S7) under current and future climate conditions, Tables S8-S11: Changes in habitat suitability areas in Sri Lanka (S8), Kenya (S9), India (S10), and China (S11) [% of loss (−) or gain (+) ] for Camellia sinensis by 2050 and 2070 under RCP 2.6, 6.0, and 8.5 using GCMs of HadGEM2-ES, CCSM4 and MIROC5, Table S12: Projected average temperature (Avg. Temp.) and rainfall (avg. rainfall) of the three GCMs compared to historical mean.