Predicting the Habitat Suitability of Melaleuca cajuputi Based on the MaxEnt Species Distribution Model

Abstract

Gelam tree or Melaleuca cajuputi (M. cajuputi) is an important species for the local economy as well as coastal ecosystem protection in Terengganu, Malaysia. This study aimed at producing a current habitat suitability map and predicting future potential habitat distribution for M. cajuputi in Terengganu based on Species distribution modeling (SDM) using the Maximum Entropy principle. Our modeling results show that for the current potential distribution of M. cajuputi species, only 10.82% (1346.5 km²) of Terengganu area is suitable habitat, which 0.96% of the areas are under high, 2.44% moderate and 7.42% poor habitat suitability. The model prediction for future projection shows that the habitat suitability for M. cajuputi would decrease significantly in the year 2050 under RCP 4.5 where the largest contraction from suitable to unsuitable habitat area is about 442.1 km² and under RCP 2.6 is the highest expansion from unsuitable to suitable habitat area (267.5 km²). From the future habitat suitability projection, we found that the habitat suitability in Marang would degrade significantly under all climate scenarios, while in Setiu the habitat suitability for M. cajuputi remains stable throughout the climate change scenarios. The modeling prediction shows a significant influence on the soil properties, temperature, and precipitation during monsoon months. These results could benefit conservationist and policymakers for decision making. The present model could also give a perception into potential habitat suitability of M. cajuputi in the future and to improve our understanding of the species’ response under the changing climate.

1. Introduction

Global warming has caused significant changes in spatial and temporal environmental patterns, and these changes also affect effort to conserve biodiversity [1]. In tropical climates, changes in temperature and precipitation can affect species habitat and plant phenology, both directly and indirectly [2]. Climate change effect also threatened Melaleuca swamp forest where the sea level rise causes saltwater intrusion and indirectly affects the Melaleuca swamp forest because of an increase in salinity [3]. Melaleuca genus is from the family of Myrtaceae that comprises 260 species, naturally distributed in Australia, Indonesia, Malaysia, Singapore, Thailand and Vietnam [4]. Melaleuca forests grow in coastal regions such as wetlands, lowlands, and peat lands, often behind the mangrove zone where pure or mixed stands may form. In Peninsular Malaysia, Melaleuca cajuputi (M. cajuputi) or known as white-paperbark tree (kayu putih) or Gelam, is the native species in Terengganu and Selangor states, but the most prominent coverage of M. cajuputi in Peninsular Malaysia is at Setiu Wetlands, Terengganu [5]. The M. cajuputi forest is one of the crucial habitats of the wetlands’ ecosystem due to its function and characteristics that can withstand flooding and dry conditions. Moreover, M. cajuputi species is a multifunction tree where the stem, leaves, and flowers are useful for the local economy [6]. However, due to anthropogenic activities, the M. cajuputi forest is degrading or converted to other land use [7]. The wide variation in botanical characteristics shows the ability of the Melaleuca genus to respond to diverse environments in a different location [4]. It was shown that the Melaleuca genus can tolerate various types of extreme conditions, whether in a wet, deep flooding area and dry land, but only some species of Melaleuca genus may adapt to increased inundation while others might not survive [8,9]. Therefore, finding out the responses between species and their environmental variables is crucial to predicting their habitat distribution.

In Malaysia, little is known about the phenology and spatio-temporal distribution of the M. cajuputi species [5]. The M. cajuputi is described as a highly resilient plant [10,11] due to its ability to adapt to various environmental conditions. However, the environmental and climate change impact on M. cajuputi species is yet to be examined, especially on its habitat sustainability and its potential to mitigate climate change through carbon stock accumulation [12]. Moreover, it was reported that the M. cajuputi forests in Terengganu (Malaysia) are degrading due to anthropogenic activities and land cover changes [13].

Predicting biodiversity distribution has become a crucial part of environmental conservation and restoration planning in recent years, and a wide range of modelling techniques were developed for this purpose. In order to understand the abiotic relationship of the M. cajuputi with its geographical distribution and climate impact, this study aims to simulate the current and future suitable distribution of M. cajuputi under different climate change scenarios using the Maxent model.

2. Materials and Methods

2.1. Study Area

The state of Terengganu (3.8779–5.8528° N, 102.3720–103.4987° E) is in the eastern part of Peninsular Malaysia (Figure 1). Terengganu has a tropical climate, and according to the Köppen–Geiger classification system, this climate is classified as Af [14]. The M. cajuputi forest in Terengganu can be found on the ridge and depression (swales) types of soil [5]. In Peninsular Malaysia, Terengganu supports the widest spread of the M. cajuputi species, covering 64.5% (147,489 hectares) of the total M. cajuputi forests [15]. The M. cajuputi forest in Terengganu is found native on soil type known as Beach Ridges Interspersed with Swales (BRIS) or locally called “tanah beris” [16].

2.2. Data Compilation

2.2.1. Melaleuca cajuputi Data

We have conducted a multi-site survey where the distribution of M. cajuputi species was recorded from 11 existing locations during the flowering season (October 2018) and dry season (April 2019). We identified the native location of M. cajuputi species based on previous studies information on the location. The unsupervised classification map complements the existing information on the location of M. cajuputi species to identify the survey locations. The species occurrences were recorded randomly in each cluster produced by the classification map. We have identified a total of 183 occurrence records from the site survey. Only three historical occurrence records of M. cajuputi are known from Marang [17], and the locations of these three records are too close to our occurrence records to be treated as independent locations. Therefore, we decided not to include the historical records and use only the current presence data.

2.2.2. Sampling Bias

Before training the model, the sampling bias was reduced by removing duplicate presence points in the same cell or grid, leaving only one occurrence points per cell. Even though the systematic sampling approach is generally more efficient in reducing the sampling bias in our case, the occurrence records are not abundant, and the bias grid method could outperform the systematic sampling approach [18]. Therefore, a bias file was created using Gaussian kernel density of sampling localities method [19] in the SDM toolbox add-on ArcMap Software [20] and the sampling bias distance was set to 3 km [21]. Then the bias file is input to the model setting to select the background points to reduce the sampling bias [22]. After the thinning process, only 44 occurrence points were retained out of 183 occurrences.

Subsequently, the bias correction was performed using the target background sites method [19]. The randomly located background sites was commonly used in training SDM, which in an unbiased manner, the species presence had an equal probability of being sampled. Nevertheless, in the case of correcting the presence bias, the same kind of bias in background sites can be included to cancel one another out. So, the maximum 10,000 points of pseudo-absence background were set, and we ensure that the points sampled have no missing data from the environmental layers. The pseudo-absence background samples were selected using the surface range envelope (SRE) approach [23].

2.2.3. Environmental Data

There are 27 environmental variables which could relate to M. cajuputi species, consist of climate, topography and soil properties that are available as spatial data within the study duration (see

Table A1. The selected environment variables for modelling the habitat suitability distribution of M. cajuputi.

Data Source	Category	Variables	Abbreviations	Units
SRTM	Topographic	Elevation	elevation	m
WorldClim	Bioclimatic	Annual mean temperature	Bio1	°C
		Mean diurnal range (mean of monthly max. temp.—in. temp)	Bio2	°C
WorldClim	Bioclimatic	Isothermality	Bio3	°C
		Temperature Seasonality	Bio4	°C
		Max temperature of the warmest month	Bio5	°C
		Min temperature of coldest month	Bio6	°C
		Temperature annual range (Bio5–Bio6)	Bio7	°C
		Mean temperature of wettest quarter (3 months)	Bio8	°C
		Mean temperature of driest quarter (3 months)	Bio9	°C
		Mean temperature of warmest quarter	Bio10	°C
		Mean temperature of coldest quarter	Bio11	°C
		Annual precipitation	Bio12	mm
		Precipitation of wettest month	Bio13	mm
		Precipitation of driest month	Bio14	mm
		Precipitation of seasonality (coefficient of variation)	Bio15	mm
		Precipitation of wettest quarter	Bio16	mm
		Precipitation of driest quarter	Bio17	mm
		Precipitation of warmest quarter	Bio18	mm
		Precipitation of coldest quarter	Bio19	mm
Department Agriculture Malaysia (DOA)	Soil properties	Lithology categories	lithology	-
		Topsoil pH	topsoilph	-
		Soil class	soilclass	-
		Topsoil Silt fraction	siltpercent	%
HWSD	Soil properties	Topsoil sand fraction	sandpercent	%
		Topsoil clay fraction	claypercent	%
		Topsoil gravel fraction	gravelpercent	%
Landsat 8 OLI	Soil salinity index	Brightness index	BI	-

for variables description). Even though the M. cajuputi is believed to be resilient to climate change [4], the response to climate change is different geographically, and the species’ adaptability in Malaysia’s climate is unknown. Bioclimatic variables are significant in defining species’ environmental niches. The 19 bioclimatic variables were downloaded from the WorldClim-Global Climate Data, including the present and future climate conditions [24]. Global climate model data is based on CMIP5 (IPCC Fifth Assessment Report, AR5) and the bioclimatic data was downloaded from WorldClim version 2.1 [25] to represent the current climate. For the future climate data, we downloaded the bioclimatic variables from three different general circulation models (GCM): CNRM-CN5, CCSM4 and MIROC5 for future prediction in year 2050. These three GCMs were selected because their total error is less than 25, suggesting a good performance for simulating temperature and rainfall in Southeast Asia [26]. The description of the GCMs total error computation was summarized in Supplementary Materials S3 and the list of performance metrics used is listed in Table S1. All the environmental data used in this study have 30 arc seconds resolution (also referred to as ~1 km spatial resolution).

The digital elevation model (DEM) data were extracted from pre-processed SRTM with 90-m resolution [27] and re-sampled into 1 km resolution by using the nearest neighbor re-sampling technique. It is expected that, during the re-sampling process, some important topographic details may be lost because of aggregation. However, the influence of DEM resolution on the topographic details depends on the surface topographic characteristics [28]. In our study area, topographic characteristics are less variable in elevation even in 90-m resolution. In areas where the species occurrences were recorded, the topography is mostly flat. The soil variables, which consist of eight types of soil properties, were downloaded, and extracted from the Harmonized World Soil Database (HWSD), where all information of the soil properties is available on a global scale [29]. Alternative soil properties data were obtained from the Department of Agriculture Malaysia for the detailed category of soil and lithology of the study area. The data obtained was from 2000 to 2009 and more recent data were not available. Moreover, due to the limitation of soil salinity map over a large area, we extracted several indices, namely salinity index-1 (SI), brightness index (BI) and normalized difference salinity index (NDSI) that can be associated with soil salinity obtained from the Landsat 8 OLI. These indices values were compared and validated with the soil salinity obtained from the field measurement. The soil salinity indices from satellite imagery are commonly used to map the soil salinity in a large area and are applied in many studies [30,31,32]. The soil salinity from the field was measured according to the standard procedure by measuring the electrical conductivity (EC) from the soil sample with the soil-to-water ratio of 1:5 [33]. Of the three salinity indices being examined, BI gives the highest correlation with the measured salinity (EC_1:5) with R² of 0.808. This index was subsequently used to generate the soil salinity index raster layer for the study area.

2.3. Ecological Niche Modelling

The maximum entropy (MaxEnt) approach was employed to model the species distribution using the ENMTools version 1.0.4 [34] package in the R platform with Maxent.jar version 3.4.4. The MaxEnt approach was chosen because of its advantages, which include: (i) it has a sound mathematical sense; (ii) use presence-only data; (iii) the model can handle both continuous and categorical environmental data; (iv) the model can function well even with low sample sizes; (v) simplicity in model interpretation, and (vi) the model provides good accuracy and prediction power [35,36,37,38,39]. All the raster layers input variables were re-sampled and clipped to the same extent and spatial resolution following the Terengganu’s state boundary with approximately 1 km spatial resolution. The preparation of raster layer data for the MaxEnt input and mapping outputs were processed using ArcGIS 10.5 with add-in SDM toolbox [20]. We trained the Maxent model using the sub-sample method, where 30% of the presence points were set aside for model testing evaluation. The model was trained by replicating each model ten times and the final output used was the average of all runs.

2.3.1. Variable Selection

To select variables for the model, first, we analyzed the variable permutation importance from the MaxEnt model run with 27 variables using the default setting. The variable importance was computed using the permutation method by the “enmtools.vip” function, which uses the vip package in R. From the “vip-vignette”, the permutation approach idea is that the training performance would be degraded if the values of an essential predictor in the training data being randomly permuted [40]. So, in the case of variable importance computed by the enmtools.vip function, the results were based on the AUC metric, and the process was repeated ten times.

Then we compared the result of the MaxEnt variables’ permutation importance and contribution ranking with the variance inflation factor (VIF) [41]. The step-wise VIF procedure was performed using the uncertainty analysis for the species distribution models package (USDM) in R [42]. Following the step-wise VIF with a threshold of less than seven, we screened the most important variables by minimizing the multicollinearity between variables. In this regard, the variables with the highest collinearity are excluded. However, to avoid eliminating the essential variables that could be informative to M. cajuputi, variables that show multicollinearity were removed. Therefore, we intercepted the variables with a permutation importance and contribution greater than 1% with the VIF variables selection. The final selected variables went through another step-wise VIF procedure with a threshold of ten. Figure 2 below represents the overview of our modelling process.

2.3.2. MaxEnt Model Complexity and Model Performance Assessment

ENMTools functions allow a user to control the model complexity by adding arguments, and for Maxent, the model complexity depends on two parameters: the model features and regularization of the parameters setting. The model regularization is necessary to control the model overfitting. It is important to use suitable regularization number and feature setting to get the optimal model using Maxent algorithm for different species and research goals [36]. Maxent model produces a model based on a set of “features” such as linear (L), hinge (H), quadratic (Q), threshold (T), and product (P), that help the model determine the ratio of the predicted probability of presence or the environmental suitability of the study area given the site’s particular environmental values [35,43] (details of each feature types function are described in Supplementary Materials, S1). In this study, we trained the model with features combination (FC) setting: “L, H, LQ, LQH and LQHPT”, and regularization multiplier (rm) value setting from 0.5 to 5 with 0.5 increments.

The optimal model selection was based on the model calibration accuracy assessment and discrimination accuracy, that is, threshold-independent metrics. The discrimination accuracy was determined by the area under the receiver operating curve (AUC). An AUC value less than 0.5 indicates that the model prediction is poor and that closer to 1.0 indicates a robust model prediction [35]. The AUC_diff is the difference between the AUC value based on the training data (i.e., AUC_train) and AUC_test (AUC_train–AUC_test). If AUC_train is smaller than AUC_test, the returned value is zero. The overfit models are expected to have high AUC_diff in the notion that the overfitting model performed well on training data but poorly on test data [44]. For many model purposes, it was strongly suggested that the calibration accuracy was more beneficial for model selection than discrimination accuracy [45]. In this study, the calibration accuracy was determined based on the expected calibration error (ECE) and maximum calibration error (MCE). We also used the continuous boyce index (CBI) as an additional assessment tool because it could provide insight into the model’s robustness and deviation from randomness [46] or simply means the model transferability to a different geographical area. More details about the CBI and other performance assessment can be found in Supplementary Materials S1. We compared the modeling results using the default setting against the optimized model. Then, the selected optimal model was projected to the future year 2050 in different climate scenarios. These climate change scenarios are based on the IPCC Fifth Assessment Report (AR5) and the representative concentration pathway (RCP) data defined by the possible range of radiative forcing values (2.6, 4.5, and 8.5 W/m²) in the year 2050. The future bioclimatic variables used were the ensemble of the selected GCMs, and for future projection, the soil properties and topography variables are assumed unchanged. The final output produced the species suitability distribution maps with habitat suitability index ranges from 0 to 1. These index values were then categorized into four levels of suitability using equal interval classification from unsuitable to high.

3. Results

3.1. Variables Selection and Model Performance

Ten variables were selected out of 27 variables after they went through the variable selection process. Bio 2, Bio 8, Bio 9, Bio 13, Bio 14, Bio 16, elevation, silt percent, lithology and BI are the variables that satisfy the selection criteria. The highest collinearity is between variables Bio16 and Bio13 with 87% positive correlation (Figure 3). Other variables have collinearity less than 65%. Figure 4 shows the variables’ importance in the final model after they went through the variable selection process. Based on the preliminary model variables’ importance and VIF step-wise procedure, variables Bio16 and Bio14 were not selected because of high multicollinearity. However, as the variable importance and permutation importance (List 1) are higher than the selected VIF variables (List 2), we included variables Bio16 and Bio14, and excluded the one with no permutation importance in the model.

The model performance with a different variation of FC and regularization multiplier can be found in Supplementary Materials, Figure S1. The results suggest that the model’s discrimination accuracy of different features type is equally good under the same number of regularization multipliers (rm). Linear features have the most significant variation from all the different features combined under different regularization multiplier values. Furthermore, the modelling results for different features type under rm = 0.5 have the lowest ECE, MCE and CBI values. However, for a good model performance, the CBI value should be at the higher side.

Overall, the model performance with different features and regularization settings in terms of discrimination accuracy were all excellent. The AUC_test scored above 0.8 and the CBI is greater than 0.85 except for the models under regularization value of 0.5.

Table 1. Performance assessment comparison of Maxent models with default setting and the optimized model parameter values, namely ML0.5 and ML1. ML1_a is the ML1 model before bias correction and ML1_b is the chosen optimized model with sampling bias correction.

Model	Default	ML0.5	ML1_b	ML1_a
Feature	Default	L	L	L
Rm	1	0.5	1	1
AUC_test	0.858	0.853	0.858	0.845
AUC_diff	−0.003	0.030	−0.003	0.048
CBI	0.931	0.838	0.931	0.911
ECE	0.347	0.279	0.313	0.262
MCE	0.066	0.056	0.059	0.058

compares the default model setting, the selected optimized model with the smallest ECE (ML0.5), the optimal model based on the overall ECE, MCE, CBI and AUC_diff score (ML1_b) and model ML1_b without employing bias file (ML1_a). The best model to be projected under future climate is selected based on the high CBI with low ECE, MCE and AUC_diff scores, which is model ML1_b with linear feature setting and regularization multiplier of one.

Analysis of the model’s variable importance suggests that ML1_b highly depends on variable Bio 2 (mean diurnal range) with 34.28% variable importance values (Figure 4). The following two essential variables are elevation and the mean temperature of the wettest quarter (Bio 8), with variable importance values of 19.22% and 13.12%, respectively. As mentioned at the beginning of this section, Bio14 and Bio16 should be excluded if the analysis is confined to the VIF selection, but the variables contribution is significant. The feature importance relies on the model error estimates, and also the collinearity between the predictor variables could affect the variable permutation importance result [47]. This statement explained why some highly correlated variables could interfere with the model performance and can be excluded. For example, Bio9 (mean temperature of the driest quarter) variable has the lowest percentage (0.94%) in the model ML1_b despite showing negatively high correlations with Bio16, Bio13 and Bio8 (see Figure 3).

3.2. Species Response and Potential Habitat Suitability Distribution

Species’ responses to changes of the selected variables are depicted in Figure 5. The response curves show how each environmental variable varies with the species locations. The curves indicate how species habitat suitability is shifted against the variable of interest (x-axis) when all other environmental variables are held constant at their average sample value.

Based on the response curves, the species habitat suitability remains high (>0.85) for the mean diurnal temperature range (Bio2) less than 7.4 °C, and from this point, the habitat suitability decreases. The species’ response toward Bio2 suggests that the most suitable habitat is achieved when the variation between the mean monthly maximum and minimum temperature is small. It is interesting to note that the habitat suitability shows a rapid increase against the mean temperature during the wettest quarter (Bio8). During the driest quarter (Bio9), the index remains high at about 0.8 for temperature between 26.4 °C and 27.0 °C during the wettest and driest quarters.

As for precipitation, the variable Bio13 (precipitation of the wettest month) significantly contributes to the model with a variable importance score of 9.3% (Figure 5). The highest habitat suitability was obtained when the precipitation during the wettest month ranges from 575 mm to 600 mm. Beyond this range (575–600 mm), the habitat suitability decreased significantly. However, for the precipitation of the wettest quarter (Bio16), the habitat suitability value remains high when the precipitation ranges from 1200 to 1450 mm, and the highest presence density is expected when the precipitation exceeds 1350 mm. Moreover, the habitat suitability value is high when the precipitation during the driest month (Bio14) is above 55 mm. Surprisingly, the highest presence density occurred when the precipitation during the driest month is less than 60 mm.

Moreover, it was reported that M. cajuputi is sensitive to salinity [4], and from the response curve of BI variable the species is mainly found in areas with salinity index from 2.0 to 2.1. However, the habitat suitability remains high in less saline soil with a salinity index from 1.6 to 2.6. The response curve of the silt percentage in the topsoil (siltpercent) remains at high habitat suitability for silt content between 10–30%. However, the species presence is most dominant in areas with a silt percentage below 15%. Furthermore, 3 out of 16 lithology classes have contributed significantly to the optimal habitat suitability, where the species presence is mainly distributed on: (i) clay and silt (marine), (ii) sand (mainly marine) and (iii) peat, humic clay and silt. Interestingly, elevation has the most significant influence, especially in low altitudes ranging from 10 to 30 m above mean sea level, and the presence density peak is at an altitude approximately 10 m above sea level. The response curves of the soil properties also fit well with the description of the habitat on the coastal and BRIS type of soil in Terengganu [13].

3.3. Changes of Habitat Suitability under the Present and Future Climate Scenarios

The predicted distributions of M. cajuputi under the current and future climate conditions in Terengganu are shown in Figure 6. The present-day habitat suitability index (HSI 2019) identifies stretches of area suitable for M. cajuputi habitat along the coastal area from Besut to Kemaman District (refer to Figure 1 for District’s location and boundary). Based on the 2019 baseline map, only 10.82% of the total area of Terengganu falls under suitable habitat. Of which, 7.42% fall under poor, 2.44% medium, and 0.96% high habitat suitability (Figure 6e). The habitat suitability areas are expected to undergo significant changes for different periods under different climate change scenarios. Compared to the current suitability prediction, it is surprising to see that the level of habitat suitability increased in the future under RCP 8.5 extreme climate scenarios for specific regions (Figure 6). However, the predicted future habitat suitability of M. cajuputi is likely to be degraded, from optimal suitability to poor in some regions under different climate scenarios, as demonstrated in Figure 6 and more details on the area suitability changes in Figure 7.

Based on the RCP 2.6 climate scenario, the projection in year 2050 reveals that the habitat suitability for M. cajuputi would decrease further in Marang, Dungun and Kemaman districts where the habitat will change from suitable to unsuitable (Figure 7a). However, the level of habitat suitability shows large increases in Kuala Terengganu, Setiu and Besut districts, as marked in light blue in Figure 7a. The largest contraction area of habitat suitability of about 3.76% is under RCP 4.5, distributed from the North to the South of Terengganu. Moreover, the habitat degradation from high/medium suitability to poor was predicted the most under climate scenario RCP 4.5 with a total reduction of 1.09%. On the contrary, a remarkable expansion of suitable habitat area is expected under the climate scenario RCP 2.6 with a total expansion of 3.62%, especially in Kuala Terengganu, Marang and Dungun districts (see the dark-blue marked area in Figure 7a). Under the extreme climate scenario (RCP 8.5), the area of suitable habitat is predicted to decrease in some parts of Setiu, Kuala Terengganu and Kemaman districts. Furthermore, the level of high suitability habitat is projected to increase in Setiu, Kuala Terengganu, and Dungun districts (Figure 6a,d). This can be seen clearly in Figure 7c, where the area that is predicted to improve from poor suitability to moderate or high is shown on the map in light blue (152.46 km²).

From all the future projections under different climate scenarios, the habitat suitability of M. cajuputi species in Setiu district is the most stable where the area of contraction is minimal, and the habitat suitability remains suitable even though the overall level of habitat suitability decreased under different climate scenarios. As for the Marang area, the habitat suitability seems unstable under different climate scenarios, as shown in Figure 7. The changes in habitat suitability level show no distinct pattern in areas previously classified as highly suitable habitat (notice the “Unchanged” areas in Marang). The suitable habitat in Dungun and Kemaman (southern area) decreased significantly in 2050 under all climate scenarios. Figure 6e highlights the changes of the unsuitable area under the current climate (2019) and in 2050 under different RCPs where the area of unsuitable habitat is predicted to increase in the future (noted the percentages in the stacked bar).

Overall, the modelling results seem to be strongly influenced by the soil properties, especially the elevation and soil salinity and bioclimatic, which can be referred to the model variable importance (Figure 4). Moreover, the variable importance indicates that the potential species habitat suitability is also controlled by the temperature and precipitation during the wettest quarter, as evidenced by the response curve where the habitat suitability is high during times of high temperature and high precipitation. This could be associated with rainfall seasons during the northeast monsoon (November to March) in Terengganu, which also coincides with the coldest months (December, January and February) [48]. The distribution of M. cajuputi forest obtained in this study is comparable with those reported by the World Wide Fund (WWF) using Landsat 5 in 2010 [12], but the map classification was not verified.

Data from the World Weather Online (WWO) website was used to derive statistics on average precipitation and temperature, which were reliably sourced from global weather satellites, the World Meteorological Organization, and a global telecommunication system [49]. The climatic patterns from 2009 to 2020 for Setiu, Kuala Terengganu (K. Trg) and Marang districts are shown in Figure S2 (Supplementary Materials S2). Generally, the climate trend obtained from the WWO agrees with the earlier study where the climate changes were analysed during the coldest months (December–January–February) that coincide with the northeast monsoon and the warmest months (June–July–August) that coincide with the southwest monsoon [48].

4. Discussion

4.1. Modelling Techniques and Performance Assessment

Many studies have used various species distribution techniques to understand the species’ ecological niche distribution and predict the potential habitat suitability. Among the species distribution models, MaxEnt model seems to be the preferred one for SDM studies [37,50,51,52,53]. The ease of use and simple steps necessary to run MaxEnt appear to have enticed many researchers to use it as a black box, despite mounting evidence that using MaxEnt with default parameter values (i.e., auto-features) does not always result in the optimal model [54,55]. Nonetheless, like any modelling approach, the quality of input data (i.e., the reliability of environmental and species presence data) and the parameterization employed for modelling significantly impact the model’s results using MaxEnt [56,57]. Issues of multicollinearity among the predictor variables that could affect the model performance were highlighted in the previous literatures [58] and reducing the number of variables could reduce the model complexity and benefits the operation time and model interpretability. In this study, we carefully selected the variables by comparing the list of variables favored by step-wise VIF and variable permutation importance. This variable selection method can avoid eliminating the importance variable for the species model. Furthermore, the final selected variables have VIF scores less than ten. Even though variables Bio16 and Bio13 have the highest collinearity, as shown in Figure 3, the permutation importance score of these two variables still differs by 2.93%. This is because Maxent can regulate redundant variable contributions, making it resistant to predictor collinearity in model training [59]. However, to avoid eliminating essential variables by simply excluding the variables with high multicollinearity, we suggest that both the variable permutation importance and VIF scores are considered as the basis for variable selection, besides the variable relation to the species.

When comparing the model complexity, the most used metrics assessments are information criterion such as Akaike information criterion (AIC) or AICc for smaller samples. However, the use of AIC as an estimate of model parsimony may lead to uncertain or confusing performance of the MaxEnt model [58]. Using AIC is acceptable for application with another model such as GLM or GAM but in the case of the MaxEnt model, the use of AIC or AICc is debatable because of a philosophical disconnection between machine learning and classical modeling [44]. We found that using the model calibration assessment, such as ECE and MCE, could help distinguish the model performance when the model performed almost the same for AUC_test and CBI under variation of FC parameters. However, the use of CBI helps reflect the propensity of the model prediction in the study area and determines the optimal model selection. Moreover, the evaluation of the threshold-independent metrics is preferable for the MaxEnt model because it uses only presence data [60]. Because of the residuals problem even after the bias correction, predictions made with clustered ecological data were not as good as compared to models using herbarium datasets [21]. In our case, by employing the bias file in the model, it improved our model prediction accuracy and correlation values. Although the value differences in AUC is very small, for CBI assessment the bias correction significantly improves the CBI scores.

The common practice for modelling potential species distribution is using the bioclimatic variables as predictors. However, in predicting plant species distribution or ecological niche modelling, the climate-only model is incomplete as most plant species are highly influenced by the soil properties [61]. Our results show that both bioclimate and edaphic predictors are equally important, which is highlighted by the model’s variable importance. It was found that combining both data as predictors has led to a better modelling result as both climate and edaphic variables could help minimize an overestimation of the species distribution range when predicted to the future [62,63]. In addition, the species that are more influenced by the soil property than climate would be less prone to climate change effect, and if the bioclimate has more influence on the species, their soil properties will make the species distribution in a more acceptable range.

4.2. Melaleuca cajuputi Potential Suitable Habitat and Climate Change Impact

The RCPs represent four alternative greenhouse gas (GHG) emissions and atmospheric concentrations, air pollutant emissions, and land use scenarios for the twenty-first century. The Special Report on Emission Scenarios (SRES) is also commonly referred to in climate change assessment studies, but RCPs which represent scenarios with climate policy cover a wider range of scenarios compared to the SRES. In terms of overall forcing, the RCP 8.5 is broadly comparable to the SRES A2/A1F1 scenario, RCP 4.5 to B1 scenario and for RCP 2.6, but there is no equivalent scenario in SRES. The RCP 2.6 represents a scenario in which global warming is limited to less than 2 °C above pre-industrial temperatures. While RCP 4.5 represents the intermediate scenario and RCP 8.5 represents a very high GHG emission scenario [64]. The most frequently used scenario in climate projection studies in Malaysia is A1B (balance across all sources of technological change in the energy system) in which there is no similar scenario for RCP.

As for M. cajuputi species, the future projection interestingly affects the species distribution differently where under the most extreme weather, the habitat suitability distribution significantly decreased especially in Setiu and Kuala Terengganu districts, and the level of suitability seems increasing toward the north, from the Dungun to Marang districts. The local climate in Setiu has significantly greater monthly average high temperature and monthly average precipitation trends compared to K. Trg and Marang during non-monsoon months (details in Supplementary Materials S2). The significant variation in temperature and precipitation for the three regions could explain why the projected future habitat suitability of M. cajuputi in Setiu remains high (Figure 6 and Figure 7). As mentioned in Section 3.2, the M. cajuputi species is sensitive to salinity, and the salinity variable ranked the highest four by contributing 9.88% in the model prediction. According to Nguyen et al. (2009) [65], a young M. cajuputi species can survive up to three months in a pot with soil salinity of 5 dS/m, while other Melaleuca species can survive up to 10 weeks with soil salinity from 2 to 60 dS/m [4]. Our results highlighted that the species habitat suitability in Terengganu are high, mainly on soil with a salinity range between 1.7 and 2.2 dS/m. Moreover, our finding in Section 3.3 suggests that the species presence localities are primarily in areas with two climate conditions: high precipitation during the wettest months and low precipitation during the driest months. This indeed agrees with the finding in Thailand that M. cajuputi species can adapt well to flooding and dry conditions and grow taller in water-logging areas than in dry areas [9]. In a climate study, the future sea level on the coast of Peninsular Malaysia was simulated, and it was demonstrated that the mean sea level would continue to increase by 0.517 m between 2040 and 2100 especially during the monsoon months [66]. Furthermore, the Malaysian Meteorological Department (2009), using GCM with average temperature and rainfall from 1961 to1999 as a baseline, projected a temperature rise for Peninsular Malaysia between 1.1 °C and 3.6 °C by 2095 [67,68]. Another study projected increases of 0.32°C per decade in the mean daily temperature for Peninsular Malaysia amounted to a total increase of 0.96°C in three decades (2019 to 2050) [69]. Because of the expected sea level rise in the future [67], saltwater intrusion will be more common and will increase the soil salinity level. Significant groups of terrestrial and freshwater species are unable to adapt fast enough to stay within the spatially changing climatic envelopes at high rates of warming.

The fast-changing climate is worrying, especially to the coastal ecosystem. M. cajuputi species was reported to have high adaptability toward climate and other environmental conditions. However, this species is not only threatened by climate change but also other factors such as land reclamation, land development, agriculture, and illegal waste dumping [70]. If no effort of conserving or rehabilitating the M. cajuputi forest, the coastal ecosystem of Terengganu is at risk of disappearing. As mentioned in Section 2.1, M. cajuputi is found native on BRIS soil and BRIS soil vegetation can easily catch fire, particularly in non-monsoon months or drought season, which can occur naturally or induced by human activities [70]. However, since M. cajuputi is a fire-adapted species which can increase its survival rates during a fire or regenerate naturally after a forest fire, the ecosystem resilience is less affected. Our modelling has simulated that by 2050, the habitat suitability area in Setiu district is predicted to remain stable with some patches of area decrease and increase in suitability level, but for Kuala Terengganu and Marang districts, the suitability area varies under different climate scenarios. Unfortunately, in the Dungun and Kemaman districts, the suitable habitat area would disappear if the current climate change trend continued as predicted under scenarios RCP 4.5 and RCP 8.5. Based on the model projection, it is worrying that the current reserve of M. cajuputi forest in Marang will continue to deplete due to climate change while in Setiu area, the threat is associated with anthropogenic activities.

The habitat suitability distribution was simulated by the variable’s response to the presence samples regardless of bias of the land use and land cover (LULC). This is the limitation of our model, which does not consider the current land use as a variable in the modelling. Consequently, the model predicts habitat site suitability based on the bioclimatic and edaphic data only. However, by knowing the potential suitable habitat of M. cajuputi in the current and future states, it helps to make a better management decision or policy for coastal ecosystem. Maintaining a healthy M. cajuputi forest will aid in the preservation of ecological services for the benefit of the coastal environment, which supports the livelihood of local communities.

5. Conclusions

In this study, we successfully modeled the habitat suitability of M. cajuputi under the current and future climate change scenarios. Our prediction showed that under the current climate, Setiu region shows a better M. cajuputi habitat suitability compared to Marang region in the long term. The predicted habitat suitability area is well distributed along the coast and highly suitable on BRIS type of soil with low percentages of silt and a high percentage of sand. The species distribution could also be affected by the temperature and precipitation during monsoon months. Over the long term, soil salinity could also affect the survival of M. cajuputi, especially when the sea level rise becomes more apparent. In the future, the modelling result could be improved using the local bioclimate data and land use data to get a better prediction of the habitat distribution under the current climate and LULC. In support to the government’s effort to promote ecological conservation, we recommend the inclusion of future climate scenarios in the existing conservation and restoration strategies. The niche model of M. cajuputi species in Terengganu could give significant insight to conservationist in their strategy to strengthen sustainable practices especially in the Setiu district, which is currently surrounded by oil palm plantations.