Spatial distribution and impacts of climate change on Milicia excelsa in Benin, West Africa

African teak (Milicia excelsa (Welw.) C.C. Berg) is an endangered multi-use species. Understanding the impact of climate change on the distribution of this species may improve the ability to anticipate or recognize its decline or expansion and to take appropriate conservation measures if necessary. Ecological niche modeling was projected in geographical space to study the current and future distribution of M. excelsa in Bénin. MaxEnt was used to estimate the potential geographic distribution of the species under two Representative Concentration Pathways (RCP). Miroc 5 summaries and two RCP 4.5 and RCP 8.5 scenarios were used as predictor variables for projections of the geographic potential of this species. The performance of the model was assessed by the area under the curve (AUC), true skill statistics (TSS) and partial receiver operating characteristics (Partial ROC). From the results, M. excelsa was more a secondary species in the Guinean climatic zone and part of the Sudanian-Guinean and Sudanian climatic zone. The projections show a significant decrease in suitable habitats for the species from the two RCP scenarios. Only a part of the Guinean climatic zone remained suitable and few protected areas will conserve in situ M. excelsa. For the sustainable conservation of M. excelsa, it is essential to strengthen the protection of sacred forests located in the Guinean climatic zone.


Introduction
Present global plant extinction rates are 100-1000 times greater than pre-human levels and could increase tenfold by the next century (Pimm et al. 1995;Ricketts et al. 2005). Native species re-introduction is a widely used method to conserve endangered species which often face population growth challenges due to dispersal limitations and transient seed banks (Thompson 1997;Clarke et al. 2007). Despite its widespread use, re-introduction success rates are low due to the difficulty in identifying suitable habitats for restoration (Godefroid et al. 2011;Drayton and Primack 2012;Questad et al. 2014;Kakpo et al. 2018). Therefore, there is a need for methods to identify suitable habitats for endangered species and evaluate environmental factors that influence their distribution (Rovzar et al. 2016).
Knowledge of the spatial and temporal distribution of a species is key for evaluating extinction risks and forecasting possible future threats from factors such as climate change (Pacifici et al. 2015;Urban 2015;Ganglo and Kakpo 2016). Ecological niche models (ENMs) are used to understand ecological requirements of species, predict geographic distributions, select areas for conservation and forecast effects of environmental change Ganglo et al. 2017;Djotan et al. 2018). Studies of geographic distributions of vector species using ENM techniques relate occurrence records and environmental characteristics across species distributions (Altamiranda-Saavedra et al. 2017).

3
The application of ENMs in biodiversity research is well-known for many economically important species in West Africa, especially in Benin. For instance, Fandohan et al. (2015) modeled the ecological niche of Thunbergia atacorensis Akoègn. & Lisowski to investigate how climate change could affect future geographical range of suitable habitats in Benin and Togo. Moutouama et al. (2016) studied the impact of climate change on the range of suitable habitats for Haematostaphis barteri Hook.f. in Benin. Idohou et al. (2016) evaluated the potential spatio-temporal range dynamics for eight economically important wild palms under present and future climate regimes across West Africa. Dotchamou et al. (2016) evaluated climate change effects on the spatial distribution of the African locust bean (Parkia biglobosa (Jacq.) R.Br. ex G. Don) for better conservation management in Benin. Adjahossou et al. (2016)  Milicia excelsa Welw. C.C. Berg, (Moraceae family) is a large tree 30-50 m in height with a diameter at maturity of 1.7-2.0 m, growing in dense forests and forest galleries as well as savannahs (Ouinsavi and Sokpon 2010). In Benin, M. excelsa is commercially known as "Iroko" and is a sacred species, especially within the villages. It is of economic importance adapted to tropical environments but critically endangered in Benin (Adomou et al. 2010). This study estimates current and future suitable areas for M. excelsa and assesses conservation gaps of the species for better conservation in Benin.

Study area
The study was carried out in West Africa with a focus on the Benin Republic located between 6°10′ and 12°50′ N and 1°-3°40′ E (Fig. 1).
The country's climatic profile shows two contrasting zones, (Guinean and Sudanian) and a transitional zone (Sudano-Guinean). The Guinean zone (between 6°25′ and 7°30′ N) is characterized by a sub-equatorial climate with four seasons (two rainy and two dry). The rainfall is approximately 1200 mm per year, beginning mid-March to July and a second period from September to November. Temperatures vary between 25 and 29 °C, and relative humidity between 69 and 97%. The Sudanian zone (9°45′-12°25′ N) has a tropical dry climate with two seasons of equal length (rainy and dry). Mean annual rainfall is often less than 1000 mm and occurs mainly from May to September. Relative humidity varies from 18 to 99% and temperatures from 24 to 31 °C. The transitional Sudano-Guinean zone (7°30′ to 9°45′ N) has two rainy seasons merging into a unimodal regime. Annual rainfall fluctuates between 900 and 1110 mm, temperatures are between 25 and 29 °C and relative humidity from 31 to 98%.

Data collection
A total of 1856 records, (georeferenced occurrences), were downloaded from the Global Biodiversity Information Facility portal (GBIF 2016). After cleaning, 209 occurrence points were used for the West African region.
Current  climate data were obtained from WorldClim version 1.4 (available at www.world clim.org/ biocl im). Future climate projections for 2055 were downloaded from AFRICLIM version 3.0 (available at www. york.ac.uk/envir onmen t/resea rch/kite/resou rces/; Platts et al. 2015). For future climatic conditions, predictions from the Interdisciplinary Research on Climate Change (MIROC5) were used. The projections were run under Representative Concentration Pathway (RCP) 4.5 and RCP 8.5 for the 2055-time horizon. By the mid-twentyfirst century, RCP 4.5 projects temperatures to rise above industrial levels by at least 1.5 °C in West Africa, with atmospheric CO 2 reaching 500 ppm (IPCC 2013). Under the more extreme RCP 8.5 projections, temperatures are predicted to rise by 2.8 °C and atmospheric CO 2 to be over 550 ppm (IPCC 2013). RCPs are third generation scenarios and are preferred to the Special Report on Emissions Scenarios (SRES) because they allow more flexibility (and reduced costs) in modeling processes (van Vuuren et al. 2011). RCPs imply collaboration between impacts, adaptation, and vulnerability research, and climate and integrated assessment modeling (van Vuuren and Carter 2014). These scenarios have been developed to explore different combinations of scenario context (demographic, socioeconomic, land use, and technology) (Moss et al. 2010). RCP 4.5 and RCP 8.5 have already been used for studying West African ecosystems (Ayihouenou et al. 2016;Gbètoho et al. 2017). These climate projections were statistically downscaled to match the bioclimatic variables using the delta method (Peterson and Nyari 2008; Ramirez-Villegas and Jarvis 2010).

Modeling of species and model evaluation
The maximum entropy species distribution model algorithm (MaxEnt, version 3.3.3 k) was used for habitat suitability modeling. During the modeling process, occurrence data were cleaned up by removing duplicate records in grids and records which did not correspond to the description of its locality (country) in order to reduce sampling bias which may result from the oversampling of some sites in the study area (Elith et al. 2006).
To optimize model performance, several parameters in MaxEnt were tested. These parameters are regularization multiplier, max number of background points and default prevalence (Fig. 2). MaxEnt models were developed using: regularization of 10,000 background points, regularization multiplier of 0.5, default prevalence of 0.7, and 10 replicates. A jackknife test was performed on the selected bioclimatic variables to determine the ones that best contribute to the model prediction. Owing to common collinearity and non-independence of climate dimensions (Zuur et al. 2010), we examined correlation patterns among variables to select those not closely correlated using environmental niche modeling tools (ENMTs). Hence, we included only a subset of variables with Pearson correlation coefficients below 0.80 (Elith et al. 2010). To avoid geographic bias, a bias surface using a kernel density estimate was generated (Galdino et al. 2016). This resulted in a raster where cells with lower values represented places with lower bias. The bias surface was used to account for sampling intensity and potential sampling bias (Jarnevich et al. 2015). All maps were converted to binary via a conservative least presence thresholding approach consisting of the lowest predicted value (0.405) corresponding to an occurrence record of the species in the calibration data set (Barve et al. 2011). All GIS-related work, (layers clipping and conversion, map composition), was carried out in QGIS 2.18.0.
The performance of the model was assessed using the area under the receiver operating characteristic curve (AUC) (Elith et al. 2006), true skill statistic (TSS) (Allouche et al. 2006;Elith et al. 2006) and Partial ROC ). The AUC is the probability that a randomly chosen presence point of the species will be ranked as more suitable than a randomly chosen absence point (Elith et al. 2006). A model is considered as having a good fit when its AUC is close to one (AUC ≥ 0.75) (Elith et al. 2006). The TSS is the capacity of the model to accurately detect true presences (sensitivity) and true absences (specificity). A model with TSS ≤ 0 indicates a random prediction, while one with a TSS close to 1 (TSS > 0.5) has a good predictive power (Allouche et al. 2006). For Partial ROC, the evaluation dataset was bootstrapped and probabilities obtained by direct count of the area under the curve (AUC) ratios falling ≤ 1 via a Visual Basic script (https ://shiny .conab io.gob.mx:3838/niche toolb 2/), with 100 iterations (Altamiranda-Saavedra et al. 2017). Model performance was assessed by dividing the occurrence data of M. excelsa randomly into calibration (70%) and evaluation (30%) subsets.

Model validation
The correlations analysis and Jackknife AUC test identified five bioclimatic variables as contributing more to the model across the study region: Isothermality, Annual Precipitation, Temperature Annual Range, Mean Diurnal Range, and Precipitation of Wettest Month (Fig. 3).
Model evaluations indicated that the model was robust (AUC = 0.979). The TSS value was 0.86 (Fig. 4) and AUC ratios were well above 1.0 (Fig. 5). Therefore, the model showed excellent performance.

Suitable habitats in present and future
The analysis revealed that about 54,009 km 2 corresponding to 47.1% of Benin's surface area (114,763 km 2 ) is currently   Fig. 6). Only a part of the Guinean climatic zone will remain suitable. The protected areas of Benin illustrate a good potential to conserve actual or current populations of M. excelsa. In fact, many of the protected areas had suitable habitat for M. excelsa in all climatic zones. In the future, according to the two scenarios, only the Lama and Djigbé forests will remain suitable for the conservation of the species.

Credibility of the model
In light of global climate changes, it is essential to identify species with high risk of extinction and to diagnose risk driver factors to reverse the decline of biodiversity (Darrah et al. 2017). Information on changes to species ecological niche is important for assessing species geographic range and their realized niche (Breiner et al. 2017). Among several programs for ecological niche modeling and habitat suitability prediction, MaxEnt is one of the most used with respect to our type of data, presence-only data. This modeling tool using presence-only data is one of the best performing algorithms among those using climate modeling approaches (Phillips et al. 2006), and is relatively robust for small sample sizes (Pearson et al. 2007). Max-Ent is a machine learning method that estimates species distribution across a study area by calculating the distribution probability of maximum entropy, subject to the constraint that the expected value of each feature under this estimated distribution should match its empirical average (Phillips et al. 2006). However, MaxEnt cannot indicate the realized niche because of its weaknesses, including uncertainties related to the projections (Elith et al. 2006;Schwartz 2012), failure to account for relevant variables or for geographic barriers , bias due to unequal sampling (Fourcade et al. 2014), overfitting, and issues related to evaluation data and methods (Radosavljevic and Anderson 2014). However, the migration of species is a determinant factor in the potential  impact of climatic change on their habitats (Elith et al. 2006). When constraints of dispersion are not considered, temperature increase could expand the distribution area of certain species (Sharma and Jackson 2008;Buse and Griebeler 2011). In the current study, indirect variables such as soil and vegetation were not considered. According to Guisan and Zimmermann (2000), direct parameters such as temperature and precipitation are more efficient when the modeling is over large areas. This is contrary to indirect parameters which are not more efficient but are more inclined to introduce errors into the model.

Ecological niche modeling and conservation of M. excelsa
Ecological niche modeling of M. excelsa showed that at present the species is more marginal in the Guinean climatic zone and parts of the Sudanian-Guinean and Sudanian climatic zones. Although the species can be found in dry zones, it reaches its ecological optimum in wetlands, mainly in dense forests and gallery forests (Ouinsavi and Sokpon 2010). In each climatic zone of Benin, at least one protected area is located in a suitable habitat of the species. In the Guinean climatic zone, there are the forests of Lama, Agrimey, Djigbé, and Sakété. In the Sudanian-Guinean climatic zone, there are the forests of Dan, Atchérigbé, Dogo-Ketou, Pénéssoulou, Ouémé Superieur, Boko, N'Dali, Sérou and Sakarou, and in the Sudanian climatic zone, there are the forests of Donga, Kilir, Soubroukou, Birni, Mékrou, and Atchérigbé. Therefore, the protected area network of Benin conserves the species.
In contrast to the present species distribution, future projections show a significant decrease in suitable habitat under RCP 4.5 and RCP 8.5. Only Lama and Djigbé forests retain suitable habitats under both scenarios. This decrease in habitat suitability may be explained by the significant changes projected for bioclimatic parameters, mainly precipitation and temperature. M. excelsa is adapted to temperatures ranging between 25 and 35 °C and precipitation between 1150 and 1900 mm per year (Daïnou et al. 2012). According to Busby et al. (2010), fluctuations in climate variables such as precipitation and temperature will have an impact on biological diversity and on the geographical distribution of suitable habitats. Doxa et al. (2017) suggest that protected areas a principal tool of in situ conservation of biodiversity. Based on future projections, this study shows that few protected areas are key to the conservation of M. excelsa. Therefore, it is important to take the action to conserve effectively M. excelsa by including sacred forests of the Guinean climatic zone in the protected area network. M. excelsa is a sacred tree and is found in many sacred forests.

Conclusions
This study estimated current and future distribution of suitable habitats for M. excelsa so as to assess conservation needs of the species. The suitability prediction will decrease in Benin area under both representative concentration pathway scenarios. At present, M. excelsa is suitable for the Guinean climatic zone and portions of the Sudanian-Guinean and the Sudanian climatic zones. Under future projections, the area of suitable habitat decreases. Among all protected areas, only the Lama and Djigbé reserves will still provide suitable habitats for the species.
Our results provide scientific rationale for planning for and implementing M. excelsa conservation. Despite the considerable potential of M. excelsa and its socio-economic attributes, the species resources are still obtained from the wild. M. excelsa has not been included in various agroforestry and conservation programmes of Benin. Currently, the species is endangered due to over-exploitation, deforestation, numerous human activities and climate change. However, for its sustainable management, forest managers should strengthen the protection of the Lama and Djigbé forests, especially the natural areas, by constructing firewalls and by organizing regular night patrols. Special attention should be given to the sacred forests located in the Guinean climatic zone by making enrichment planting of the species and classifying them as protected area. In addition, the species can be used in agroforestry systems during re-forestation campaigns in the Guinean climatic zone.
Although ecological niche modeling can estimate habitat suitability, it cannot predict species response in term of adaptation, dispersion, migration, or extinction. Thus, suitable habitats inside a protected area does not imply that species will migrate there naturally because the reaction of competing species is not considered. Future research on the conservation of M. excelsa should consider its phenotypic diversity since species conservation implies the conservation of all its genetic diversity.