Using ecological niche models to predict the impact of global climate change on the geographical distribution and productivity of Euterpe oleracea Mart . ( Arecaceae ) in the Amazon

We assess the impact of climate change on the geographic distribution and productivity of Euterpe oleracea (Arecaceae), commonly called açaí. To construct the ecological niche model of E. oleracea, we used 95 points of occurrence, fi ve bioclimatic variables in current and future climate scenarios and Maxent software. Th e Akaike Information Criteria (AIC) was used to rank the ability of the models (considering ecological, socioeconomic and spatial variables) to explain the variation in productivity of E. oleracea among 200 municipalities. Th e maps showed that regions with higher environmental suitability for E. oleracea were concentrated in northern and northeastern Brazil, which was similar to the spatial pattern of productivity data of E. oleracea. Future climate conditions tend to promote an increase in the geographical distribution of this species, even though the new regions are in the arch of Amazon deforestation. Only space and the environmental suitability (indicated by the ecological niche model) were important for explaining the productivity of E. oleracea. Th at is, municipalities that are more productive are located in more suitable environmental regions. Th erefore, it is important to use niche models to explain demographic changes and to estimate species demographic attributes.


Introduction
Global climate change has aroused much interest and concern in the mass media and academia, which refl ects a growing number of papers published on the topic in recent years (Nabout et al. 2011).Recent studies have investigated the impact of global climate change on biodiversity, food security, economic impacts (e.g.Tol 2009), and socioenvironmental vulnerability (e.g.Torres et al. 2012).Th ese studies have shown that global climate change may reduce the area that is climatically suitable for focal species or shift their geographical distributions (mainly studies in Brazil, see Diniz-Filho et al. 2009a), which may also lead to economic losses and increases in socio-environmental vulnerability (see Tôrres et al. 2012).In parallel, some studies have focused on developing mitigation strategies, such as the selection of conservations areas for future climate scenarios (e.g.Loyola et al. 2013) and international agreements for reducing greenhouse gases (for example, the Kyoto protocol).
In the case of plants, global climate change can affect geographic distribution, life cycle, biological interactions, and many others factors (see Thuiller et al. 2008).Moreover, global climate change can affect the productivity of many plant species, especially species that undergo harvesting and extractive practices (Nabout et al. 2011;2016).Such species may suffer progressive changes in their production in coming years, particularly in the Brazilian Amazon region, which is a very important part of the local and rural economy (Homma 2012).
Among the extractive products with large natural stocks, the fruit of the açaí palm (Euterpe oleracea) is considered particularly important (Homma 2012).Belonging to the family Arecaceae, E. oleracea has a predominantly tropical and subtropical distribution.In 2014, E. oleracea fruit became one of the most important non-wood extractive activities, with highest value of production (in Real) among other plants (Brasil 2014).The majority of E. oleracea fruit production (about 80%) has its origin in extractive harvesting.Only 20% comes from the areas that are managed and grown in the lowland and upland regions (IBGE 2016).
Experimental, empirical, or mathematical models can been used to evaluate the impact of global climate change on biodiversity (Nabout et al. 2011).Recently, the use of Ecological Niche Models (ENM) has increased in the scientific literature and has proved a very effective tool in various areas of science (Peterson et al. 2011;Vaz et al. 2015).For example, ENM has been used to predict the spatial abundance, genetic variability, spatial distribution, and extinction of species (Nabout et al. 2011;2016;Peterson et al. 2011;Tôrres et al. 2012;Diniz-Filho et al. 2015;Soares et al. 2015).
Besides climate (obtained by ENM), socio-economic factors may be able to explain the variation in productivity, abundance, or genetic variability among sites (Diniz-Filho et al. 2009b;Soares et al. 2015).For plants of economic interest, it is expected that high productivity regions occur in municipalities with lower socioeconomic development (see for example, Nabout et al. 2012) and better environmentally protected regions (Nabout et al. 2016).
Here, we investigate the effects of climate change on the geographic distribution and productivity of E. oleracea in 200 Brazilian municipalities.In addition, we also determined the relative importance of ecological (ENM) and socioeconomic factors in the total harvested production of E. oleracea in 200 Brazilian municipalities.

Productivity data
Productivity data for E. oleracea were obtained from the Brazilian Institute of Geography and Statistics (IBGE, www.ibge.gov.br) in 2013, using their website's automatic recovery system (SIDRA) and "Vegetal Extraction" tab.We identified a total of 200 Brazilian municipalities from all five geographic regions in which the production of E. oleracea had been recorded over ten years period (2002 -2011).For each municipality, we calculated the average annual production, and as these municipalities have different geographic areas, we also measured the average production of each municipality by its geographic area, according to IBGE.Thus, the average production of E. oleracea fruits per municipalities is given in tons per hectare (ton.ha - ).

Socioeconomic and ecological data
For each of the 200 E. oleracea producing municipalities, we obtained the follow socio-economic data: i) Gross Domestic Product (GDP) per capita (Brazilian Real; R$), ii) Municipal Human Development Index (HDI), iii) Percentage of Conservation Units (%CONSERV) by municipalities, and iv) Percentage of Rural Population (%RURAL).We then used these data to build statistical models to explain the average productivity of E. oleracea in each municipality.We expected that the HDI and GDP variables would have a negative relationship with productivity (i.e., a higher productivity rate in poorer municipalities).In contrast, we expected that the %CONSERV and %RURAL variables would show a positive relationship (i.e., higher production in the best preserved regions, comprising a larger proportion of the rural area).
We obtained GDP data per capita (R$), Municipal HDI, and total population data in the city and rural population (for the year 2007) from the "Download channel electronics" tab IBGE database.We obtained the %CONSERV and %RURAL variables for each municipality from the database of the Socio-Environmental Institute and Monitoring Program of Protected Areas and GIS Laboratory (Brasil 2016).
We also used spatial eigenvectors as a proxy of other environmental and/or socioeconomic processes that have geographic patterns and were not explicitly included in this study.We used the Principal Coordinates of Neighbor Matrices (PCNM) method to generate the spatial eigenvectors (Griffith & Peres-Neto 2006).PCNMs were based on the distance among municipality centroids, with a truncation distance equal to 307.1 km.A total of 53 PCNMs were generated, however, only one (PCNM2) was used in further analyses.To select the spatial eigenvector, numerous linear regressions (combining all eigenvectors) with the response variable (productivity) were carried out and we selected an eigenvector that captured the spatial structure (indicated by lower spatial autocorrelation in residual regression).PCNMs were performed in SAM (Spatial Analysis in Macroecology v 4.0) (Rangel et al. 2010).
To obtain the ecological suitability of E. oleracea for each municipality we used the Ecological Niche Model approach.We compiled occurrence data for E. oleracea by searching the scientific publications available on the "Centro de Referência em Informação Ambiental" (CRIA 2013) website (http:// splink.cria.org.br) and other publications known to us (Fig. 1).A total of 95 unique occurrence points were obtained.
We used five climatic variables in the model: Annual Mean Temperature (BIO1), Temperature Annual Range (BIO7), Precipitation of Wettest Month (BIO13), Precipitation of Driest Month (BIO14), and Precipitation of Warmest Quarter (BIO18).The same variables were used for the future scenario, derived from a pessimist (A2a) scenario of the CCCma global climate model (Canadian Centre for Climate Modelling and Analysis) (see Karl & Trenberth 2005).The climate scenarios (current and future) were obtained from Worldclim (www.worldclim.org), with values projected to the year 2080 for the future scenario (Hijmans et al. 2005).All climate variables were converted into a grid with 0.0417 degrees of resolution.We selected these variables using the Jackknife criteria (considering all 19 climatic variables available in the Worldclim database) (see, for example, Nemésio et al. 2012).We generated potential species distribution models for South America using the MaxEnt algorithm in MaxEnt 3.3.3(Phillips et al. 2006).The input parameters followed the program default choices; except that the number of iterations was set at 1000 and duplicates were removed.We tested the efficiency of the model using the area under the ROC curve method (AUC; see Elith et al. 2006).

Data analysis
We used the Akaike Information Criterion (AIC) to select the best model (with set of predictors variables) that explains the productivity of E. oleracea in Brazilian municipalities.We used the following variables as predictors: Present Environmental Suitability (PES), GDP per capita (R$), Municipal HDI, %CONSERV, %RURAL, and PCNM.We generated 63 models, however, we presented only models the best models for productivity (delta AIC<2; Burnham & Anderson, 2002).The AIC was performed in the Spatial Analysis in Macroecology (SAM V.4) software (Rangel et al. 2010).The variables PES, GDP, and productivity were log-transformed (logX+1) while HDI, %CONSERV, and %RURAL were arcsine-square root transformed.

Results
The majority of the 200 municipalities that had 10 years' of productivity data for E. oleracea were concentrated in the northern and northeastern regions of Brazil.Therefore, the average productivity (i.e., 2002-2011) of E. oleracea is higher in these regions.A total of 92 municipalities of northern and northeastern of Brazil (36% of all investigated municipalities) have average productivity ranging between 0.9886 and 0.0102 ton.ha −1 .There were a total of 6 municipalities with the highest average productivity values (5.46 to 11.98 ton.ha -1 ); these municipalities were located in the state of Pará.
The state of Pará has 94 municipalities that currently produce E. oleracea (46% of the investigated municipalities), and these municipalities together are responsible for 86% of the total national productivity of E. oleracea.The municipality with the highest average yield was Limoeiro do Ajuru (Pará-PA), which recorded 11.9 ton.ha -1 , followed by the municipality of São Domingo do Capim (PA), with 10.0 ton.ha -1 .
The potential geographical distribution of E. oleracea indicated that under the present climate scenario, the species is distributed mainly in the northern and the northeastern part of Brazil (Fig. 2A).The area under the curve (AUC) generated by the MaxEnt algorithm was equal to 0.941, indicating an efficient model.Under the future climate change scenario, the limits of the distribution range underwent few changes; moreover, there was no size reduction in the geographic distribution (Fig. 2B).While, some regions that were considered poorly suited became appropriate; other regions with high climate suitability became unsuitable under the future climate scenario.
More than 60% of the municipalities that produce E. oleracea may gain favorable climatic conditions in the future, according to future climate change predictions (Fig. 3).The municipalities that are expected to gain suitability in the future are in the state of Amazonas (Tapauá, Jutaí, and Tefé), whose suitabilities will change from 0.33, 0.44 and 0.53 to 0.61, 0.66 and 0.74, respectively.The municipalities that are expected to reduce the greatest amount of suitability are in the state of Pará (Alenquer, Medicilândia, and Uruará), whose suitabilities will change from 0.53, 0.65 and 0.52 to 0.37, 0.52 and 0.39, respectively.
We generated 63 models to explain the average productivity of E. olereacea by municipalities.By using AIC, we found that the best model to explain the average yield was the model that combines Environmental Suitability and Spatial Filters, followed by the model using only spatial filters (Tab.1).It is noteworthy that the top five models were good because the AICcs were lower than two (Burnham & Anderson 2002).Moreover, the PES presented a positive relationship with the productivity of E. oleracea, which means that productivity tends to be higher in municipalities found in regions with higher environmental suitability (predicted by ENM) (Tab.2).Another important variable was the spatial eigenvector, however, in the present paper this variable indicates the spatial structure of productivity, thus geographically close municipalities tended to present similar productivity values (Tab.2).

Discussion
Here, we have project the potential geographical distribution of E. oleracea in Brazil under both current and future climate scenarios, and highlighted potential areas for its sustainable harvest.E. oleracea presented a geographic distribution in northern and northeastern of Brazil, covering only some states (see above).Thus, E. oleracea is highly dependent on the moisture that restricts its natural geographical distribution, since it requires wellirrigated soils, such as those found in naturally occurring floodplains (Muller 2016).
This study has demonstrated the relevance of ENM for predicting the productivity of E. oleracea, along with the potential impact of climate change on the geographic distribution of this species.Numerous studies have shown the importance of ENM as ecological predictors.In fact, most studies that use ENM have used these models to understand (sometimes visually) the geographic distribution of species under different climate scenarios (e.g.Collevatti et al. 2013).However, only recently have studies used models as predictors of population/genetic information (see Diniz-Filho et al. 2009b;2015;Wal et al. 2009;Nabout et al. 2011;2012;2016;Tôrres et al. 2012).In all cited papers, niche models were good to evaluate the positive relationship between habitat suitability and productivity or density.Although, various algorithms (i.e.techniques) exist for the construction of niche models (e.g.MaxEnt, GARP, distances approaches; see Franklin 2009), with variations in the relationship with the biology of studied species (see Tôrres et al. 2012, Diniz-Filho et al. 2015), generally all techniques of niche models have found positive Pearson correlation coefficients (see Diniz-Filho et al. 2015).Thus, different techniques of niche models were able to predict productivity or other demographic variables (see Diniz-Filho et al. 2015).Here, we used the MaxEnt algorithm, which is one of the most commonly used methods (Vaz et al. 2015), it is a machine and statistical learning method (Franklin 2009), which associates presence-only species records with bioclimatic variables, according to the principle of maximum entropy (see Elith et al. 2011).Despite the wide availability of methods, MaxEnt is among the top-performing methods in terms of prediction accuracy (Elith et al. 2006).
The positive relationship between suitability and productivity shows the relevance of ENM for the management of native species of economic importance.Thus, by using maps generated by niche models we can indicate current and future geographic regions for sustainable harvesting and suggest areas for the conservation of species of economic importance.Therefore, the integration of science and policy is essential to make decisions, such as the selection of areas for conservation and the planning of methods to stimulate rural economies.Particularly for the Amazon, which has undergone rapid changes to its landscape (see Sawyer 2008), where niche models might be an important tool to understand the distribution of species in a still little known biome.
The other predictors used in this study (%CONSERV, %RURAL, GNP and HDI) were not important in explaining the variation of E. oleracea productivity between the studied municipalities.In fact, other studies have shown that socioeconomic variables are not important in explaining the productivity of mangaba, but are important in explaining the productivity value (see Nabout et al. 2016).Thus, municipalities with higher GPD tend to have higher mangaba productivity value (Nabout et al. 2016).
Most plant ENMs have shown that there will be a reduction in the potential geographical distribution of such species under current predictions of future climate scenarios (e.g.Diniz-Filho et al. 2009a;Mendoza-Gonzalez et al. 2013;Simon et al. 2013).However, E. oleracea has showed a potential increase in its geographic distribution; although relatively small.Moscoso et al. (2013) conducted a study with 12 native Amazon species; six woody species and six species of palm (E.oleracea included) under current climatic conditions, and the results of the generated models were similar to those in the present study.Of course, this is good for the preservation of biodiversity, and for the local economy as it is a species of economic interest.Therefore, it is important to use this approach, especially to stimulate other municipalities that will have areas of potential climate suitability to invest in long-term conservation planning.
Although our study suggests there will potentially be new occurrences of this species, it is important to consider that this scenario is optimistic because other variables, such as land use, change of pollinators and pests were not considered in the niche model (Hannah et al. 2002).In other words, although climatically suitable according to ENM, new regions may not be appropriated for the presence of E. oleracea at the present moment.In fact, the future distribution model for E. oleracea indicates regions that are currently known as the arc of deforestation (see Malhi et al. 2008), especially in the state of Tocantins.
Finally, based on the results presented in this paper, we indicate some guidelines for future studies: i) assess whether current protected areas include the future distribution of E. oleracea, particularly for regions in the arc of deforestation; ii) assess the economic loss and gain of this species in the municipalities that currently produce it; iii) investigate appropriate management strategies, preventing its overexploitation; and, iv) encourage new municipalities to produce this species, especially those that are currently in regions where the future climatic conditions will be suitable.
2036/2013).The work by JCN has been supported by CNPq grant (306719/2013-4).ULV was supported by CAPES scholarships.We thanks Levi Carina Terribile and anonymous reviewer for suggestions that improved previous version of this manuscript.

Figure 1 .
Figure 1.The locations of the 95 occurrence data-points for Euterpe oleracea used in the analyses.

Figure 2 .
Figure 2. Environmental suitability map for Euterpe oleracea, with its occurrence in the present (A) and future (B) scenario.Please see the PDF version for color reference.

Figure 3 .
Figure 3. Frequency distribution of Δ values of climate suitability of Euterpe oleracea -producing municipalities.

Table 1 .
Top five models generated and AIC numbers found for the combinations of variables: Present Environmental Suitability (PES), Human Development Index (HDI), Proportion of Rural Population (%RURAL), Percentage of Conservation Unit Areas (%CONSERV), and Spatial eigenvectors filters (PCNM)

Table 2 .
Importance of each variable and the angular coefficient indicating the relationship with the dependent variable (i.e., productivity of Euterpe oleracea).The importance of considering the AIC across all models in the set in which this variable occurs.AIC ranges from 0 to 1, with values closer to 1 indicating that the variable occurred in the best models.