Potential and Future Geographical Distribution of Eremanthus erythropappus (DC.) MacLeish: a Tree Threatened by Climate Change

Eremanthus erythropappus is a commercially-important tree which has a long history of exploitation in the Brazilian State of Minas Gerais. The knowledge on the potential geographical distribution of E. erythropappus is therefore critical for the species sustainability. Thus, the aim of this study was to estimate and map current and future ecological niche for E. erythropappus in Minas Gerais. We used the Random Forests algorithm to model the ecological niche for current and future climates scenarios. Our predictions indicate Espinhaço, Mantiqueira, and Canastra mountain ranges as core areas of distribution and forecast drastic reductions in potential areas under all climate scenarios. Based on our results, we highlight that the continual harvesting of naturally-occurring E. erythropappus populations will not be sufficient to supply the market demand. Silviculture practices would likely serve as an economically viable and ecological sustainable alternative to harvesting natural populations.


INTRODUCTION
Eremanthus erythropappus (DC.) MacLeish, (commonly-known as Candeia), is a commercially important tree which occupies transition areas between the Brazilian Atlantic Rain Forest and Savannah Biomes. The bulk of the naturally-occurring populations of the species in South America are within the state of Minas Gerais, Brazil, particularly in upland areas (> 9 00m) (Silva et al., 2008;Clark et al., 2011).
Eremanthus erythropappus trees are light-demanding ecotone specialists often found in open grasslands and edge of forests, and rarely under the forest canopies (Oliveira & Fluminhan, 1999;Araújo et al., 2017). Habitats with a high dominance of E. erythropappus (≥ 70% frequency) trees are locally known as Candeais. These habitats often represent a transitional vegetation zone between closed forest and open savannah habitats (Oliveira & Fluminhan, 1999;Araújo et al., 2017).
This tree species has a long history of exploitation, either for use as fence posts, firewood, and on a larger economic scale, for essential oil extraction (Clark et al., 2011;Donadelli, 2012;Scolforo et al., 2016a). In the case of essential oil extraction, the targeted component is alpha-Bisabolol (Scolforo et al., 2016a), which is a key ingredient for skincare and global cosmetics industries (Vieira et al., 2012;Scolforo et al., 2016b). Despite the availability of synthetic substitutes of alpha-Bisabolol (Clark et al., 2011) and also the lengthy period of tree growth, which spans at least 12-15 years (Pérez et al., 2004;Mori et al., 2009), naturally-occurring populations of E. erythropappus are still important sources of alpha-Bisabolol for the global market (GFA Consulting Group, 2006), comprising 80% of the alpha-Bisabolol exported (Clark et al., 2011).
Before regulations on E. erythropappus extraction, the species was often exploited in an unsustainable fashion. In 2004, regulations were placed on E. erythropappus exploitation in the Minas Gerais State. In 2007 further regulations were implemented on the quantity of E. erythropappus that may be harvested and types of exploitation management systems (Scolforo et al., 2012).
Nevertheless, alpha-Bisabolol is still being purchased in the European market without any requirement of traceability. To date, only two companies, Atina and Citróleo, possess the Forest Stewardship Council certification regarding suitable management of E. erythropappus populations (Donadelli, 2012). Also, the increasing global demand for natural alpha-Bisabolol has not been matched by adequate or sufficient mechanisms of regulation and governance. Thus illegal extraction and unsustainable harvesting practices of the species is still rampant (Clark et al., 2011;Donadelli, 2012).
Due to these factors, E. erythropappus is increasingly being threatened by human activities. Additionally, there is the risk, as yet unexplored, of climate changes negatively affecting the future distribution and survival of the species. It is well-known that each species has its climatic niche, and species within their climatic niche will have a better chance of remaining healthy, productive and able to maximize their ecological and/or economic value under changing climate (Wang et al., 2016). Species growing outside their climatic niche are therefore likely to have limited productivity and to be more vulnerable to disturbances (Hamann & Wang, 2006;Kurz et al., 2008).
The knowledge of the potential and future geographical distribution of E. erythropappus is therefore critical for the sustainability of native populations, as well as to inform silvicultural practices to ensure productivity of E. erythropappus plantations. The aims of our study are, therefore, to: (1) model and map the current ecological niche of E. erythropappus in the state of Minas Gerais, identifying the core areas of E. erythropappus occurrence; (2) predict and assess changes in the potential distribution of E. erythropappus under scenarios of climate change by the years 2050 and 2070; and (3) access the adequacy of conservation units for protecting the current and future populations of the species.

Study region
The study region is the Brazilian State of Minas Gerais, with a land area of over 586,000 km 2 , equivalent to the area of countries such as France and Spain ( Figure 1). Minas Gerais harbors the largest occurrence of E. erythropappus populations, and the highest number of extractive industries of alpha-Bisabolol (Donadelli, 2012). The region also encompasses a wide altitudinal range (between 40 and 2,600 meters) and five climate classes according to the climate classification system of Koppen (Martins et al., 2018). This variability gives rise to the occurrence of different vegetation biomes (Savannah, rainforest and semi-arid woodland) and the diverse phytogeography within the region. Savannah covers 57% of State (west/northwest), with warmer and wetter climate during summer, and a pronounced dry period. Rainforest, locally known as Atlantic forest, dominates the south and east regions, occupying over 41% of State, and features wetter and milder climates, particularly in the south. Finally, semi-arid woodland comprises only 2% of the state, predominantly in the north, and is characterized by semi-arid and sub-humid climates, with higher temperatures and lower rainfall throughout the year. Oxisols represent the predominant soil type in the state, with areas of Cambisols and Quartzarenic/Litholic Neosols (Curi et al., 2008).

Species distribution modeling
We used an ecological niche modeling (ENM) approach to model the potential distribution of E. erythropappus under climate change scenarios. This approach has been widely used for tree species, either to verify the vulnerability and impacts of climate change in order to maintain the species sustainability (Pouteau et al., 2012;Turchetto-Zolet et al., 2016), or as a tool for guiding recovery plans and indicating potential plantation areas (Coelho et al., 2016;Carvalho et al., 2017). The ecological niche models are based on Hutchinson's fundamental niche and its corresponding abiotic component known as the ecological niche (Sóberon & Peterson, 2005). The modeling process allows us to use abiotic conditions observed in the known geographic distribution of our target species (realized niche) to predict potentially adequate environments where the species is theoretically viable. By projecting our model results onto a map, we can then visualize regions with potential distribution of the species at present or in the future.
Several techniques are available for ENMs. Among them the Random Forests algorithm (RF) stands out for its performance and simple parametrization (Cutler et al., 2007;Wang et al., 2016). The algorithm is considered the best-performing statistical approach for mapping vegetation and for ecological niche model building (Garzón et al., 2006;Lorena et al., 2011;Wang et al., 2012), and is being adopted by various countries for biodiversity modeling and for informing government politics on global warming (Wang et al., 2016;Prasad et al., 2016).

Occurrence data
Presence and absence observations for E. erythropappus in Minas Gerais state were obtained from several database providers. Presence data was obtained from 324 sustainable management plans of native E. erythropappus populations. These data were collected at environmental agencies within the state. To complement this dataset, we obtained an additional 372 and 294 observations of presence data points from the online databases speciesLink and Global Biodiversity Information Facility (GBIF), respectively. We also used field data from 197 fragments of native vegetation, where E. erythropappus were present in 36 and absent in 161 fragments. In total therefore, we compiled 1,187 points (1,026 presences observations and 161 absences). These points were rasterized at the spatial resolution of 0.008333 arc min (approximately 1 km), which was the same resolution we used to rasterize climatic variables. Many of these observations were within 1 km of each other, belonging to the same raster pixel, and therefore we engaged in a process of data "cleaning". To accomplish this, each pixel was assigned a presence/absence status of E. erythropappus, resulting in 452 pixels with presence and 160 pixels with absence observation status ( Figure 1). We then used this final 612 presence and absence observations to extract environmental variables for input into the subsequent modeling process.

Environmental data
The historical and future climate data were sourced from the WorldClim online database (Hijmans et al., 2005) using the finest spatial resolution of 30 arc-seconds (~1 km) available. For modelling, we used seven climatic variables (Table 1) relating to temperature and rainfall and corresponding to the current period   (Carnaval & Moritz, 2008): BIO1 -Annual Mean Temperature; BIO4 -Temperature Seasonality; BIO10 -Mean Temperature of Warmest Quarter; BIO11 -Mean Temperature of Coldest Quarter; BIO12 -Annual Precipitation; BIO16 -Precipitation of Wettest Quarter; and BIO17 -Precipitation of Driest Quarter. According to Bucklin et al. (2014), climate predictors alone are an effective and efficient approach for initial assessments of environmental suitability, whereas additional predictors have relatively minor effects on the accuracy of predictions. Nevertheless, we incorporated an additional three topographical and soil features as predictive variables (Table 1), as these variables have been documented as environmental drivers of E. erythropappus distribution within the region (Pérez et al., 2004;Scolforo et al., 2012). Soil texture with original scale of 1:500,000 was categorized into three classes (fine, medium and coarse). Information of slope was obtained from the Shuttle Radar Topography Mission (SRTM) with 90 m of resolution, and was classified into four classes: flat or gently sloping, sloping, moderately or strongly sloping, and steep or very steep. These variables were derived from the Ecological Economic Zoning of Minas Gerais (Curi et al., 2008) and rasterized at 1 km spatial resolution. Altitude (meters above sea level) was also obtained from WorldClim at the same spatial resolution (~1 km).  We adopted two representative concentration pathways (RCPs) from IPCC Fifth Assessment Report. The RCP4.5 assumes growth in the greenhouse gas emissions trajectory is limited through several initiatives, while RCP8.5 corresponds to a worst-case scenario, with the highest estimated magnitude of greenhouse gas emissions where no climate specific mitigation targets are set (McIntyre et al., 2017). These scenarios were also obtained from the WorldClim using the spatial resolution of 30 arc-seconds (~1 km). We configured our models to generate predictions for the years 2050 (average for 2041 to 2060) and 2070 (average for 2061 to 2080). For each future period, there are four different scenarios of climatic conditions (2 RCPs × 2 GCMs). For ease of interpretability and in order to incorporate the uncertainties related to projections of climate change (Wang et al., 2012;Zhang et al., 2015), we averaged the results of these four projections outputs at each period ( Figure 2).

Habitat suitability modeling and evaluating
We used the Random Forests (RF) package (Liaw & Wiener, 2012) in the R platform (R version x64 3.4.1) (R Core Team, 2017) to model the relationship between abiotic variables and current known locations of E. erythropappus occurrence. We then used future climate data in the years 2050 and 2070 to generate predictions of the future E. erythropappus habitat niche.
Usually, when the database has unbalanced data (e.g. more data points for presence than absence), the results may be overestimated for the frequency class. To overcome this problem, we used the methodology of Wang et al. (2016) and employed an ensemble of 100 RF models, each of which was built with randomly sampled data points for presence, while the data points for absence remained the same. The final prediction was the most voted class or mean of probabilities of all RF models.
After tests, we set 500 decision trees in each forest and used the default, square-root of the number of predictors at each split. We assessed the predictive RF accuracy by the percentage of correctly classified data (overall accuracy), and the probability of presence and absence being correctly predicted (sensitivity and specificity) through out-of-bag data (OOB). We also used the importance values generated by RF (based on Gini Index) to identify the abiotic factors important for determining the ecological niche of E. erythropappus (Menze et al., 2009). Subsequent to training, we inputted the final RF models with the environmental data of the study area to generate predictions for the current period and projections for future periods (2050 and 2070).
Our ENM predicted areas currently occupied by the species and also areas environmentally suitable for the species, but with no known records of occupancy due to various factors, including physical barriers, seed dispersal, migration, and/or human interference (Wang et al., 2012). This should be interpreted as the quality of the site considering the species preferences, and we interpreted here as habitat suitability (HS). Because the classification of probability in binary data comparisons of range area changes for different scenarios, it is required a threshold to classify probabilities in presence and absence. We set three thresholds for classifying the habitat suitability of E. erythropappus, considering the species as at a given grid if the habitat suitability (HS) is above or equal: 0.3 -broadly distributed; 0.5 -moderately distributed; and 0.7 -restricted distribution.
To track changes in latitudinal distributions, we compared geographic centroids of predictions and projections species range (Zhang et al., 2015). The geographic centroid weighted by probability was calculated using the mean center function in ArcGIS10.1 (ESRI, 2010). This function calculates the average coordinates weighted by probability in the study area through the Formulas 1 and 2, where i x and i y are the central coordinates for grid i, n is equal to the total number of grids, and i w is the probability at the grid i.

Assessment of ecological niche and important abiotic variables
Our environmental niche models for Eremanthus erythropappus presented reliable results, with an overall error rate average of 17.2%, and presence/absence error rate of 17.5% and 17.2%, respectively. As expected, the OOB error rate was similar for presence and absence after the sampling rate was balanced for both classes through bootstrapping. Amongst the ten environmental input variables, precipitation of Wettest Quarter (BIO16), mean temperature of coldest quarter (BIO11), annual mean temperature (BIO1), and altitude were the four most influential factors according to RF importance values (decrease impurity values of around 20%; Table 1). In contrast, the categorical variables of slope and texture were the least influential environmental variables (Table 1).

Current potential suitable habitat for Eremanthus erythropappus
According to our environmental niche models, E. erythropappus is currently distributed in the central-south portion of the state, with a particular concentration in mountainous environments (Figure 3).
The habitat suitability thresholds revealed that the most suitable areas for the occurrence of E. erythropappus is the Espinhaço, Mantiqueira and Canastra mountain ranges, which we regard as core areas for the species. The area covered by habitat suitability greater than 0.3, 0.5 and 0.7 corresponded respectively to 158,569 km 2 , 61,467 km 2 and 24,188 km 2 (Table 2).

Ecological niche of Eremanthus erythropappus under climate change scenarios
We found varying degrees of net losses and gains in ecological niche for E. erythropappus under different modelled habitat suitability (HS) thresholds. The loss in suitable habitat is particularly drastic in 2050 and 2070 under habitat suitability (HS) thresholds of 0.5 and 0.7 (Figure 4b, c, e, f; Table 2). Under HS ≥ 0.5, suitable habitat for E. erythropappus in 2070 is found only in core regions, such as Espinhaço and Mantiqueira (Figure 4b, e), and has vanished in Canastra mountain range. On the other hand, presence probabilities for E. erythropappus under the ≥ 0.7 threshold shows serious decline in the state, covering only a small area in the Espinhaço and Mantiqueira mountain ranges in 2050 (Figure 4c) and practically extinguished in 2070 (Figure 4f).
The substantial gain in areas of suitable habitat in 2050 and 2070 in relation to the current conditions appeared only when we adopted the HS ≥ 0.3 threshold. The average area for HS ≥ 0.3 was 158,569.00 km 2 (2018), decreasing to 123,155.00 (2050) and increasing to 192,608 km 2 in 2070 (Figure 4a, d).
We also detected potential changes in distance and direction of E. erythropappus occurrence mean centers (Figure 4g and h) where all mean centers were located in the Espinhaço mountain range. Under the HS ≥ 0.3 threshold, there was a shift in suitable habitat towards to the west region. But under HS ≥ 0.5, the shift was towards the south region. We did not calculate the mean center for 2070 under HS ≥ 0.7 due to the scarcity of remaining grids under this scenario.
Minas Gerais State has 103 conservation units encompassing a wide range of vegetation structural types. However, only half of them (58 units) are suitable E. erythropappus habitat (HS ≥ 0.5 threshold). These conservation units encompass a land area of 4,087.37 km 2 of protected suitable habitat for the species (Figure 5). We observed a 74.3% reduction in suitable habitat area for the species within these conservation units in the year 2070, with only 1,049.96 km 2 (Table 3) of suitable habitat remaining.

DISCUSSION
The bioclimatic variables BIO16 (precipitation of wettest quarter), BIO11 (the temperature of coldest quarter), BIO1 (Annual Mean Temperature) have the strongest bearing on our ecological niche models. According to Condit et al. (2000) and Brenes-Arguedas et al. (2011), this rainy period (BIO16) is crucial for plant photosynthesis, growth and resource allocation to reproduction, which affects the formation and quality of E. erythropappus seeds (Feitosa et al., 2009). The mean temperatures, annual and the coldest quarter, indicate that E. erythropappus favors mild temperatures but does not tolerate ground frost in the winter (Magalhães et al., 2008;Siqueira & Pinto, 2009).
Altitude is the fourth most important variable and is of great importance for E. erythropappus trees in our model. Consistently, the upland areas have Table 2. Area (km 2 ) of current and future suitable habitat of Eremanthus erythropappus under three different habitat suitability (HS) thresholds: ≥ 0.3 -broad distribution; ≥ 0.5 -moderate distribution; and ≥ 0.7 -restricted distribution. Where "Loss" represents the contraction area of the ecological niche of the species in the future scenario in relation to the current scenario; "Gain" represents the ecological niche expansion area in the future scenario in relation to the current scenario; and "Stable" represents the area covered by the ecological niche of the species both in the current scenario and in the future scenario. Total area (%) was calculated through the ratio of the area predicted in the future under the area predicted in the present multiplied by 100.  (Pérez et al., 2004;Silva et al., 2008;Clark et al., 2011;Scolforo et al., 2012). It is worth noting that lowland or valley sites usually are not dominated by the tree species, where inter-species competition is likely the main reason for its exclusion. However, E. erythropappus occasionally may colonize these fields after deforestation (Scolforo et al., 2012;Pádua et al., 2016) and bushfires during the warmer months of the year (Oliveira & Fluminhan, 1999;Araújo et al., 2017).
The current and future suitable habitat of E. erythropappus is crucial for informing policies on the next steps in the conservation and sustainable management of the species. Eventually, these concerns will need to be addressed if climate change starts to affect the productivity of the species in certain areas, resulting in reduced population or lowered productivity. We therefore identified three core areas with high potential for the conservation of the species, namely areas within Mantiqueira, Espinhaço and Canastra mountain ranges. It is also noteworthy that these mountain complexes represent some of the most ancient open vegetation areas in eastern South America, and are of high conservation priority (Silveira et al., 2016). In fact, it is recognized that upland areas are less affected by climate changes and can still guarantee the species habitat. According to Gwitira et al. (2014), high-altitude areas will remain as biodiversity hotspots of savannah domain. Forest loss in warmer regions tends to accelerate the movement of species upslope, especially in the tropics, where both climate and habitat loss stressors will increase the risk of net lowland biotic attrition (Guo et al., 2018).
Based on the most optimistic scenario of our models (HS ≥ 0.3), these three core areas could maintain areas that will remain connected as a huge network. However, the pessimistic scenarios (HS ≥ 0.5) for 2050 and 2070 show limited gains in suitable area and the loss of continuous distribution of the species between the three core mountain areas. Most worryingly, we detected a huge reduction in suitable habitat area the species in 2070 under the most pessimistic scenario (HS ≥ 0.7), where there will be practically no suitable habitats for the species in Minas Gerais.  Table 3. Predicted extent of suitable habitat area for Eremanthus erythropappus within conservation units in Minas Gerais, Brazil at a threshold HS ≥ 0.5 (Total protected area). The classes of loss in area (%) were estimated based on the contraction area of the ecological niche of the species within each conservation unit in relation to the total area of the conservation unit, and represents de number of conservation units for each class of loss. "Gain" corresponds the number of conservation units which presented expansion area of suitable habitat for the species.  In contrast to losses in suitable area, we did detect some gains in habitat suitability in the west of the state, in the Canastra mountain. However, this happened only when we adopted a HS ≥ 0.3. These gains may be attributed to increasing precipitation of wettest quarter (BIO16) in this region, which is forecasted in climate scenarios ( Figure 2).
Currently, conservation units within Minas Gerais state cover a considerable area of the ecological niche of E. erythropappus. Several of our future predictions have shown how deep the negative impacts on nature reserves are. The losses in suitable area and connections may lead to reduced genetic diversity in the species, thereby amplifying the risk of the species' extinction. Genetic studies indicate that even though E. erythropappus possess high genetic variability, trees typically have significant co-ancestry, where individuals located in close proximity to each other are genetically similar (Pádua et al., 2016). For this reason, conservation and restoration strategies efforts should be done at protection units around the three core distribution areas.
It is therefore timely to consider strategies for exploitation use of E. erythropappus that will be socially acceptable, economically viable, and ecologically sound. Since climate changes will impact negatively natural populations of the species, silviculture practices would likely serve as an economically viable and ecological sustainable alternative to harvesting natural populations.
Several studies have already been published regarding the sustainability of plantations for oil extraction from E. erythropappus (Feitosa et al., 2009;Silva et al., 2014;Scolforo et al., 2016a, b). Practically, plantations of E. erythropappus in the coming decades may be most feasible in the western parts of the Canastra mountain ranges, where we predicted some potential gains in suitable habitat for the species.
Additionally, E. erythropappus trees are easy to be propagated even in many soil types and conditions due to its ecological preferences (Amaral et al., 2015;Meira et al., 2015;Pereira et al., 2015). We recommend therefore that management initiatives spearheaded by the Brazilian government should encourage land owners to cultivate this species for commercial use or forest restoration areas. At present, Minas Gerais has 25,000 km 2 of rural areas and 4,000 km 2 of legal reserves under the Brazilian Forestry Service (SFB, 2019) which could benefit from reforestation programs incorporating this tree species (February 2018).
In conclusion, Eremanthus erythropappus is an ecologically and economically important tree species in Brazil, but our projections on how climate change couldaffect it is of concern for the sustainability of future natural populations of the species in the State of Minas Gerais. We suggest that it may be prudent to engage in activities that will ensure sustainable harvesting of E. erythropappus in Minas Gerais State, and as an insurance against the potential negative impacts of climate change on the species.