Maxent modelling for habitat suitability of vulnerable tree Dalbergia latifolia in Nepal

Dalbergia latifolia Roxb., commonly known as rosewood, is one of the highly valuable tropical timber species of Nepal. The tree species was widely distributed in the past, however, over-exploitation of natural habitat, deforestation, forest conversion for agriculture, illegal logging and the invasion of alien species resulted in the classification of this species as vulnerable by the IUCN (International Union for Conservation of Nature) category. So, the prediction of habitat suitability and potential distribution of the species is required to develop restoration mechanisms and conservation interventions. In this study, we modelled the suitable habitat of D. latifolia over the entire possible range of Nepal using a Maxent model. We compiled 23 environmental variables (19 bioclimatic, 3 topographic and a vegetative layer), however, only 12 least correlated variables along with 43 spatially representative presence locations were retained for model prediction. We used a receiver operating characteristic (ROC) curve to assess the model’s performance and a Jackknife procedure to evaluate the relative importance of predictor variables. The model was statistically significant with an area under the curve (AUC) value of 0.969. The internal Jackknife test indicated that elevation was the most important variable for the model prediction with 71.3% contribution followed by mean temperature of driest quarter (9.8%). The most (>0.6) suitable habitat for the D. latifolia was 235 484 hectares with large sections of area in two provinces whereas, the western most provinces were not suitable for D. latifolia as per Maxent model. The information presented here can provide a framework for nature conservation planning, monitoring and habitat management of this rare and endangered species.


Introduction
Geographic position, altitudinal variation and climatic conditions create different ecosystems and habitats that dictate the distribution of vegetation. In addition, soil fertility, topography and irradiance play vital role in habitat distribute in local/small scale. Habitat which is combination of the space and all abiotic environmental factors and other organisms (Yi et al. 2016), plays a significant role in the distribution and existence of species in different location. Forest are the most vital and valuable renewable natural resource of the globe, and acts as a prime repository for a significant portion of terrestrial biological diversity (Bozzano et al. 2014). Furthermore, forests play a dynamic role in the protection of the fragile biodiversity of mountain ecosystems, and they support the livelihoods of forest-based communities (Thomson 2000). In this century, due to the combined impacts of urban expansion, global warming, forest fires, over-exploitation, and the invasion of alien species, our forest resources are in imminent danger (Seto et al. 2012), resulting in twenty percent of the plant species in the risk of extinction (Nic Lughadha et al. 2020) globally. Recently, in-situ and ex-situ conservation approaches have been adopted by conservationists to conserve rare and endangered wild flora and fauna under challenging conditions. So, from this conservation perspective it is crucial to monitor and assess potential suitable habitat of threatened species in order to rehabilitate them in their ecosystem (Gaston 1996).
Different environmental characteristics have significant roles in limiting the distribution of species (Wiens 2011). So, the relationships between species and the environment are the most used parameters for species distribution modelling and mapping under prevailing conditions. Models for predicting the potential geographic distributions can be used for habitat suitability mapping and conservation plan preparation of globally important species (Graham et al. 2004;Guisan and Thuiller 2005). Numerous species distribution modeling methods are available to simulate the spatial distribution of plant species along with the study of spatial patterns of species diversity (Wisz et al. 2008;Graham et al. 2006) and impacts of climate change (Saran et al. 2010). However, to assess the habitat suitability for rare and endangered plant species, very few predictive models are developed because of comparatively small sample size (Pearson et al. 2007). The most frequently used species distribution models are Maxent, random forest, boosted regression trees, generalized additive models, and multivariate adaptive regression spines, which have similar predictive performances (García-Callejas and Araújo 2016). In this study, we used the maximum entropy algorithm (Maxent) because it requires only species presence (or occurrence) data along with environmental information, and it deals effectively with limited occurrence data of small sample sizes (Merow et al. 2013;Elith et al. 2011). It is based on the maximum-entropy principle, where the best approximation to an unknown distribution is that which maximizes entropy while satisfying any known constraints (Jaynes 1957). Because of its high performance with both spatially biased data (Loiselle et al. 2008) and limited occurrence records (Pearson et al. 2007), it has been one of the most frequently used  modelling techniques (Ortega-Huerta and Townsend Peterson 2008). In addition, this model can use both continuous and categorical environmental variables as input variables and incorporates the interactions between the variables ) and best suited for countries like Nepal with diverse geographical and environmental conditions. Despite this, Maxent has some limitations such as sampling bias of occurrences, selection of features and the region used for background sampling (Elith et al. 2011;Kramer-Schadt et al. 2013).
Dalbergia latifolia Roxb. (Fabaceae) is a slow growing, predominantly single-stemmed deciduous tree that reaches over 40 m in height and found commonly within Dalbergia sissoo Roxb. (Jackson 1994) and other deciduous forests (Nath et al. 2011). This species occurs naturally up to 1000 m altitude (Jackson 1994) in India, Nepal, Indonesia (Soerianegara and Lemmens 1994) and Malaysia (Praciak et al. 2013) , however large scale plantations have been established across Kenya, Myanmar, Nigeria, Philippines and Vietnam (Orwa et al. 2009). D. latifolia is one of the most valuable timber trees in the Indian sub-continent (Praciak et al. 2013), where bark and leaves are locally used as medicine (Jain et al. 2005;Padal et al. 2010). Moreover, bark is also used as ethno-veterinary medicine in India (Selvaraju et al. 2011). Due to its high value for multiple uses (Troup 1921), it is over-exploited in its natural range (Thapa 2004), which has resulted in a ban on the harvest, transportation and export for commercial purposes of this species (GoN 2001) in Nepal. D. latifolia is protected under the Indian forest act, 1927 and also in Forest act, 2019 of Nepal and reported as "Vulnerable" by International Union for Conservation of Nature (Damaiyani and Prabowo 2019). Despite the enormous threats, there have been very limited concentrated efforts to address conservation concerns, such as mapping of the distribution pattern in Nepal. Over the last few decades, the majority of the studies were focused on the genetics and ex-situ conservation of D. latifolia, however, the availability of suitable habitat for large scale plantation has not been determined. Information on the potential habitat availability of this vulnerable tree species can be important conservation tool, in order to permit conservation initiatives, reintroduction and restoration in its preferred habitat.
In this study, we used a Maxent model to predict the potential distribution of D. latifolia across Nepal and to identify the environmental factors associated with this distribution. Due to varying topographic condition and the formation of micro-climatic zones in small altitudinal gradients, environmental factors alter more frequently from east to west than north to south of the country. For this, we hypothesized that D. latifolia has a distribution range from eastern to mid-western parts of Nepal up to 1000 m altitude and that bioclimatic and topographical factors, and companion vegetation are the major influencing variables determining its distribution.

Study area
The study was carried out in Nepal, which has an area of 0.147 million square kms and lies between India and China in the latitudes of 26°22´N to 30°27´N and longitudes of 80°04´E to 88°12´E (Fig. 1). The variation in physiography associated with the range in elevation from 60 m to 8848 m (Mount Everest) resulted in a variety of biomes (McColl 2014) in five distinct physiographic regions: High Himalaya, High Mountain, Middle Mountain, Siwalik, and Terai (Uddin et al. 2015). The average annual rainfall is 1600 mm, but it varies by eco-climatic zone, which form different micro-climatic conditions within broader physiographic zones. According to land cover mapping, both forest and other wooded lands together cover 6.61 million hectares, 44.74% of total area of the country with altogether 443 tree species belonging to 239 genera and 99 families (DFRS 2015). Even though Nepal covers a relatively small area, it harbors 35 forest types (Stanton 1972), which are often categorized into ten major groups (MoFSC 2014) and 118 ecosystems (Dobremez 1970), where numerous endemic flora and fauna species of global importance are recorded.

Presence Record and Spatial Filtering
The information regarding the occurrence of D. latifolia was obtained from meetings with Division Forest Offices (DFOs) of Terai and Inner-Terai regions of Nepal. After first level of screening, focus group discussions were organized with related stakeholders and local people in order to get the occurrence information. Field surveys were conducted to determine the occurrence locations and to record the phenotypic information in different locations representing all climatic zones, physiographic regions and management regimes. Altogether, 230 presence locations were recorded mostly from the tropical mixed hardwood forest and sub-tropical region of the country with higher concentration in the eastern part (Fig. 1). These records were from 11 different districts (Bardiya, Kapilvastu, Palpa, Nawalparasi, Makawanpur, Bara, Parsa, Rautahat, Sarlahi, Sunsari and Jhapa).
Presence records collected from the field may pose sampling bias due to the variety of disturbances from the connected road networks, urbanization and by focusing the particular vegetation type/ forest stand (Kadmon et al. 2004). To minimize the sampling bias the data need to be collected using spatial extent of the region with planned sampling regime (Barry et al. 2006). However, having zero spatial error would be only possible in the ideal data set. And in order to minimize this bias, we managed spatial filters of grid sizes 2 km × 2 km by retaining only one presence data in each grid cell. Thus, only 43 presence data were used for model prediction in Maxent.

Environmental Variables
We used different environmental layers of climatic, topographic, and non-climatic data as potential predictors of D. latifolia habitat distribution. These were primarily chosen on the basis of their biological relevance to species distribution and other habitat modelling studies (Kumar et al. 2006;Pearson et al. 2007;Xu et al. 2019). For climatic layers, we used 19 bioclimatic variables (Table 1) developed by Hijmans et al. (2005), produced as interpolated climatic surfaces for global lands other than Antarctica using 1950-2000 climatic data (Fig. 2). Moreover, the topographic factors, such as slope, aspect and altitude, also play significant role for the growth of a species (Karami et al. 2015). Digital Elevation Model (DEM) data was obtained from worldclim website (https:// www.worldclim.org/) and it was used to get slope and aspect using R 4.0.0.
To minimize correlations between the predictor variables, multicollinearity among all of the variables was examined by cross-correlations (Pearson correlation coefficient, r). Among different highly cross-correlated variables (r > ±0.75, alpha = 0.05), only one was included based on the potential biological relevance to the distribution and for the ease of interpretation (Kumar et al. 2006). Thus, predictor variables were reduced to 12 from 23 for mapping of potential suitable habitat.  (Table 1). These spatial products were derived from 19 original (modification was not made) bioclimatic variables developed by Hijmans et al. (2005) (posted on the Bioclim website; https://www.worldclim.org/). Scaling factor (see Table 1) is used in the units of the studied variables. (Note: "°C" is degree centigrade, "mm" is millimeter and "%" is percentage).  (Table 1). These spatial products were derived from 19 original (modification was not made) bioclimatic variables developed by Hijmans et al. (2005) (posted on the Bioclim website; https://www. worldclim.org/). Scaling factor (see Table 1) is used in the units of the studied variables. (Note: "°C" is degree centigrade, "mm" is millimeter and "%" is percentage).

Mapping and modelling
We used Maxent software, version 3.4.1 (https://biodiversityinformatics.amnh.org/open_source/ maxent/) to predict the potential suitable habitat of D. latifolia in Nepal. Maxent is a maximum entropy-based machine learning program that takes presence only data as an input, as well as a set of environmental covariates of a defined geography to provide a probability distribution model as a final output (Merow et al. 2013). After completion of the distribution modelling and mapping, we assessed the statistical significance of the model. In the majority of modelling work, independent datasets are generated for the purpose of model validation by splitting the field data into training and test/validation (Fielding and Bell 1997;Guisan and Zimmermann 2000). However, for the modelling of endangered and threatened species, splitting of data sets for validation purposes may not be possible because of the very limited number of samples. So, in this study, we used a Jackknife procedure to examine the importance of individual variables for model prediction as proposed by Pearson et al. (2007). Jackknife predicts the importance of the environmental covariates under the basis of area under curve (AUC) gains for three different scenarios (without variable, with only one variable and with all variables). Beside this, we used the receiver operating characteristic (ROC), and the AUC method to evaluate the model's goodness of fit. The model with the highest AUC value was considered as the best performer. AUC measures the predictive performance of the model by  (Table 1). These spatial products were derived from 19 original (modification was not made) bioclimatic variables developed by Hijmans et al. (2005) (posted on the Bioclim website; https://www. worldclim.org/). Scaling factor (see Table 1) is used in the units of the studied variables. (Note: "°C" is degree centigrade, "mm" is millimeter and "%" is percentage).
comparing the ability of the model to random prediction, and values range from 0 to 1 ). The final suitability map had values ranging from 0-1, which was further regrouped into four classes for the purpose of interpretation: 0-0.2 is unsuitable; 0.2-0.4 is low; 0.4-0.6 is moderate; and 0.6-1 is high suitability (Ansari and Ghoddousi 2018).

Model validation and key environmental variables
The Maxent model predicted the potential suitable habitat for D. latifolia with satisfactory statistical accuracy, with an AUC value of 0.969. The Maxent model's internal Jackknife test indicated that the variables with higher importance in predicting habitat suitability were 'Elevation' and 'Mean temperature of driest quarter' (Fig. 3). These two environmental variables presented the highest gain in comparison to the rest of the others and hence contained the most valuable information. Similarly, the results for the response of the model to each variable showed consistent results with the internal Jackknife test. The most highly contributing variables were Elevation (71.3%); Mean temperature of driest quarter (9.8%) and Slope (6.4%), whereas, other remaining variables contributed less than 6% to the final model (Table 2). Fig. 3. The Jackknife test for evaluating the relative importance of environmental variables for Dalbergia latifolia Maxent model on habitat suitability. (Note: "ASP" is aspect; "Bio14 is Precipitation of driest period; "Bio15" is Precipitation seasonality; "Bio17" is Precipitation of driest quarter; "Bio18" is Precipitation of warmest quarter; "Bio19" is Precipitation of coldest quarter; "Bio2" is Mean diurnal range; "Bio3" is isothermality; "Bio9" is Mean temperature of driest quarter; "ELE" is elevation; "FT" is forest types and "SLO" is slope).  Fig. 4A indicates that the distribution of D. latifolia with a > 0.5 probability of presence was limited to between approximately 60-300-meters in altitude. Outside of this window the distribution of the species sharply declined with increasing altitude, with very low rates of presence above 1000 meter. The probability of presence increased sharply with the rise of mean temperature of Mean diurnal range/ mean of monthly max. and min. temp. (°C) 0.8 100 Bio14

Relationship between major environmental variables and geographical distribution
Precipitation of driest period (mm) 0 100 Bio17 Precipitation of driest quarter (mm) 0 100 Bio3 Isothermality 0 100 the driest quarter and the highest probability of presence was obtained within the temperature range of 22-25 °C (Fig. 4B). Similarly, the probability of presence increased with rise in monthly precipitation from 100 to 125 mm and thereafter sharply decreased (Fig. 4C). Whereas, a negative relationship was found between precipitation of coldest quarter and probability of presence of D. latifolia (Fig. 4D).

Potential distribution of Dalbergia latifolia in Nepal
The Maxent model predicted that Central as well as Eastern parts of the Terai region of Nepal were the most suitable habitats for D. latifolia with very low habitat suitability in the Western region.
The model result showed that an area of 222 403 hectares had high potential, a 111 712-hectare area showed moderate potential while a 269 583-hectare area showed low potential for the growth of D. latifolia. Similarly, an area of approximately 7.2 million hectares was predicted to be unsuitable (Fig. 5). The map shows that the Chure foothills and the surrounding area of eastern Nepal has high potential habitat suitability for D. latifolia and the majority of the area falls in the Parsa National Park and its surroundings. It was found that more than half of the suitable habitat is located only in province 2 ( Table 3). The eastern most province (Province 1) and the province with capital city (Bagmati Province) both have similar areas with higher potential (Table 3). Whereas, western most provinces (Karnali and Sudur Paschim) were found least suitable for growing D. latifolia (Fig. 5).

Discussion
Among the 12 variables considered for predicting the habitat suitability of D. latifolia, elevation was the major driver of its distribution in Nepal. Moreover, mean temperature in the driest quarter, slope and precipitation seasonality also contributed to its distribution. It was found that central and eastern regions with elevation 60-300 m have a higher chance (p > 0.5) of high habitat suitability in natural condition. The two eastern most provinces (Province 1 and 2) represent more than 75% of total suitable habitat for growing D. latifolia in Nepal.
The response curves (Fig. 4) for elevation showed that the highest habitat suitability lies in the range of 200-1000 meters. This is similar to Jackson (1994) who reported D. latifolia to have a widespread distribution in the East Asia-Indian subcontinent up to the altitude of 1000 m in evergreen or deciduous forests with deep, well-drained and moist soils (Soerianegara and Lemmens 1994; Krishnamurthy et al. 2010). Moreover, Praciak et al. (2013) concluded that in India, D. latifolia was wide spread up to an altitude of 1350 m, ranging from the Himalayan tract to southern India. However, the area with the least amount of potential habitat is mostly represented by the forests of higher altitude regions where mean temperature fluctuates in winter and summer (Jackson 1994). Beside this, temperature and precipitation related variables were also found to be important predictors in defining the suitable habitat of D. latifolia. It tolerates minimum temperatures as low as 0-6 °C (Tchoundjeu and Atangana 2011) and the response curve for mean temperature of the driest quarter was consistent with the prediction of best suited temperature as 20-25 °C. Additionally, studies have shown that the annual rainfall in its natural habitat ranges from 750-5000 mm (Louppe et al. 2008) and it can tolerate up to 6 dry months per year with only 40 mm of mean monthly rainfall (Tchoundjeu and Atangana 2011). Our results also showed that the optimum occurrence probability of D. latifolia would be obtained with an average monthly precipitation of around 120 mm and the average minimum precipitation during the coldest quarter was found to be 50 mm. This shows that D. latifolia prefers tropical environments to grow and spread, therefore, the elevation, mean temperature and precipitation of the coldest quarter were the limiting factors for areas suitable for D. latifolia. The low relative importance of non-climatic variables (Forest type) as shown by internal Jackknife procedure and variable importance test indicate that D. latifolia has not any specific preference for the trees it grows in association with and that it can grow in different forest types. Kadambi (1954) reported that it commonly grows with Tectona grandis L.f, Terminalia species, Anogeissus latifolia Roxb. ex DC. (International Plant Names Index 2019) and bamboos in Nepal, however further studies are needed to figure out the effect of forest density and invasion of exotic species for promoting plantation in different densities, particularly in a degraded landscape. The model showed that the forests of the Central region of Nepal, particularly the range of the Parsa National Park, was the most suitable habitat of D. latifolia, whereas the abundance of high-quality habitat decreases from east to west. The density of D. latifolia was reported to be good in the predicted potential habitat in the past, however due to the multipurpose values of the tree for veneers, flooring, moldings, guitar industry, carvings and paneling (Orwa et al. 2009), its occurrence is now becoming very low in its natural habitat. This study only predicted potential habitat of D. latifolia in forest lands and there might also be large areas of potential habitat outside forest lands as well.
Since, Maxent predicts the fundamental niche (different from occupied niche) rather than the realized niche, the predicted potential habitat through Maxent modelling almost always appears over-estimated (Pearson et al. 2007;Kumar and Stohlgren 2009). However, in reality the species might have failed to disperse from its geographic range due to the various anthropogenic and natural barriers (Yang et al. 2013). Although the environmental variables used in this study were limited to only climatic, topography and forest type, they were found able to predict the potential habitat with a high degree of precision. This potential habitat distribution map for D. latifolia can be one of the guiding tools for the forest resource managers and land use planners in order to set priorities for restoration in its natural habitat as well as to initiate different conservation interventions. Beside this, the predicted potential habitat may not be occupied by the species presently, thus the identified habitat can be the best candidate areas to be considered for conservation prioritizations with ex-situ management.

Conclusion
We used a maximum entropy model and based on three sets of environmental variables we assessed the potential suitable habitat of D. latifolia in Nepal. Our results revealed that the environmentally suitable area for D. latifolia in Nepal are limited to the lowland Terai and less elevated Churia of the Central and Eastern regions. Elevation and the Mean temperature of driest quarter were found to be the most important environmental variables to define the habitat suitability for D. latifolia followed by Slope and Precipitation characteristics. The habitat suitability map might be more realistic if we were able to add more environmental covariates, such as soil factors and anthropogenic threats, which will be strongly recommended for future researchers. The study could be useful to optimize D. latifolia distribution by conducting conservation activities in areas of high and moderate suitability where the species is currently lacking.