Evaluation of Climate Change Impacts on the Potential Distribution of Styrax sumatrana in North Sumatra, Indonesia

: This study aims to assess the impact of climate change on the distribution of Styrax sumatrana in North Sumatra by applying the maximum entropy (MaxEnt) model with biophysical factors (elevation, slope, aspect, and soil), climatic factors (19 bioclimate data sets for 2050 and 2070), and anthropogenic factors (land use land cover (LULC) changes in 2050 and 2070). The future climate data retrieved and used are the output of four climate models from Coupled Model Intercomparison Project Phase 5 (CMIP5), namely, the CCSM4, CNRM-CM5, MIROC5, and MRI-CGCM3 models, under the Representative Concentration Pathways (RCPs) 4.5 and 8.5 scenarios. The MaxEnt modelling results showed the importance of the mean temperature of the coldest quarter and the LULC variables. Styrax sumatrana rely on environmental conditions with air temperatures ranging from 13 to 19 ◦ C. The potentially suitable land types for Styrax sumatrana are shrubs, gardens, and forests. The future predictions show that the suitable habitat for Styrax sumatrana is predicted to decrease to 3.87% in 2050 and to 3.54% in 2070 under the RCP4.5 scenario. Under the RCP8.5 scenario, the suitable area is predicted to decrease to 3.04% in 2050 and to 1.36% in 2070, respectively. The degradation of the suitable area is mainly due to increasing temperature and deforestation in future predictions. The modelling results illustrate that the suitable habitats of Styrax sumatrana are likely to be reduced under future climate change scenarios or lost in 2070 under the RCP8.5 scenario. The potential future extinction of this species should alert authorities to formulate conservation strategies. Results also demonstrated key variables that should be used for formulating ex situ conservation strategies.


Introduction
Indonesia is a country with massive forest biodiversity. More than 20,000 flora species are found in Indonesia, of which 40% are endemic. However, global changes, such as climate, land conversion, and natural disasters, threaten the tree species of Indonesia. Currently, 240 plant species are endangered, 36 of which are trees [1]. Recently, the impact of climate change has directly or indirectly affected tree species. For example, climate change significantly impacts tree growth. It has been shown that the circular rings of trees in Thailand decrease as a result of long-term climate change [2]. In China, tropical acacias that have been planted are destroyed by cold damage in inland areas and typhoon damage in coastal areas [3]. The pine species in Southeast Asia will face threats from increasing temperatures in lowland areas [4]. Thus, there is preliminary evidence indicating that climate change will impact tree species that spread around tropical areas, such as Indonesia, where the high biodiversity is threatened.
Among all of the regions in Indonesia, North Sumatra is one of the regions with a massive forest area that has been threatened by climate and land-use change [5]. North Sumatra has experienced large-scale development in the last three decades, and industrialization Species distribution modelling has also been applied in Indonesia, for example, to estimate bird distribution [24] and to predict the endangered tree species distribution in Borneo [25]. Moreover, alien invasive species distribution has been modelled in a national park using the MaxEnt and GLM approaches [24,26]. Furthermore, the model is applied to analyze the threat of invasive species in primate conservation by predicting anthropogenic risk from human intervention [27] and to formulate conservation strategies for proboscis monkeys in Kalimantan [28]. However, SDM applications are still limited in Indonesia, and most of the research uses SDM approaches, such as the GLM, MaxEnt and survey analysis, with limited environmental variables. Moreover, the impacts of climate change on species distribution have not been reported in any region of Indonesia by using SDMs.
The hypothesis in this study is straightforward: climate change will impact species distribution in the future. Therefore, this study aims to address the impacts of climate change on tree distribution and predict the future distribution of Styrax sumatrana by considering climate change and LULC change for 2050 and 2070. To understand the potential area suitable for Styrax sumatrana in the North Sumatra region, MaxEnt was applied in this study. The available presence locations and environmental data were collected and used in the MaxEnt modelling. Then, the future predictions of the Styrax sumatrana distribution in 2050 and 2070 were calculated by using the future climate projection results of four general circulation models (GCMs), namely, the CCSM4, CNRM-CM5, MIROC5, and MRI-CGCM3 models, from the coupled model intercomparison project phase 5 (CMIP5) data set under two RCP4.5 and RCP8.5 scenarios. Moreover, the LULCs for 2050 and 2070 are considered in the MaxEnt modelling by predicting the LULC changes with the artificial neural network-based cellular automata method.

Study Area
The study area covers the whole North Sumatra region, as shown in Figure 1. The region is situated between 0 • 34 6.39" S-4 • 18 16.59" N and 97 • 3 32.7"-100 • 25 27.85" E. The total area is approximately 7,106,868 ha. The location is surrounded by the Indian Ocean to the west and Malaysia to the east. Moreover, to the north, the region borders the Aceh region, and to the south, it borders Riau and West Sumatra. The North Sumatra region is divided into 26 regencies and 7 cities.
North Sumatra has temperature ranged from 13.4 to 33.9 • C with humidity around 78-91%. Annual precipitation notes around 800 to 4000 mm. Climate change phenomena has recorded in North Sumatra by the increasing of temperature and fluctuating precipitation in the last five years [5]. The altitude ranges from 0 to 2400 m above the sea level. The mountain slope is spread in the western part while declivous area was lengthways to the eastern part.
The total forest in North Sumatra covered approximately 42.9% of the total area in 2014. However, the forest cover has been gradually decreasing for years. Currently, the total amount of forest cover area is only 26.1% of the total area. The preservation area was 39.5% from the forest cover and 1.6% was conservation area [5]. The number is small rather than the forest cover that use for the log-forest product plantation.

Data Set
MaxEnt requires various data sets that represent the location of the species occurrence and the environmental information where the species is present (Phillips, Anderson and Schapire, 2006). The data sets consist of the species presence location, elevation, aspect, slope, soil, climate, and anthropological factors, and they are regrouped into five categories: the species distribution sample, biophysical factors, climate factors, anthropogenic factors, and future climate. The presence locations of the species are used as the species distribution sample input in the MaxEnt model. The elevation, aspect, slope, and soil are grouped as biophysical factors. The climate factors include 19 bioclimate variables collected from WorldClim [31,32]. The anthropogenic factor is represented by the LULC. Future climate data are from climate change projections for 2050 and 2070. Generally, the spatial data of environmental variables must be in raster format and have the same geographical coordinate system over the same spatial extent. The resolution used in this study is 30 m with a spatial extent of −0.568° S-4.305° N, 97.059°-100.424° E. The forecast data use the same spatial extent as the current environmental data. Each category is described in detail in the following subsection. Table 1 depicts a summary of the data set used in this study.  North Sumatra known as the region that produce Sumatra benzoin from the Styrax sumatrana trees which is important livelihood for the farmer [10]. Styrax sumatrana is mainly distributed in six regions, namely, Tapanuli Utara, Tapanuli Tengah, Toba Samosir, Humbang Hasundutan, Pakpak Barat, and Dairi [29]. However, small populations have been noted in two other regions: Karo and Simalungun [30].

Data Set
MaxEnt requires various data sets that represent the location of the species occurrence and the environmental information where the species is present (Phillips, Anderson and Schapire, 2006). The data sets consist of the species presence location, elevation, aspect, slope, soil, climate, and anthropological factors, and they are regrouped into five categories: the species distribution sample, biophysical factors, climate factors, anthropogenic factors, and future climate. The presence locations of the species are used as the species distribution sample input in the MaxEnt model. The elevation, aspect, slope, and soil are grouped as biophysical factors. The climate factors include 19 bioclimate variables collected from WorldClim [31,32]. The anthropogenic factor is represented by the LULC. Future climate data are from climate change projections for 2050 and 2070. Generally, the spatial data of environmental variables must be in raster format and have the same geographical coordinate system over the same spatial extent. The resolution used in this study is 30 m with a spatial extent of −0.568 • S-4.305 • N, 97.059 • -100.424 • E. The forecast data use the same spatial extent as the current environmental data. Each category is described in detail in the following subsection. Table 1 depicts a summary of the data set used in this study.  [30]. The locations were in a restricted area which was used as the preservation and conservation forest. Small number were located in farmers' garden and used to produce the benzoin resin as non-timber forest product. The total number of samples included 63 locations. However, it was reduced to 34 locations while re-gridding over a unit area (1 km 2 ). In this re-gridding process, several samples within a unit area were noted as one location on the map because they were too close to each other. The presence locations of Styrax sumatrana in North Sumatra region are presented in Figure 1b. Recently, the species distribution of Styrax sumatrana in North Sumatra was calculated using the MaxEnt model by Saputra et al. [33]. This research showed that the current distribution of Styrax sumatrana was approximately 8.91% or approximately 663,221.94 ha in the North Sumatra region ( Figure 2).  [30]. The locations were in a restricted area which was used as the preservation and conservation forest. Small number were located in farmers' garden and used to produce the benzoin resin as non-timber forest product. The total number of samples included 63 locations. However, it was reduced to 34 locations while re-gridding over a unit area (1 km 2 ). In this re-gridding process, several samples within a unit area were noted as one location on the map because they were too close to each other. The presence locations of Styrax sumatrana in North Sumatra region are presented in Figure 1b. Recently, the species distribution of Styrax sumatrana in North Sumatra was calculated using the MaxEnt model by Saputra et al. [33]. This research showed that the current distribution of Styrax sumatrana was approximately 8.91% or approximately 663,221.94 ha in the North Sumatra region ( Figure 2).

Biophysical Factors
Biophysical factors are mostly derived from the Advanced Space-borne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) and Food and Agriculture Organization (FAO). The ASTER GDEM is one of the most widely used DEM data sets. The elevation was obtained from the DEM. The ASTER GDEM V2 data set has a 30-m spatial resolution in the GeoTIFF image format with decimal degrees and WGS84 datum [34]. To cover the North Sumatra region, 18 tiles of the ASTER DEM were downloaded. The aspect and slope of the study region were also estimated from the DEM by applying spatial analysis. The soil map and soil type were retrieved from the FAO and refer to the FAO soil classification [35] at a 30-m resolution. The biophysical factors as the raster map are displayed in Figure 3.  Tree distributions are sensitive to climate, and small changes in climate can affect the morphology of trees [36]. Quantifying the effects of climate change on past, present, and future species distributions can be achieved by examining climate over time. To evaluate the climate over time, the USGS provides bioclimatic predictors that highlight climate conditions most related to species physiology [32]. The bioclimatic predictors consist of 19 bioclimate variables, such as annual trends (e.g., the mean annual temperature, annual precipitation), seasonality (e.g., the annual range in temperature and precipitation), and extreme or limiting environmental factors (e.g., the temperature of the coldest and warmest months and the precipitation of wet and dry quarters). These bioclimatic predictors

Climate Factors
Tree distributions are sensitive to climate, and small changes in climate can affect the morphology of trees [36]. Quantifying the effects of climate change on past, present, and future species distributions can be achieved by examining climate over time. To evaluate the climate over time, the USGS provides bioclimatic predictors that highlight climate conditions most related to species physiology [32]. The bioclimatic predictors consist of 19 bioclimate variables, such as annual trends (e.g., the mean annual temperature, annual precipitation), seasonality (e.g., the annual range in temperature and precipitation), and extreme or limiting environmental factors (e.g., the temperature of the coldest and warmest months and the precipitation of wet and dry quarters). These bioclimatic predictors are commonly used as climate factors in SDM analysis [13,32,37]. The present climate factors are considered over the period 1950-2000 and were retrieved from WorldClim data set [31]. The list of 19 bioclimate variables and their descriptions are presented in Table 2.

Anthropogenic Factors
The LULC of the study region is considered an anthropogenic factor. The changes in land cover can be interpreted as human interactions with their environment [38]. Land-use change also affects the distribution of species due to development and human interest [20,39]. The LULC provides information on the land types where species are located. This information allows SDMs to estimate the distribution of species during changes in land use [20].
Saputra and Lee [40] simulated and predicted the changes in LULC in North Sumatra, Indonesia, by using an artificial neural network-based cellular automation (ANN-CA) model. In their prediction, five criteria (altitude, slope, aspect, distance from the road, and soil type) were considered exploratory data in the learning process of the ANN to determine their effects on LULC changes between 1990 and 2000. After the learning process, the 2010 LULC was predicted and compared with the real satellite-driven 2010 LULC, resulting in high agreement, with a Kappa index of 0.83 and a percentage of correctness of 87.28%. Then, the ANN-CA model was applied to predict LULC changes in 2050 and 2070. The LULC predictions for 2050 and 2070 demonstrated high increments in plantation area of more than 4%. Meanwhile, forest and crop areas were projected to decrease by approximately 1.2% and 1.6%, respectively, by 2050. The forest and crop areas were further reduced by 1.2% and 1.7%, respectively, by 2070, implying humans influence LULC changes, such as a change from forest and cropland to plantation. Saputra and Lee [40] showed that the prediction of LULC changes using the ANN-CA model produced reliable results. The predicted North Sumatra LULC data sets for 2050 and 2070 by Saputra and Lee [40], as shown in Figure 4, were utilized in the species distribution modelling for Styrax sumatrana.
Sustainability 2021, 13, x FOR PEER REVIEW 8 of 23 [40] showed that the prediction of LULC changes using the ANN-CA model produced reliable results. The predicted North Sumatra LULC data sets for 2050 and 2070 by Saputra and Lee [40], as shown in Figure 4, were utilized in the species distribution modelling for Styrax sumatrana.

Future Climate and Scenarios
The prediction of species distribution in 2050 (average of predictions for 2041-2060) and 2070 (average of predictions for 2061-2080) was conducted by using future climate data under the two representative concentration pathways (RCPs) for carbon dioxide, i.e., the RCP4.5 and RCP8.5 scenarios [41]. The RCP4.5 scenario indicates an atmospheric mean temperature increase of 1.4 °C and 1.8 °C in 2050 and 2070, respectively, while RCP8.5 projects a temperature increase of 2 °C in 2050 and 3.7 °C in 2070. The temperature of the Earth is predicted to increase by 1.4-5.8 °C from 1990 to 2100, and precipitation is estimated to increase by up to 1.0% in mid-and high-latitude regions and by 0.3% in tropical zones [42]. The climate factors for future prediction utilized the climate projection results from four GCMs in the CMIP5 archive. The four GCMs used were CCSM4, CNRM-

Future Climate and Scenarios
The prediction of species distribution in 2050 (average of predictions for 2041-2060) and 2070 (average of predictions for 2061-2080) was conducted by using future climate data under the two representative concentration pathways (RCPs) for carbon dioxide, i.e., the RCP4.5 and RCP8.5 scenarios [41]. The RCP4.5 scenario indicates an atmospheric mean temperature increase of 1.4 • C and 1.8 • C in 2050 and 2070, respectively, while RCP8.5 projects a temperature increase of 2 • C in 2050 and 3.7 • C in 2070. The temperature of the Earth is predicted to increase by 1.4-5.8 • C from 1990 to 2100, and precipitation is estimated to increase by up to 1.0% in mid-and high-latitude regions and by 0.3% in tropical zones [42]. The climate factors for future prediction utilized the climate projection results from four GCMs in the CMIP5 archive. The four GCMs used were CCSM4, CNRM-CM5, MIROC5 and MRI-CGCM3 based on their higher spatial resolutions among the CMIP5 GCMs [42].
MaxEnt is a machine learning method that uses the presence location as the input for analysis and collects information related to the environmental conditions at the exact point where the location of the sample exists [43]. The goal of MaxEnt is to estimate a target probability distribution by finding the probability distribution of maximum entropy, subject to a set of constraints that represent our incomplete information about the target distribution [18]. The perfect formula for SDMs with both presence and absence data relies on Bayes rule [37] and can be represented as follows: where P(y = 1|x) is the probability that the species is present at site x (y ranges from 0 to 1), P(x|y = 1) is the present observation or realized distribution at area x annotated as π(x), P(y = 1) is the probability of presence, and P(x) = 1/|X| is the all-area probability from site X. Similarly, However, MaxEnt uses presence-only data, which means it needs to be transferred into a spatial distribution first to calculate the probability of species presence [18]. Regarding the spatial information that needs the distribution form, Bayes estimation rule cannot be implemented. Rather than estimating Bayes estimation rule, MaxEnt uses the Gibbs distribution from the sets of features f 1 , . . . , f n and weights λ = (λ 1 , . . . ., λ n ). Features or predictors are the information available about target distributions (precipitation, temperature, etc.) or environmental variables. The weights or constants are the expected value of each feature that should match its empirical average (the average values of each sample point taken from the target distribution). The value of the weight by default depends on the features and number of samples [18,44]. This distribution calculated the realized distribution (π(x)) from Bayes rule. The formula was defined as follows [37]: is an exponential distribution parameterized by a vector of feature (f ) and weight (λ), and Z λ is a normalization constant ensuring that the values of q λ (x) sum to one over the area.
Since the MaxEnt objective is to estimate the probability of presence by maximum entropy, we include the formula delivered from information theory by Shanon [37]. The formula is computed as follows: where H is the maximum entropy, and q λ (x) is a Maxent distribution from Equation (3).
After obtaining an estimate of q λ , sufficient information is gained to calculate the probability distribution P(y = 1|x), as indicated by Equation (5) where q λ is the estimated probability of presence with the maximum entropy of π, and H is the entropy of q λ . The probability distribution is arranged from 0 to 1. A higher number means there is a higher probability the species will occur. The probability number was classified into five classes based on FAO classification [45]. This classification in Table 3 was used to underline the suitable and unsuitable areas for Styrax sumatrana distribution. The areas with a probability greater than or equal to 0.4 are considered as suitable and others are unsuitable. The maximum number of iterations was set to 500 to allow the model to have adequate time for convergence. The convergence threshold was 1 × 10 −6 . The maximum number of background points for sampling was kept at 10,000. The logistic output from the model indicated the probabilities (0-1) but was scaled up in a nonlinear way for easier interpretation.

Baseline Condition
To compare the current and future changes in tree species distributions and to investigate the climate change effects on the future distributions, the baseline condition was defined since the data set was collected from several different timelines. The presence location data set of Styrax sumatrana was surveyed from different years: 2000,2004,2010,2012, and 2015. The information was gathered by the Environment and Forestry Research and Development Agency of the Aek Nauli team. The result consisted of 34 location data with biophysical and environmental condition information. The elevation and soil type maps were collected from the data for 2000. The aspect and slope variables were made based on the DEM with spatial analysis using ArcGIS 10.1 (Esri, Redlands, California, USA). Meanwhile, the 19 bioclimate data were downloaded from WorldClim version 2.1 with an average value from 1997 to 2000 [31]. The 2000 LULC was used for the baseline condition.

Future Prediction
Future predictions were mainly based on LULC changes and future climate conditions. Additionally, biophysical factors such as altitude, slope, aspect, and soil type used the same data sets as the baseline condition. The LULC prediction maps in North Sumatra were obtained from Saputra and Lee [40] and utilized as the input for future predictions as described. The input for climate conditions was collected from 4 GCMs under two RCP scenarios, 4.5 and 8.5. Most of the studies of climate change impacts with SDMs used climate conditions from fewer than four climate models [4,13,46,47]. The four GCMs were chosen to encompass the full range of variation in the models in the multimodel ensemble CMIP5 [42]. Therefore, a total of 16 MaxEnt models by 2 RCP scenarios and 4 GCM combinations for 2050 and 2070 were conducted for future predictions. The overall MaxEnt modelling process is shown in Figure 5.

Validation
The MaxEnt models were evaluated using the area under the curve (AUC), whose values were calculated from the receiver operating characteristic (ROC) curve. The ROC curve is a graph showing the performance of a classification model at all classification thresholds [43]. The ROC curve consists of sensitivity on the y axis and 1-specificity on the x axis for all possible thresholds [18]. Sensitivity explains how well the model predicts presence, whereas specificity explains how well the model predicts absences or, in this case, pseudoabsence [48]. The 1-specificity explains the probability of the species presence, while in practice, there is absence. When the sensitivity is near 1 and the specificity is near 0, then the AUC has the highest rank in the ROC graph.
The AUC value estimates the probability of the randomly selected grid cell in a correctly adjusted model and defines the success of the model with all possible thresholds. If the AUC is greater than 0.5, the model performs better than a random estimation, and the better model is presented by the higher AUC in the ROC graph. If the AUC value becomes higher and closer to 1, then the separation is better, and the model is more sensitive and descriptive [18].

MaxEnt Modeling
With respect to the input number of presence-only sample locations, sensitivity tests were performed with 63 locations and 34 locations. The test results showed that the MaxEnt results with 63 and 34 sample locations did not present significant differences. Therefore, the MaxEnt results with 34 presence-only locations are discussed hereinafter.

Validation
The MaxEnt models were evaluated using the area under the curve (AUC), whose values were calculated from the receiver operating characteristic (ROC) curve. The ROC curve is a graph showing the performance of a classification model at all classification thresholds [43]. The ROC curve consists of sensitivity on the y axis and 1-specificity on the x axis for all possible thresholds [18]. Sensitivity explains how well the model predicts presence, whereas specificity explains how well the model predicts absences or, in this case, pseudoabsence [48]. The 1-specificity explains the probability of the species presence, while in practice, there is absence. When the sensitivity is near 1 and the specificity is near 0, then the AUC has the highest rank in the ROC graph.
The AUC value estimates the probability of the randomly selected grid cell in a correctly adjusted model and defines the success of the model with all possible thresholds. If the AUC is greater than 0.5, the model performs better than a random estimation, and the better model is presented by the higher AUC in the ROC graph. If the AUC value becomes higher and closer to 1, then the separation is better, and the model is more sensitive and descriptive [18].

MaxEnt Modeling
With respect to the input number of presence-only sample locations, sensitivity tests were performed with 63 locations and 34 locations. The test results showed that the MaxEnt results with 63 and 34 sample locations did not present significant differences. Therefore, the MaxEnt results with 34 presence-only locations are discussed hereinafter.
MaxEnt output provides the result of the AUC for validating the performance of the model. In Figure 6, the MaxEnt result shows an AUC value of 0.949, which is higher than 0.9 and thus classified as the best model performance [18]. MaxEnt output provides the result of the AUC for validating the performance of the model. In Figure 6, the MaxEnt result shows an AUC value of 0.949, which is higher than 0.9 and thus classified as the best model performance [18]. In the baseline condition (2000) result in Figure 2, the percentage of suitable area for Styrax sumatrana was approximately 8.91% of the total North Sumatra region. The contribution rates of each input variable for the baseline condition are presented in Table 2. The three main variables with the highest contribution rates were the mean temperature in the coldest quarter (28.5%), the altitude (26.3%), and the LULC (15.4%). The changes in the mean temperature of the coldest quarter and the LULC variables in the future projection model mainly determined the species distribution in 2050 and 2070 for each scenario since the altitude was considered constant. The effects of each variable in the prediction of future projection are described in the following subsections.

Effect of the Mean Temperature of the Coldest Quarter
The mean temperature of the coldest quarter, shown in Figure 7, provides the seasonal distribution of Styrax sumatrana. The distribution of the mean temperature of the coldest quarter in the baseline condition is indicated in Figure 7a, where the temperature in North Sumatra ranged from 11.9 to 26.6 °C. The presence location shown by the red colour has a temperature of approximately 13-23 °C. Figure 7b indicates the Styrax sumatrana requirement at low temperature. The range of the low temperature was from 13 to 23 °C. The graph shows the optimum mean temperature was between 13 and 19 °C. The probability of Styrax sumatrana distribution will decrease when the temperature surpasses the peak (approximately 19 to 26 °C). In the baseline condition (2000) result in Figure 2, the percentage of suitable area for Styrax sumatrana was approximately 8.91% of the total North Sumatra region. The contribution rates of each input variable for the baseline condition are presented in Table 2. The three main variables with the highest contribution rates were the mean temperature in the coldest quarter (28.5%), the altitude (26.3%), and the LULC (15.4%). The changes in the mean temperature of the coldest quarter and the LULC variables in the future projection model mainly determined the species distribution in 2050 and 2070 for each scenario since the altitude was considered constant. The effects of each variable in the prediction of future projection are described in the following subsections.

Effect of the Mean Temperature of the Coldest Quarter
The mean temperature of the coldest quarter, shown in Figure 7, provides the seasonal distribution of Styrax sumatrana. The distribution of the mean temperature of the coldest quarter in the baseline condition is indicated in Figure 7a, where the temperature in North Sumatra ranged from 11.9 to 26.6 • C. The presence location shown by the red colour has a temperature of approximately 13-23 • C. Figure 7b indicates the Styrax sumatrana requirement at low temperature. The range of the low temperature was from 13 to 23 • C. The graph shows the optimum mean temperature was between 13 and 19 • C. The probability of Styrax sumatrana distribution will decrease when the temperature surpasses the peak (approximately 19 to 26 • C).  The effect of the mean temperature of the coldest quarter for the baseline condition is shown in Figure 8a. Three different temperature ranges were notably distributed in North Sumatra. The green colour represents a suitable temperature for Styrax sumatrana distribution within 13-19 °C. The yellow colour indicates a suitable temperature, yet it gradually declines from 19 to 23 °C. The red colour presents the unsuitable temperature as the Styrax sumatrana variable required to survive.
The mean temperature of the coldest quarter in 2050 under the RCP4.5 scenario from the CNRM-CM5 model is illustrated in Figure 8b. The prediction uses a moderate scenario that presented an increment in temperature of approximately 1 °C. In this state, the green colour, which shows a suitable temperature for Styrax sumatrana, was decreased from the baseline condition. The green area will change into a yellow area, illustrating the increment of area for the mean temperature of the coldest quarter.
The mean temperature from the same model with the projection in 2070 in the same scenario is demonstrated in Figure 8c as that in Figure 8b. Based on the timeline in the     The effect of the mean temperature of the coldest quarter for the baseline condition is shown in Figure 8a. Three different temperature ranges were notably distributed in North Sumatra. The green colour represents a suitable temperature for Styrax sumatrana distribution within 13-19 • C. The yellow colour indicates a suitable temperature, yet it gradually declines from 19 to 23 • C. The red colour presents the unsuitable temperature as the Styrax sumatrana variable required to survive.
The mean temperature of the coldest quarter in 2050 under the RCP4.5 scenario from the CNRM-CM5 model is illustrated in Figure 8b. The prediction uses a moderate scenario that presented an increment in temperature of approximately 1 • C. In this state, the green colour, which shows a suitable temperature for Styrax sumatrana, was decreased from the baseline condition. The green area will change into a yellow area, illustrating the increment of area for the mean temperature of the coldest quarter.
The mean temperature from the same model with the projection in 2070 in the same scenario is demonstrated in Figure 8c as that in Figure 8b. Based on the timeline in the same scenario, the mean temperature of the coldest quarter will increase, and the suitable mean temperature of the coldest quarter in which Styrax sumatrana can survive will decrease.
The mean temperature of the coldest quarter in 2050 under the RCP8.5 scenario is presented in Figure 8d, where the temperature increases almost twice that under the RCP4.5 scenario. The area with the mean temperature of the coldest quarter of approximately 13-18 • C will decrease dramatically from the baseline condition. The conditions in the RCP8.5 scenario in 2050 are similar to those in the RCP4.5 scenario in 2070. The RCP8.5 scenario accelerated the increase in temperature from the 2050 conditions to the 2070 conditions. The mean temperature of the coldest quarter in 2070 under RCP8.5 is shown in Figure 8e. In this year, the suitable temperatures under which Styrax sumatrana can survive were lost, and only small scales in the northern North Sumatra region remained. The other three GCMs depicted the same trend as CNRM-CM5. The mean temperature of the coldest quarter showed the shifting of the temperature from the baseline conditions to 2050 and 2070. The low temperature declined while the higher temperature spread across the region. The RCP8.5 scenario indicated very large changes from RCP4.5 such that the total area with a higher temperature was twice that of RCP4.5. The condition causes Styrax sumatrana to go extinct due to global warming.
The jackknife test was applied to determine the significance of environmental variables in the creation of the model. Based on the jackknife analysis in Figure 9, the mean temperature of the coldest quarter variable could explain the potential distribution of Styrax sumatrana without considering the other variables. This variable had the highest impact on the model when the model used only the mean temperature of the coldest quarter variable. Figure 9 shows the gain for the model when using only variables greater than 0.8. This value was higher than another variable under the same conditions. This result indicates that the mean temperature of the coldest quarter provides the highest contribution to the modelling result. same scenario, the mean temperature of the coldest quarter will increase, and the suitable mean temperature of the coldest quarter in which Styrax sumatrana can survive will decrease.
The mean temperature of the coldest quarter in 2050 under the RCP8.5 scenario is presented in Figure 8d, where the temperature increases almost twice that under the RCP4.5 scenario. The area with the mean temperature of the coldest quarter of approximately 13-18 °C will decrease dramatically from the baseline condition. The conditions in the RCP8.5 scenario in 2050 are similar to those in the RCP4.5 scenario in 2070. The RCP8.5 scenario accelerated the increase in temperature from the 2050 conditions to the 2070 conditions.
The mean temperature of the coldest quarter in 2070 under RCP8.5 is shown in Figure  8e. In this year, the suitable temperatures under which Styrax sumatrana can survive were lost, and only small scales in the northern North Sumatra region remained. The other three GCMs depicted the same trend as CNRM-CM5. The mean temperature of the coldest quarter showed the shifting of the temperature from the baseline conditions to 2050 and 2070. The low temperature declined while the higher temperature spread across the region. The RCP8.5 scenario indicated very large changes from RCP4.5 such that the total area with a higher temperature was twice that of RCP4.5. The condition causes Styrax sumatrana to go extinct due to global warming.
The jackknife test was applied to determine the significance of environmental variables in the creation of the model. Based on the jackknife analysis in Figure 9, the mean temperature of the coldest quarter variable could explain the potential distribution of Styrax sumatrana without considering the other variables. This variable had the highest impact on the model when the model used only the mean temperature of the coldest quarter variable. Figure 9 shows the gain for the model when using only variables greater than 0.8. This value was higher than another variable under the same conditions. This result indicates that the mean temperature of the coldest quarter provides the highest contribution to the modelling result.

Effect of Altitude
Elevation contributed 26.3% of all variables in the Styrax sumatrana distribution. Most Styrax gardens are found above 600 m in North Sumatra. High altitude is known to be suitable for the phenotypic performance for Styrax sumatrana growth [49]. In the projection of this study, the altitude was considered constant in 2050 and 2070. Although it exists at a wide range of elevations in North Sumatra, the optimum high Styrax sumatrana distribution ranged only from 600 to 1400 m.

Effect of LULC
The jackknife test results in Figure 9 show that the bar chart of LULC was less than the other variables when the LULC was omitted from the model. The model would not be

Effect of Altitude
Elevation contributed 26.3% of all variables in the Styrax sumatrana distribution. Most Styrax gardens are found above 600 m in North Sumatra. High altitude is known to be suitable for the phenotypic performance for Styrax sumatrana growth [49]. In the projection of this study, the altitude was considered constant in 2050 and 2070. Although it exists at a wide range of elevations in North Sumatra, the optimum high Styrax sumatrana distribution ranged only from 600 to 1400 m.

Effect of LULC
The jackknife test results in Figure 9 show that the bar chart of LULC was less than the other variables when the LULC was omitted from the model. The model would not be reliable if the LULC variable was not considered in the model. This result means that there was a value of the LULC that was not mentioned in the other variables. Thus, LULC is important in the model and increases the accuracy of the model. LULC category distribution in the baseline condition in Figure 10a shows that crops and forests have higher portions of the area in North Sumatra among all categories. However, the species was mostly found in four types of LULC categories, namely, shrubs, gardens, forests, and crops. ever, the species was mostly found in four types of LULC categories, namely, shrubs, gardens, forests, and crops.
In the response chart of the LULC in Figure 10b, a high probability was found in shrubs and gardens followed by forests. The probability of shrubs was the highest, with a value greater than 0.73. Meanwhile, gardens had a probability of only approximately 0.62, followed by forests with a probability of approximately 0.46. These three categories were classified into suitable areas for Styrax sumatrana.  In the response chart of the LULC in Figure 10b, a high probability was found in shrubs and gardens followed by forests. The probability of shrubs was the highest, with a value greater than 0.73. Meanwhile, gardens had a probability of only approximately 0.62, followed by forests with a probability of approximately 0.46. These three categories were classified into suitable areas for Styrax sumatrana.
Future climate prediction data were collected from four GCMs. Each model presented the different values of climate variables for both scenarios and times. The comparison of the suitability among the 4 GCMs of the climatic factor for the different years under the RCP4.5 and RCP8.5 scenarios is illustrated in Figure 11.
Future climate prediction data were collected from four GCMs. Each model presented the different values of climate variables for both scenarios and times. The comparison of the suitability among the 4 GCMs of the climatic factor for the different years under the RCP4.5 and RCP8.5 scenarios is illustrated in Figure 11. The four models indicated the area changes compared to the baseline condition. All models showed the degradation of the future potential area suitable for Styrax sumatrana ( Figure 11). The MRI-CGCM3 model predicted the highest change in suitable area for Styrax sumatrana between the projection and the baseline condition. The CCSM4 model predicted the second highest changes. Meanwhile, a small number of changes were presented by the CNRM-CM5 and MIROC5 models. The one of the main differences among those models was in the atmospheric grid resolution. CCSM4 has the smallest grid by 0.94 latitude and 1.25 longitude. Moreover, MRI-CGCM3 has the second smallest grid by 1.12 × 1.12 in latitude and longitude, respectively. CNRM-CM5 and MIROC5 have higher grid values of 1.4 × 1.4 in latitude and longitude, respectively. The grid resolution represents the same temperature and precipitation value for the climate model in one certain area. The four models indicated the area changes compared to the baseline condition. All models showed the degradation of the future potential area suitable for Styrax sumatrana ( Figure 11). The MRI-CGCM3 model predicted the highest change in suitable area for Styrax sumatrana between the projection and the baseline condition. The CCSM4 model predicted the second highest changes. Meanwhile, a small number of changes were presented by the CNRM-CM5 and MIROC5 models. The one of the main differences among those models was in the atmospheric grid resolution. CCSM4 has the smallest grid by 0.94 latitude and 1.25 longitude. Moreover, MRI-CGCM3 has the second smallest grid by 1.12 × 1.12 in latitude and longitude, respectively. CNRM-CM5 and MIROC5 have higher grid values of 1.4 × 1.4 in latitude and longitude, respectively. The grid resolution represents the same temperature and precipitation value for the climate model in one certain area. The average value of these different models was applied to create a higher validation rather than using the separate result from each model. This approach minimizes the bias from different grid intervals of the models and provides precise results.
However, the value of the model from the MIROC5 climate model for the RCP4.5 scenario in 2070 showed the opposite result as that of the other models. MIROC5 had a positive value for the potential area that was suitable for the species in 2070 under the RCP4.5 scenario. In MIROC5, the prediction in 2070 showed a more optimistic value for the mean temperature of the coldest quarter than did the other models. In MIROC5, the area with a high temperature in 2070 was smaller than that in the other models for the same year and scenario. Even though the value was positive with MIROC5 climate input, the averaged values from all models showed the degradation of the Styrax sumatrana potential distribution area. Thus, overlay analysis was used with the same weight for all the model results.
The MaxEnt model predicted the average distribution probability of Styrax sumatrana from four GCMs (Figure 12). The potential area decreased from 2050 to 2070 in scenarios RCP4.5 and RCP8.5. The higher temperature in the future prediction will decrease the sufficient area fulfilling the Styrax sumatrana requirement for growth. The potential area of Styrax sumatrana is presented in Table 4. This result illustrated that by 2050 and 2070, the total suitable area for Styrax sumatrana will decline. In the RCP8.5 scenario, the model showed that Styrax sumatrana will be almost extinct with a suitable area of 1.36%. The average value of these different models was applied to create a higher validation rather than using the separate result from each model. This approach minimizes the bias from different grid intervals of the models and provides precise results. However, the value of the model from the MIROC5 climate model for the RCP4.5 scenario in 2070 showed the opposite result as that of the other models. MIROC5 had a positive value for the potential area that was suitable for the species in 2070 under the RCP4.5 scenario. In MIROC5, the prediction in 2070 showed a more optimistic value for the mean temperature of the coldest quarter than did the other models. In MIROC5, the area with a high temperature in 2070 was smaller than that in the other models for the same year and scenario. Even though the value was positive with MIROC5 climate input, the averaged values from all models showed the degradation of the Styrax sumatrana potential distribution area. Thus, overlay analysis was used with the same weight for all the model results.
The MaxEnt model predicted the average distribution probability of Styrax sumatrana from four GCMs (Figure 12). The potential area decreased from 2050 to 2070 in scenarios RCP4.5 and RCP8.5. The higher temperature in the future prediction will decrease the sufficient area fulfilling the Styrax sumatrana requirement for growth. The potential area of Styrax sumatrana is presented in Table 4. This result illustrated that by 2050 and 2070, the total suitable area for Styrax sumatrana will decline. In the RCP8.5 scenario, the model showed that Styrax sumatrana will be almost extinct with a suitable area of 1.36%.  Climate change has a great impact on Styrax sumatrana distribution with increasing temperature. When the Styrax sumatrana requirement was the mean temperature of the coldest quarter between 13 and 19 °C, the future prediction from the four climate models showed that the temperature would increase and be less suitable for the species to survive.

Discussion
The increase in temperature in future scenarios indicates that the distribution of Sty-  Climate change has a great impact on Styrax sumatrana distribution with increasing temperature. When the Styrax sumatrana requirement was the mean temperature of the coldest quarter between 13 and 19 • C, the future prediction from the four climate models showed that the temperature would increase and be less suitable for the species to survive.

Discussion
The increase in temperature in future scenarios indicates that the distribution of Styrax sumatrana in North Sumatra will change. A small number of suitable areas that meet the survival requirements of the species in the future projection will remain in North Sumatra. The model predicts the possible area in the region that satisfies the Styrax sumatrana distribution requirement.
As shown in Figure 12, most of the possible area for Styrax sumatrana distribution will be less than that in the baseline condition. In RCP4.5, the species will decline by more than half compared to the baseline. In 2070, under RCP8.5, the suitable area is predicted to completely disappear. Thus, the results from different scenarios imply that the different countermeasures are necessary for preventing the suitable area from being reduced in future predictions.
Two main variables that contributed to the Styrax sumatrana distribution in this study were the mean temperature of the coldest quarter and the LULC change. Tree species cannot survive in areas with high temperatures outside of the shrub, forest, and garden in LULC categories.
The jackknife analysis indicated the importance of the mean temperature of the coldest quarter compared to the other variables. The temperature changes determined the suitable area for Styrax sumatrana. The CNRM-CM5 model predicted a smaller area that had the required temperature for Styrax sumatrana in the future than in the baseline condition. The baseline condition has 30% of total area that suitable for Styrax sumatrana, while only 16-24% of total area that satisfied the temperature required for the species in the future.
The predicted LULC changes in the future contributed to the degradation of suitable area. This LULC change from forest to open area, for example, increased the air temperature in the future and decreased land allocation for tree species, including Styrax sumatrana. The LULC change will degrade the required area for Styrax sumatrana and limit the area affected by human intervention. The land-use change due to grazing pasture increases land degradation in tropical forests around 0.11% and in secondary vegetation around 0.19%. The same situation occurred where forests and gardens with heterogenic species were converted into crop or plantation areas that were more monogenic [20].
Styrax sumatrana produces resin in granular or chunky shapes. The product of this species supplies a prominent income for farmers of approximately 47.64-79.23% [9,10,50]. The reduction in Styrax sumatrana habitat will also decrease the production of resin and hence the income of farmers. Thus, farmers will change Styrax gardens into other LULC types for a more valuable commodity such as rubber, coffee, or others that provide seasonal income. This LULC conversion will negatively affect the sustainability of Styrax sumatrana in its current locations. Therefore, resin production will decline in the market. Furthermore, Styrax sumatrana may disappear in the future, as predicted in 2070 under RCP8.5. Furthermore, such changes in LULC from economic activities by farmers might cause regional or local warming that could accelerate the degradation of the Styrax sumatrana distribution.

Conclusions
Ecological niche modellings were carried out to evaluate the impacts of climate change on spatial distribution of Styrax sumatrana in North Sumatra, Indonesia, by using MaxEnt with biophysical factors, anthropogenic factors, and climatic factors in 2050 and 2070. Elevation, slope, aspect, and soil types are used as biophysical factors and considered constant in target years. LULC is considered as anthropogenic factor, and its changes in 2050 and 2070 are also addressed in MaxEnt modelling. For future climatic factors in the modellings, four GCMs results from CMIP6 dataset are considered for 2050 and 2070 under two RCP 4.5 and 8.5 scenarios.
The modelling results indicate that the mean temperature of the coldest quarter, altitude, and LULC types are the most sensitive factors for the presence of Styrax sumatrana. Styrax sumatrana requires low temperatures (13-19 • C) in the highland to survive in North Sumatrana. This species is mostly suitable for habitation in shrub, garden, and forest types. The results with future climate conditions demonstrated that the suitable areas for Styrax sumatrana distribution in 2050 and 2070 for both the RCP4.5 and RCP8.5 scenarios is predicted to decrease. In particular, the suitable area in 2070 under RCP 8.5 was predicted to almost disappear in North Sumatra.
The warming climate and LULC change potentially increase the probability of Styrax sumatrana extinction in 2050 and 2070 in the RCP8.5 scenario. These factors will also impact farmers, as their income will decrease when the suitable area for Styrax sumatrana declines. Furthermore, the potential areas will be reduced as LULC types change into categories that provide more economic value than Styrax sumatrana.
Conservation strategies are needed to prevent for the species from extinction. The potential suitable area for Styrax sumatrana will be limited. However, ex situ conservation strategy can be implemented outside its natural habitat. This approach preserves the genetic resources of the species to prevent extinction in the presence area. Results of this study reveal the significant variables for species distribution and survival that are also important as the input in decision-making process for in situ and ex situ conservation strategies.
Finally, some limitations of this study, such as limited number of GCMs input for future climatic factors, limited numbers of presence locations of the species and spatiotemporal changes of LULC in baseline condition, have to be carefully considered in the interpretation of the results and to be improved in future works.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.