Current presence of species

A total of 542 georeferenced presence points were obtained for S. purpurea and S. mombin from three main sources: National Forest and Soil Inventory (^{CONAFOR, 2009}) with 23 % of the collected points, national and international herbarium collections reported by ^{Cruz-León et al. (2010)} with 42 % and an international biodiversity information server which contributed 35 % (^{Global Biodiversity Information Facility [GBIF], 2015}). Data were homogenized by eliminating points found outside the national territory, which did not contain collection information or that would have been reported before 1980, in order to avoid possible georeferencing or identification errors. Those points located in urban use sites were also eliminated as the aim was to model the natural distribution potential of the species (^{Soberón & Nakamura, 2009}). As a result of the debugging, 316 points for S. mombin and 189 for S. purpurea were used. It was observed that Yucatán, Chiapas and Campeche accounted for 44 % of the references used in the modeling.

Biogeographic variables

Through a literature review, the environmental variables reported to have the most influence on plum distribution were selected. Twenty-two predictor variables were considered (Table 1): 16 bioclimatic, two edaphic, one latitudinal, one geological and two topographic. Edaphic, geological and topographic information was obtained from official information sources in Mexico (^{Instituto Nacional de Estadística, Geografía e Informática [INEGI], 2005}, ^2013a, ^2013b). On the other hand, detailed mapping of the bioclimatic variables was generated using the methodology described by ^{O’Donnell and Ignizio (2012)}. These variables have the advantage of being biologically significant by combining temperature and precipitation throughout the year, allowing analysis of seasonal patterns (^{Hijmans, Cameron, Parra, Jones, & Jarvis, 2005}).

To obtain the bioclimatic variables of the reference scenario, the information of ^{Gómez et al. (2008)}, which is based on the concept of areas of climatic influence, defined as homogeneous in annual precipitation and temperature, was used. Climate change models were obtained and processed from the UNIATMOS platform (^{Fernández, Zavala, & Romero, 2009}). Bioclimatic variables with climate change scenarios were processed with the same methodology of ^{O’Donnell and Ignizio (2012)}. The HADGEM2_ES (English), GFDL_CM3 (American) and Ensamble REA models were considered. These models were selected due to their good ability to reproduce the observed climate, and the fact they have been widely used in Mexico (^{Conde, Estrada, & Martínez, 2011}). It has been reported that on average for Mexico the HADGEM model projects increased temperature and decreased precipitation, while the GFDL forecasts a general increase in both variables (^{Gómez-Díaz et al., 2007}).

Within the climate scenario groups, representative concentration pathway (RCP) 8.5 W∙m^-2 was used. An RCP measures the amount of energy added to the atmospheric system and that generates effective heating (^{van Vuuren et al., 2011}). The potential distribution of plum was projected to the horizon 2075-2099 to know the possible long-term impacts. The combination of the models and radiative forcing resulted in three climate change scenarios for both species of plum: 1) HADGEM RCP 8.5 W∙m^-2, 2) GFDL RCP 8.5 W∙m^-2 and 3) REA RCP 8.5 W∙m^-2. Each represents an alternative of how the climate could change in the future (^{Conde et al., 2011}).

Ecological niche modeling

Potential distribution was estimated using the maximum entropy algorithm of Maxent software, which has been widely used in Mexico with good results for modeling the distribution of flora and fauna (^{Miller & Knouft, 2006}; ^{Trejo et al., 2011}). Entropy, in this context, is defined as the measure of the amount of information considered in the selection of an event (^{Phillips, Anderson, & Schapire, 2006}). The Maxent algorithm analyzes the conditions where the species has been reported by contrasting the georeferenced presence points with environmental variables. With this, a set of conditions under which the species can live ideally, what is known as an ecological niche, is defined (^{Soberón & Nakamura, 2009}). By regression functions, the software estimates the probability of species presence at a site, which can be geographically translated into a potential distribution (^{Elith et al., 2011}).

Unlike other species distribution models, Maxent has the advantage of using presence-only species records, which facilitates the use of the information available in the databases consulted. The output format selected in Maxent was the logistic one, as the aim was to determine the probability of presence under the typical conditions of the species (^{Elith et al., 2011}). Jacknife resampling was implemented to measure the degree of importance of the environmental variables in the potential distribution. Also, response curves showing the probability of plum presence in light of the variation in the different environmental factors were obtained.

Validation of the distribution models

For statistical validation, 30 % of the total presence points were randomly extracted, retaining the remaining 70 % for modeling, as proposed by ^{Naoki, Gómez, López, Meneses, and Vargas (2006)}. The Maxent software provides an internal evaluation based on the area under the curve (AUC) parameter; however, four independent indices were used to assess the overall performance of the current species distribution model.

1) The sensitivity index is defined as the proportion of correctly predicted presences, where a high value indicates low error of omission, or type I error. 2) The specificity index is designed to determine the proportion of correctly predicted absences, where a high value indicates low error of commission, or type II error. 3) The positive predictive power (PPP) is the proportion of correctly predicted presences in relation to all locations. 4) The Kappa index is a global index that measures the degree of agreement between two estimators, where a high value indicates substantial agreement (^{Parra, Graham, & Freile, 2004)}).

All indices required obtaining a confusion matrix, which was generated and classified in accordance with the methodology of ^{Landis and Koch (1977)}. A validation was also performed by geographical overlapping of the modeling results with the S. purpurea producing regions identified by ^{Avitia-García et al. (2000)} and with the potential distribution of S. mombin proposed by ^{Pennington and Sarukhán (2005)}.

Mapping current and future potential suitability

Four suitability classes were defined as a function of the probability of species occurrence in the territory, as proposed by ^{Monterroso-Rivas, Gómez-Díaz, and Tinoco-Rueda (2013)}. The 0-25 % probability range was called“unsuitable,” 25-50 %“marginal,” 50-75 %“suitable” and 75-100 % “optimal.” For purposes of this study, suitability areas were defined as those having between 50 and 100 % probability of presence. These thresholds base their division on the arithmetic concept of quartiles and have been used in studies where the categorization of values is required. By using a geographic information system, information resulting from the Maxent software was processed and classified. The maps generated were intersected with Series V land-use and vegetation mapping (^{INEGI, 2013a}).

Finally, every time the current scenarios for the species were obtained, 16 variables were replaced by those obtained under the aforementioned climate change models. The software was re-run and the new suitability areas were obtained and compared with the reference scenario.

Spondias purpurea

Mexican plum would find suitable sites for growing in 3.8 % of the national territory, according to the reference scenario, equivalent to about 7,400,00 ha. Of that area, 11 % is in the optimal suitability category and 89 % is in the suitable one. The average January temperature, the minimum temperature of the coldest period and thermal seasonality were the main variables that helped explain 63 % of the distribution of the species in the country. By also considering the kind of rock, elevation and soil unit, 81 % of the species’ current presence can be explained. The importance of temperature for the species coincides with the report by ^{Avitia-García et al. (2000}), who state that cold is harmful to the species and that it only tolerates low temperatures for short periods of time, but they do not specify limits and time. In this sense, the high contribution of thermal seasonality in the modeling also has a reference in the literature, as it has been established that the Mexican plum requires that the difference between the hottest and coldest months does not exceed 10 °C (Avitia- García et al., 2000).

The presence of S. purpurea was found to be strongly related to climate group Aw, and its distribution proved to be predominant on the Pacific Ocean slope, in central Veracruz and the Yucatán Peninsula (Figure 1). ^{Cruz-León et al. (2010}) also located it primarily on the West Coast. Comparing the results obtained in this study with the distribution of plum producing areas in Mexico reported by ^{Avitia-García et al. (2000)}, most coincide with suitable and optimal suitability categories found by the modeling.

Figure 1 Current potential distribution for Spondias spp. in Mexico, a) S. purpurea and b) S. mombin. 

Of the nine producing regions reported in the literature, five overlap in optimal areas (Nayarit, Jalisco-Colima, Oaxaca-Chiapas, Veracruz and Yucatán), two in suitable and marginal areas (Morelos and Guerrero) and two did not overlap (Sinaloa-Sonora and Tabasco). Regions reported as producers that do not match the present study were the border region between Sonora and Sinaloa and central Tabasco. The existence of S. purpurea in very wet or very dry climates is closely related to anthropogenic effects, as ^{Cruz-León et al. (2010}) have established that the presence of this species in very wet climates can only be due to human action.

In total, 53 % of the suitability areas were located in primary and secondary vegetation of low-height deciduous forest, in secondary arboreal vegetation of medium-height sub-deciduous forest and in rainfed agriculture and cultivated pastures. This is consistent with several studies that associate the presence of S. purpurea with low-height deciduous forest (^{Ruenes-Morales, Casas, Jiménez-Osornio, & Caballero, 2010}; ^{Rzedowski, 2006}). On the other hand, ^{Avitia-García et al. (2000)} describe it as a species that forms part of the secondary vegetation derived from tall or medium-height sub-deciduous or deciduous forests. Cultivated pasture has been one of the major changes that humans have made in low-height deciduous forest areas for livestock activities according to ^{Castelán, Ruiz, Linares, Pérez, and Tamariz (2007)}, so a reference was found for the important presence of this species in this land-use class.

Land-use classes solely related to human activity make up 31 % of all suitability areas. However, natural vegetation classes associated with the presence of the species account for 40.9 % of suitability areas (10 % more than those related to human activities). This suggests that despite human intervention it has not been fully domesticated and still prefers its wild niche. This is consistent with the observations of ^{Hernández-Bermejo and León (1992)}, who highlighted the semi-domesticated character of the species.

Regarding the results obtained with climate change, the main areas of change in terms of suitability would be recorded in Sinaloa, Nayarit, northern Veracruz, and the northern strip of the state of Puebla. According to the spatial distribution reported under climate change models at 8.5 W∙m^-2, the area of S. purpurea would decrease by 10.3 % and 5.1 % in the optimal and suitable classes respectively, considering the average of the three models. Adding up the suitable and optimal areas, an overall 5.7 % decrease in the suitability area for the species is expected. Figure 2 shows the potential distribution of both species with climate change scenarios.

Figure 2 Potential distribution with climate change for Spondias spp. at a distant horizon under RCP 8.5 W∙m^-2: 1) S. purpurea and 2) S. mombin. Scenarios: A) GFDL_CM3, B) HADGEM2_ES and C) Ensamble REA. 

All models agreed in predicting a decrease in optimal areas. The most optimistic was GFDL in finding a 7.7 % decrease, the HADGEM was the most pessimistic by considering a 14.6 % drop and the REA was at an intermediate point by finding a 8.5 % decrease. The climate change scenarios showed agreement in predicting an increase in marginal areas, with 5.7, 3.4 and 8.1 % for the REA, HADGEM and GFDL models, respectively. The increase in this category is because suitable areas would cease to be so and would be located in lower suitability categories, according to the climate change scenarios.

Spondias mombin

According to the baseline scenario, S. mombin could be potentially distributed in 6.6 % of the national territory, equivalent to 12,900,000 ha. However, of that area, only 1 % is in the optimal category, while the rest is in the suitable one. Although the distribution of a species is associated with environmental features and ecological interactions, collecting a larger number of specimens at certain sites may be influenced by human interest. In this regard, ^{Ruenes-Morales et al. (2010}) have documented the presence of the genus Spondias in the Yucatán Peninsula, which could have led to an oversampling in that area.

The variables that contributed most to modeling the presence of the species were the minimum temperature of the coldest period, average monthly rainfall in April and annual thermal oscillation, contributing 66 % to the model. If the influence of the tropics and the type of rock is added, 85 % of the distribution of the species is explained. In this regard, ^{von Carlowitz et al. (1991)} state that the average temperature should be 25.5 °C and the minimum 16 °C. On the other hand, ^{Orwa, Mutua, Kindt, Jamnadass, and Anthony (2009)} specify that S. mombin is severely damaged by low temperatures, and although the limit is not specified, it matches this study. The distribution of the species was influenced by Am-group climates, so an increased probability of presence was recorded in the humid areas of Veracruz, Tabasco and Chiapas, in addition to covering almost the entire Yucatán Peninsula. For their part, Quintana Roo, Campeche and Yucatán account for 65 % of the areas suitable for S. mombin. If Chiapas is added to the above states, 81 % of the total national area is reached. It has been documented that S. mombin is distributed in warm humid areas (^{Pennington & Sarukhán, 1998}).

The three land-use types that account for most of the suitable area for the production of Spondias spp. are the secondary tree vegetation of the medium-height sub-deciduous and deciduous forest and cultivated pasture, adding up to 52 % of the national total. These results are consistent with studies reported for the species, where it has been documented that it is found associated with secondary vegetation of tall or medium-height evergreen and sub-evergreen forests and of medium-height sub-deciduous forests (^{Avitia-García et al., 2000}). ^{Cruz-León et al. (2010}) describe it as common in secondary forests and present in primary forests.

The climate change scenarios project an average overall decrease of 13.3 % in suitable areas for S. mombin. Of the three models, the GFDL was the most pessimistic by predicting a 16.9 % decrease, the REA was more optimistic by forecasting an 8 % drop and the HADGEM model projected a 14.8 % decline. Considering climate change scenarios (Figure 2), it is expected that the main areas of change, according to the GFDL model, relative to the baseline scenario, will be Yucatán, Quintana Roo, Campeche and Veracruz, whereas the HADGEM model predicts changes in Yucatán, Quintana Roo, Veracruz and Chiapas and the REA model reports potential changes basically in the Yucatán peninsula.

Figure 3 shows the change trends in the suitability categories by species. S. mombin would increase its suitability areas by 65.8 %, which is the average of the three models. The suitability areas tend to increase under climate change scenarios due to the decreased susceptibility of Spondias spp. to low temperatures, a limiting factor which has been widely reported by various studies (^{Avitia-García et al., 2000}; ^{Cruz-León et al., 2010}; ^{von Carlowitz et al., 1991}). However, this increase should be viewed with caution since it represents a very limited area (135,000 ha).

Figure 3 Changes in suitability categories according to climate scenarios, a) S. mombin and b) S. purpurea. 

Considering only the optimal areas, the REA model is the most optimistic by predicting an increase of 107.4 %, while the GFDL projects an increase of 40.7 %. On the other hand, areas with the suitable category provide a differentiated decrease of 17.8, 15.9 and 9.9 %, according to the GFDL, HADGEM and REA models, respectively. The average change decrease in this category was 14.5 %, according to the projections by the three models. The area defined as marginal would increase by an average of 3.8 %, with the GFDL predicting the greatest change (5.8 %). Although the variation in unsuitable areas is projected to average 0.7 % in terms of area, there would be 1,131,000 ha with very poor conditions for the development of S. mombin.

Final comments

According to the climate change scenarios, the GFDL model proved to be the most severe in terms of suitability by projecting a 12.2 % decrease in the suitable area, considering both species. The HADGEM is at an intermediate point by predicting a 9.5 % decrement in suitability areas. The REA model was the most optimistic in finding a 6.8 % drop. Areas with some degree of suitability (sum of optimal and suitable categories) would decrease on average by 5.7 and 13.3 % for S. purpurea and mombin, respectively, which means the latter species is the more vulnerable of the two studied. Figure 4 provides a comparison of the suitability classes according to each climate change scenario assessed.

Figure 4 Percentage changes in area with suitability for Spondias spp. under climate change scenarios and RCPs. 

Further research is needed to broaden our understanding of the impact of the main factors affecting the distribution of these species. This, in turn, will enable identifying, proposing and evaluating measures aimed at adapting to or mitigating climate change within the context of rescuing plant genetic resources native to Mexico.