Study area

The study area includes the southern United States of America, Mexico and Central America (Figure 1), a region that covers 3 946 177.4 km² and is culturally known as the Mesoamérica and Aridoamérica regions (^{Kirchhoff, 1943}).

Figure 1 Study area for Phaseolus lunatus. 

Occurrence data

A matrix of occurrence data of Phaseolus lunatus was integrated from georeferenced passport data based on information from the following entities: International Center for Tropical Agriculture (CIAT), National Institute of Ecology, University of Guadalajara, National Research Institute Forestales, Agrícolas y Pecuarias (INIFAP), ^{Fraitag and Debouck (2002}), as well as private data provided personally by Dr. Rogelio Lépiz Ildefonso, Professor-Researcher, CUCBA, University of Guadalajara.

Climatic variables

Normal monthly images of maximum temperature (12), minimum temperature (12), average temperature (12) and precipitation (12) were used, in addition to images of the 19 bioclimatic variables, all derived from the Worldclim system and corresponding to the period 1960-2000. (^{Hijmans et al., 2005}). In addition, 12 normal monthly images of annual potential evapotranspiration (ETP) were used, which were obtained from the CGIAR global geospatial soil water balance database (Trabuco and Zomer, 2010). The longitude, latitude and altitude of the presence sites were also added as analysis variables. For this, the altitude image of the USGS/NASA SRTM site was taken (^{Jarvis, 2008}). With the monthly ETP and precipitation data, 12 monthly humidity index images (IH) were generated using the IH = P/ETP expression. All this information was worked with a resolution given by a cell size of 30 arc seconds. In total, 94 variables were analyzed.

Climatic characterization of presence sites of Phaseolus lunatus

The geographic coordinate information of the lima bean occurrence sites and the climatic and bioclimatic images of the Worldclim system were used to characterize the climatology of the distribution sites of this species. For this, the information extraction procedure of the ArcGis program was used (^{ESRI, 2010}). The extracted information was incorporated into the P. lunatus occurrence data matrix in the Microsoft Excel spreadsheet. Finally, the climatic type of the presence sites was characterized using the system proposed by ^{Ruiz et al. (2018}) and that illustratively refer ^{Sánchez et al. (2018}), in which 36 environments are theoretically possible, product of the combination of six hydric zones, determined by intervals of a humidity index (IH= annual precipitation/annual potential evapotranspiration) and annual average temperature intervals (Tma), respectively (Table 1).

Statistical analysis

From the climatically characterized site matrix of occurrence, a multiple correlation analysis was performed to refine the list of climatic variables. The correlation coefficient criterion ≥ 0.95 was used for the elimination of variables. With the refined list of variables, a principal component analysis was carried out and the variables with the highest statistical weight in the presence of lima beans were identified. The NTsys 3.3.0 program was used for this analysis (^{Rohlf, 2009}).

Climate descriptors

From the climatic variables selected for the analysis and considering the climatic characterization of the occurrence data matrix, the ecological descriptors of lima beans were deduced, for which the minimum, maximum and average values of each variable were considered.

Modeling of potential niches of distribution of Phaseolus lunatus

Prior to modeling potential distribution niches, a randomness test was applied to determine if the sites of occurrence of Phaseolus lunatus are randomly distributed (^{Bivand et al., 2014}). This test was run using a program implemented in R. The Maxent model (^{Phillips et al., 2006}) was used to model the potential distribution of lima beans under the current or reference climate scenario and under the climate scenarios year 2050 and year 2070.

For the modeling of the potential distribution of P. lunatus in future or climate change scenarios, the data derived from the general circulation models were used: HadGEM2ES (European Network for Earth System-MetOffice Hadley Center) and CCSM4 Community Earth System Model (CCSM4) under two routes of concentration of greenhouse gases (RCP) 4.5 and 8.5. The first model, developed in the United Kingdom, incorporates elements that allow for adequate modeling of continental vegetation as it includes various types of vegetation: broadleaf trees, needle trees, C₃ and C₄ grasses, as well as shrubs (^{Bellouin et al., 2007}).

The second model is developed in the United States of America by the National Center for Atmospheric Research (NACAR), which incorporates vegetation cover as one of its four coupled components. These models have been statistically evaluated by ^{Conde et al. (2011}), conclude that the CCSM Model is one of those that best describes the climate of the Mexico and Central America region, while HadGEM has already been used for studies of the vegetation of Mexico with good results (^{Instituto de Ecología, 2009}).

The Maxent model has been described as being especially efficient in handling complex interactions between response variables and predictor variables, as it is a robust model even in small data samples (^{Fourcade et al., 2014}). Maxent uses the principle of maximum entropy on species presence data to estimate a group of functions that relate environmental variables with habitat suitability to approximate the potential geographical distribution of the species (^{Phillips et al., 2006}). Maxent represents the distribution of a species; through a probability function P on a set X of sites in the study area. The probability function for i environmental variables is expressed as follows.

𝑃(𝑥)=𝑒𝜆∙𝑓(𝑥)/𝑍𝜆. Where: P(x)= probability function; λ= vector of weighting coefficients; f= corresponding vector of functions of environmental variables; Z= normalization constant to ensure that P(x)= 1. The P(x) values thus obtained represent relative suitability values for the presence of the species, thus constituting the basis of a potential distribution model.

In the modeling of the potential distribution of Phaseolus lunatus with Maxent, the default parameters were used, using the replication execution crossvalidate with 1 000 repetitions. 50% of the data was used for training and 50% for testing. In climate change scenarios, this configuration was maintained. Finally, a presence/absence binary map of Phaseolus lunatus was generated using the cutoff threshold method, the one with an omission rate close to 0% but with the highest cutoff value (^{Sánchez et al., 2018}).

Climatic adaptation of Phaseolus lunatus

In Figure 2, the distribution of Phaseolus lunatus in the different agroclimatic zones of the study area can be observed. As can be seen on the map in Figure 2 and Table 2, this species is present in 17 of the 26 agroclimatic zones in the study area. However, of the 17 agroclimatic zones in which P. lunatus occurs, two correspond to arid zones.

Arid areas can hardly accommodate the development cycle of a legume such as P. lunatus, since the annual precipitation levels in these environments is much less than 500 mm, a value indicated as the minimum rainfall threshold for adaptation of this species (^{Baudoin, 1989}). The presence of P. lunatus in arid zones in the study area is most likely due to having access to the water resource, but not through rainfall.

Figure 2 Agroclimatic regions and distribution of Phaseolus lunatus. 

Table 2 shows that the greatest number of sites of occurrence of Phaseolus lunatus occurs in the warm humid-sub-humid agroclimatic zone, indicating that this is the climatic environment that is closest to the optimal climatic conditions for this plant.

Without considering the arid hydric zone, P. lunatus is able to distribute itself in four hydric zones (from semi-arid to humid) and four thermal zones (from tempered to very warm) although it seems to prefer the sub-humid-humid to humid interval and the interval of semi-warm to very warm. The distribution of this legume in semi-arid and tempered conditions is definitely very small (Table 2).

Ecological descriptors

Table 3 shows the ecological descriptors of Phaseolus lunatus from the most representative climatic variables. As can be seen, even though this legume has a preference for semi-warm, warm and very hot environments (Table 2), its presence at altitudes of 0 to 2 386 m (Table 3) indicates a wide adaptability, since it can be distributed thermally from areas tempered (minimum value of 13.2 °C annual mean temperature, Table 3) even in tempered environments (Table 2). The altitude range found coincides approximately with that reported by ^{López et al. (2005}) from 0 to 2059 m and differs slightly from that indicated by ^{Baudoin (1989}) from 800 to 2 000 masl, as well as from that reported by ^{Debouck (2020}) from 50-2 800 m. These differences are surely due to the different geographical areas considered by these authors.

According to Table 3, the growth season interval in which Phaseolus lunatus occurs (83 to 365 days) coincides with the development cycle that ^{Debouck (2020}) reports from 90 to 365 days for P. lunatus.

This interval of 83 to 365 days of growing season, allows us to deduce that this plant adapts to environments with a short rainy season as well as to rainy environments practically all year round, with accumulated precipitation in the growth period of 250 to 4 250 mm, so that some populations of P. lunatus present in these extreme humidity conditions could be indicating the existence of populations carrying genes resistant to these conditions, as suggested by ^{Maxted et al. (2013}).

Similarly, in the case of temperature, P. lunatus adapts to growing seasons with an average minimum temperature of 8 °C as well as growing seasons where the maximum temperature can average even 37 °C. The populations representative of these thermal values may also be of interest due to their possible adaptation genetics to extreme temperatures. Furthermore, according to Table 3, there are populations of Phaseolus lunatus that experience the warmest month with maximum temperatures above 40 °C and others that have the coldest month with temperatures below 0 °C, especially those forms of the species of long biological cycle, of indeterminate growth habit (^{Debouck, 2020}).

It is interesting to note that the reports on P. lunatus indicate a thermal range of 16 to 26 °C (^{Debouck, 2020}), in the present study the interval obtained is broader, 13.2 to 29.2 °C on an annual scale and 13.2 to 29.9 °C at the level of the growing season (Table 3). These values constitute thermal adaptation limits not reported to date for P. lunatus, which could indicate that this species has had the ability to explore more thermally extreme areas over time.

Likewise, the present investigation obtained an upper limit of annual accumulated precipitation of 4 250 mm, an amount that exceeds previous reports such as that of ^{Debouck (2020}) of 2 800 mm. Finally, according to Table 3, the average humidity index during the production cycle indicates as lower and upper limits of adaptation of P. lunatus, 0.51 and 3.09, respectively. There is no previous reference to these descriptors, because variables such as evapotranspiration, moisture balance and humidity index are not typically parameters considered in the description of crop ecology.

Modeling the potential distribution of Phaseolus lunatus

The randomness test reported that the points of occurrence of Phaseolus lunatus have an agglomerated distribution and not a random distribution, therefore the distances between the points of occurrence were estimated and the data were re-sampled using the ArcGis program (ESRI, 2010), using a mesh constructed from geographic distances to obtain randomly distributed samples. With these samples, the modeling of potential distribution areas of P. lunatus was performed.

Of the total of 94 variables that were considered in the present study, the multiple correlation analysis and the principal components analysis allowed reducing the list to 29 variables with the greatest significance in the distribution of Phaseolus lunatus. Table 4 shows the Eigenvector value of the 29 climatic variables with the greatest significance in the distribution of P. lunatus in the first three main components.

It is observed then that in the first component variables that have to do with the humidity factor (precipitation) are grouped, in the second component, variables that have to do with the energy factor (temperature) and in the third component, other types of variables that are derived from the interaction of both factors.

When running the Maxent analysis with the 29 variables, a model was obtained whose operational curve or ROC curve reported an area under the curve (AUC) of 0.94, which according to ^{Araujo et al. (2005}) represents an excellent fit of the model and therefore allows a good discrimination of suitable areas against areas not suitable for Phaseolus lunatus (^{Thuiller et al., 2009}) when preparing a binomial map. The Jackknife analysis allowed identifying the variables with the greatest contribution to the modeling of the Phaseolus lunatus distribution. Table 5 shows the contribution of the 11 variables that contribute the most in this regard. Both Bio04, Prec09 and Bio16 remained as the three variables with the greatest contribution to the distribution of this legume in all the analyzed climatologies, both in the reference one and in the climate change climatologies.

The analysis of the evolution that these three variables would have in the climatologies of climate change is presented in Figure 3. There it can be seen that Bio04 tends to increase in climate change scenarios. Prec09, on the other hand, tends to decrease, with a percentage of -1% in the 2050 RCP4.5 CCMS4 model, to 15% in the 2070 RCP8.5 model HADGEM-ES.

Regarding Bio16, the CCMS4 model estimates a reduction both from 2050 to 2070 and when changing from RCP 4.5 to 8.5, while the HADGEM-ES model estimates a stable value of Bio16.

Figure 3 Expected evolution of three variables in various climate change scenarios. 

The binary map for the reference scenario revealed that, on a large scale, the current distribution of Phaseolus lunatus coincides with the potential distribution simulated by the Maxent model (Figure 4). However, on a smaller scale potential distribution area are visible that lack sites of current occurrence of this legume. This indicates the presence of suitable sites for future exploration and possible collection of new populations of P. lunatus.

Figure 4 Current distribution (green color) and potential (orange color) of Phaseolus lunatus in various climatologies: A) 1961-2010; B) HADGEM-ES RCP4.5 2050; C) HADGEM-ES RCP4.5 2070; D) HADGEM-ES RCP8.5 2050; E) HADGEM-ES RCP8.5 2070; F) CCMS4 RCP4.5 2050; G) CCMS4 RCP4.5 2070; H) CCMS4 RCP8.5 2050; I) CCMS4 RCP8.5 2070. 

In climate change climatologies, the HADGEM-ES model in the two RCPs indicates a decrease in the potential distribution surface of P. lunatus mainly in the Yucatan Peninsula and Central America (Figures 4B, 4C, 4D and 4E) that goes from -8 to -10% in RCP4.5 at a range of -13 to -37% (Figure 5) in RCP8.5. These effects on the Yucatan Peninsula are particularly important, since this area represents one of the most important centers of diversity of the Lima bean MII genetic group.

From the 2050 RCP4.5 scenario (Figure 4B), a decrease in potential areas for this plant is observed, in a direction from south to north of the Yucatan Peninsula, this is northern Guatemala and Campeche, which is accentuated in the following scenarios (Figures 4C, 4D and 4E) covering an increasing area of this territory, until the potential areas are reduced to a small area in the northeast corner of the Peninsula, in the state of Yucatan (Figure 4E).

In the case of the CCMS4 model, Figures 4F, 4G, 4H and 4I show the sequence of maps of potential distribution areas of the lima bean. In it, it can be seen that RCP4.5 represents favorable environmental conditions, which would cause an increase of 7% in 2050 and + 2% in 2070 (Figure 5), occurring mainly in the north-central and western region of Mexico, in addition to the Caribbean (Figures 4F and 4G), which could in the future give rise to the presence of new populations of lima beans in these regions and with it a new configuration of the distribution of this species in the study area. While RCP8.5 would lead to a reduction in the potential surface of lima beans in the Yucatan Peninsula and Central America (Figure 4) from -0.17% in 2050 to -22% in 2070 (Figure 5). These positive and negative results for P. lunatus are consistent with that reported for Tithonia diversifolia in Mexico by ^{Duran et al. (2020}).

Figure 5 Effect of climate change on the potential distribution surface of P. lunatus.