Loss and Gain in Potential Distribution of Threatened Wild Cotton Gossypium thurberi in Mexico under Future Climate

Abstract

Gossypium thurberi, a threatened wild cotton species native to northern Mexico and southwestern USA, is globally important because its agronomic traits can be introgressed into cultivated species to improve fiber quality and resistance to biotic and abiotic stressors. However, studies on the current and future potential distribution of the species are scarce. The objectives of this study were (1) to develop a distribution model of G thurberi using a Geographic Information System platform, (2) determine environmental factors that influence the current distribution of the species in Mexico, and (3) estimate the potential distribution of the species under current and future climates. We analyzed the following variables: Annual Available Soil Water (mm year⁻¹, AASW), Flowering Growing Degree Days (FGDD), absolute minimum temperature (°C, Tminabs), and altitude (amsl, ALT). Results showed that the current potential distribution of G. thurberi in northern Mexico, estimated at 112,727 square kilometers, is projected to be drastically reduced by 77 and 86%, considering a possible increase in temperature of 1.5 °C and 2 °C in near-future (2021–2040) and mid-future (2041–2060) climates, respectively, and a 100 mm reduction in average annual precipitation under both climates. The greatest reduction will be in areas in Sonora (Mexico) adjoining Arizona (USA), where the largest populations of the species are currently reported. AASW, FGDD, and ALT jointly influence the distribution of G. thurberi, with AASW as the dominant factor under future climate change. The areas that may continue to harbor populations of G. thurberi under future climate will present AASW of 0.2–55.6 mm year⁻¹, FGDD of 242–547, and ALT between 550 and 1561 amsl. The projected future potential distribution in the country includes new suitable areas, including one in the Trans-Mexican Volcanic Belt, that may serve as refuge areas. The findings can contribute to the design of more precise collection efforts and conservation strategies to prevent species extinction.

1. Introduction

Cotton is the world’s most important textile fiber crop and also a significant source of edible oil and animal feed [1,2]. Grown in over 90 countries, it provides income to millions of farmers, many in developing countries [3]. However, with the use of only a few related genotypes of cotton (only four out of 53 species are widely cultivated) [4,5] and the adoption of transgenic cultivars [6,7], the genetic base of the crop has narrowed, making the crop vulnerable to biotic and abiotic stressors. Wild Gossypium species can be exploited to infuse novel genetic variations into the gene pool of cultivated cotton to protect it against agro-environmental challenges [4].

Gossypium thurberi Todaro is a wild cotton species native to Mexico and parts of southwestern USA [8]. It is the only cotton species naturally occurring in the United States, found in Arizona below 34° latitude [9]. In Mexico, it is distributed in the northern and northeastern regions of Sonora and in small patches of western Chihuahua [9,10,11]. The species is globally important because its agronomic traits can be introgressed into cultivated species to improve fiber quality (fineness, strength and elongation) [5,12], resistance to diseases, such as Fusarium and Verticillium wilt [13,14], and pests like silverleaf whitefly, Bemisia argentifolii (Hem., Aleyrodidae) [8], and resistance to frost, hard freezes, heat, and drought [15,16,17,18]. Among its other breeding traits are prolific boll bearing, high ginning out turn, and better spinning [19]. It has also been used as breeding material to reduce calyx dust causing black lung disease [20]. The successful use of G. thurberi in cotton crop improvement in various countries, such as the USA [12], China [14], India [21], Uzbekistan [13], and Bulgaria [22], has resulted in new high-yielding, early-maturing, and disease-resistant commercial varieties [14].

Global and regional studies have shown that climate change can cause geographical redistribution of species, thereby changing ecosystems [3,23,24,25]. Collection trips in Arizona in 2015 and 2017 by Frelichowski et al. [26] showed less G. thurberi populations than expected. They were found in rangelands where grazing by livestock poses a risk, but were absent in some areas where their presence had been reported (e.g., Prescott and Chiricahua Mountains), indicating possible losses of populations. In Mexico, the loss of important genetic diversity has been attributed to climate change and other factors, such as the use of pesticides, genetically modified crops, and other advances in modern agriculture as well as deforestation, touristic developments, and urbanization [27,28,29]. Several wild cotton species are now critically in danger of extinction [28,29]. G. thurberi is among those in the threatened category [30].

To prevent the devastating loss of valuable genetic resources, extinction risk assessments and prioritized conservation strategies must be established [29,31,32]. Species distribution models (SDMS) are commonly used in conservation biogeography to estimate environmental suitability for a species in geographic space by relating species presence or abundance to environmental variables [33,34]. Species distribution modeling has come to include ecological or environmental niche modeling, habitat suitability modeling, and climate envelope modeling, although they differ in their approaches [35]. The integration of statistics, biology, and geography in SDMS has resulted in more complex methods and greater reliance on various software to support research [33,35]. A 2015 survey showed ArcGIS as the most commonly used general software [35]. Other geostatistic tools include remote sensing and global positioning systems (GPS) [36].

A wide range of SDM methods has been developed to achieve various research objectives, with SDM algorithms based on machine learning now dominating the field [36,37,38,39]. Maximum entropy theory modeling (MaxEnt) [40] is currently the most commonly used software worldwide for modeling species distribution [34,38,41]. In Latin America, 73% of the analyses used MaxEnt [37]. Other commonly used methods are the genetic algorithm for rule-set prediction (GARP) [33], the generalized linear models (GLM) [42], Bioclimatic Prediction and Modeling System (BIOCLIM) [43], Random Forests [44], and generalized additive models (GAM) [42].

The majority of recent SDM research has focused on the impacts of climate change in order to assess the influence of the changing environment on species distributions and to make future predictions to guide mitigation programs [35,38]. The field of conservation biogeography now has practical methods for incorporating climate change within conservation studies [34]. Climate databases and global and regional climate models have been used in this regard. GCMs (General Circulation Models) simulate historical and projected climate variables for climate change studies [45]. Global climate models help determine climatic changes across the whole world, such as in the case of the Intergovernmental Panel on Climate Change (IPCC) studies, while regional climate models are usually used to study specific ecosystems [46,47,48].

Knowledge of the current state of wild cotton species in Mexico, including their geographic location and potential distribution, is scarce. Additional studies using precise spatial models are needed to provide essential potential distribution information for diverse applications in ecology and conservation [33]. Continental projections need to be complemented with smaller-scale studies to better understand geographic patterns and develop regional conservation strategies [23]. Some studies have contributed distribution maps of G. thurberi [9,26,27,29,49]. However, most of these do not include the possible effects of future climate on species distribution. Although cotton is adapted to growing in temperate, subtropical, and tropical environments in the world, climate change may affect its populations, both positively and negatively [3,32,50]. Species distribution models must thus be developed to predict future changes in the potential geographical range of cotton species and their effects on regional ecology and biogeography scenarios [25,51].

To help address these gaps, the current study aimed to investigate the distribution of G. thurberi in Mexico, to provide details on its potential habitats, and to estimate the effects of future climate change on its geographic distribution. To our knowledge, this is the first study of the potential distribution of G. thurberi in Mexico considering both current and future climates. The findings can be used in the ex situ and in situ conservation of this valuable species, especially in designing more precise collection and conservation strategies. The results will also be useful to other researchers performing related studies.

2. Materials and Methods

2.1. Area under Study

The study focuses mainly on northern Mexico (32°32′2.33″ N, 117°7′41.35″ W and 22°13′7.20″ N, 98°33′16.34″ W), which comprises the states of Sonora and Chihuahua, where populations of the species G. thurberi are currently concentrated, as well as the nearby states of Baja California, Baja California Sur, Coahuila, Nuevo León, and Durango (Figure 1). The rest of the country was included in the analysis of potential areas for future distribution.

2.2. Plant Data

The study used 76 collection sites of G. thurberi from the database of CONABIO-SNIB (National Commission for the Knowledge and Use of Biodiversity-National Biological Information System) [30]. In addition, from the same database, 85 sites of three cotton species that did not co-exist with G. thurberi in the presence data sites were used as pseudo-absences. Records of the presence and absence of the species are usually used in distribution studies [52]; however, when there are no absence data, pseudo-absence data or “false absences” may be used [53,54,55].

The occurrence data of G. thurberi were partitioned using a one-time data-splitting approach [56], whereby the data were split randomly, with one portion used for model construction and calibration and the other for validation [57]. While many studies have used the 70–30 ratio for splitting data for training and validation, there has been no consensus as to the best ratio to be applied, and data partitioning tends to be arbitrary [58,59]. In the present study, considering the small geographic range of G. thurberi and the compactness of the area where populations of the species G. thurberi are currently concentrated, we decided that the use of 15 occurrence sites (20%) of G. thurberi would be adequate for model building and calibration. The use of a small number of occurrence samples, while not considered ideal, has been useful in the modeling of rare or endangered species [59,60,61]. A study evaluating sampling size showed that the ability to model species effectively is strongly influenced by species ecological characteristics independent of sample size, and the modeling of species with compact spatial distributions tends to be easier than modeling species widespread in geographic and environmental space [60]. In this case, the geo-referenced G. thurberi collection sites are mostly located in the Mexican state of Sonora, and only a few sites are reported in the neighboring state of Chihuahua (Figure 1).

We also considered that for building a model with good generalization ability, validation will require sufficient data to reliably estimate model predictive performance [56,62]. Validation is particularly vital in model construction when spatial or temporal extrapolation is one of the objectives [62]. In this case, we needed a model with a high level of precision that may be extrapolated to other regions of the country presenting the suitable conditions required by the species and may be used to determine future distribution. Hence, we validated the model using datasets of 61 G. thurberi occurrence sites (80%) and 85 pseudo-absence sites.

The quality and accuracy of the plant data were analyzed as follows: ArcGIS 10.3 software of Esri, Inc.^® (1999–2014) was used to project the collection sites onto thematic layers with relevant location information [63,64]. The following information was obtained for each site: state, municipality, locality, latitude, and longitude. Then, the information was compared with that reported in the CONABIO-SNIB database [30] to check for any disparity which would require further analysis. No errors were observed in the data.

2.3. Environmental Parameters

Machine learning has made it possible to consider a wide range of parameters in SDM studies, which sometimes leads to duplication in parameterization and the use of irrelevant predictors. Several authors recommend the careful selection of essential and ecologically relevant variables for analysis [65,66]. Some studies have demonstrated good predictive ability even with only three variables [67]. For this study, we considered the following variables because of their importance to the cotton plant [2,15,68,69,70,71]: annual available soil water (mm year⁻¹, AASW), Flowering Growing Degree Days (FGDD), absolute minimum temperature (°C, Tminabs), and altitude (amsl, ALT).

To characterize the G. thurberi sites, the climate vector dataset (minimum and average temperature and monthly and annual averages of precipitation) with a scale of 1:1,000,000 of INEGI [72] was used. Likewise, the altitude thematic layer of INEGI [73] was used at a resolution of 1:10,000. From this information, it was possible to generate the following variables:

Annual Available Soil Water: this is the difference between annual precipitation and real annual evapotranspiration (AASW = Ppt − ETR) [74], where AASW is the annual available water in the soil (mm year⁻¹), Ppt is the mean annual precipitation (mm) and ETR is the real annual evapotranspiration (mm). ETR is calculated with the Turc method [75,76] as follows:

$$ETR=\frac{Ppt}{\sqrt{0.9+\frac{Pp{t}^2}{L{c}^2}}}$$

(1)

where ETR is the real annual evapotranspiration (mm year⁻¹), Ppt is the mean annual precipitation (mm), Lc is the calculated temperature factor (300 + 25 t + 0.05 t³) and t is mean annual temperature (°C).

Absolute Minimum Temperature: the minimum temperature variable (Tminabs), considered important in determining the distribution of species [68,69], is the theoretical minimum temperature (°C) in a 20-year period, calculated using the following equation of Prentice et al. [77]:

Tminabs = 0.006 Tc ² + 1.316 Tc − 21.9

(2)

where Tc in °C is the average minimum temperature of the coldest month of the year.

Flowering Growing Degree Days: the variable Growing Degree Days (GDD), also known as Heat Units (HUs) or thermal time, reflects the relationship between growth rate and temperature and refers to the specific amount of heat units that plants require to develop from one stage of development to another [15,77,78,79,80]. In the present study, we considered Flowering Growing Degree Days (FGDD), which refers to the necessary GDD to cover the flowering phase, the prime physiological stage of the cotton plant that determines the final productivity [2,71]. It is calculated as follows:

FGDD = ∑ max{0, (T m − Tb)}

(3)

where Tm is the mean daily temperature (°C) obtained from the monthly averages by sine curve interpolation between mid-months to obtain quasi-daily values [81] and Tb is the base temperature for flowering considering a value of 18 °C [82,83].

Altitude: altitude (ALT) is the vertical distance of a point on Earth from sea level. A gradual decrease in average temperature occurs from low to the highest altitude, with a change of less than 0.65 °C estimated for every 100 m of altitude in temperate and subtropical areas [84].

2.4. Model Construction and Calibration

The G. thurberi model was built using a Geographic Information System (GIS) platform. GIS tools make it possible to combine different types of data sources, carry out complex analyses, and generate maps [85,86].

First, we characterized six randomly chosen occurrence sites of G. thurberi to obtain information on factors influencing its distribution [87,88,89,90]. During this process, the values of AASW, FGDD, ALT, and Tminabs were obtained for each site through the Extract Multi Values to Points Tool of ArcGIS software. This provided the initial approximation of the profile of the species. The generated data matrix with the characteristics of the collection sites was used as a base in the succeeding stages of the study.

Figure 2 outlines the actions carried out during characterization and other phases of the study.

We derived the intervals formed by the minimum and maximum values of the variables AASW, FGDD, ALT, and Tminabs and used these to build the initial G. thurberi model. The Reclassify Tool in the Spatial Analyst Tools of ArcGIS was used to generate images with data falling within the range defined for each variable. From the overlay of these images, the model was built using the Map Algebra in the Spatial Analysis Tools.

The model was then calibrated, a process defined by Richardson and Berish [91] as the adjustment of parameters to improve the capacity of the model to predict the dataset used for verification. The objective was to reduce omission and commission errors. Using a confusion matrix, also known as contingency matrix or error matrix [52], we classified errors as follows: Omission errors are when the species was found, but the model predicted absence (false negatives); commission errors are when the species was not found, but the model predicted presence (false positives). To decrease the omission errors, the initial values of each variable for the presence data sensu lato (i.e., off-model sites that were within 1 km of the model) were verified in order to determine the variables and the ranges to be modified. Nine additional presence sites that had not been previously used were included in the calibration process.

2.5. Model Validation

During the validation stage, we compared the model with a different dataset in order to verify if the model had been developed as expected [91]. For validation methods requiring both presence and absence records, sampling pseudo-absences from the study area in place of real absence data has been applied in some studies [53,55]. We used 61 presence data of G.thurberi and, as pseudo-absences, 85 data of three other cotton species, namely, G. hirsutum (51 sites), G. turneri (12 sites), and G. armorianum (22 sites) from CONABIO-SNIB [30]. The data had not been used in model construction.

We considered that the conversion of data layers from vector to raster format might have changed the shape and dimension of some areas and that, during validation, the comparison of areas in the model with points (presence and absence data) might have caused some errors of omission. Hence, we measured the distance between each off-model presence site and the model. For the distance measurement, a mesh with a spacing of 1 km between centroids was generated using the ArcGIS Data Management Tools. Sites that were within 1 km were classified as presence sites sensu lato, while those located within the model, directly obtained with the Extraction Spatial Analysis Tool of ArcGIS, were classified as presence sites sensu strictum. Accuracy and concordance statistics were obtained with both the sensu lato and sensu strictum sites. We tested each model, using seven non-parametric statistical tests described in Section 2.8.

2.6. Sensitivity Analysis

Sensitivity analysis involves the study of the variation in model output, how it can be apportioned, qualitatively and quantitatively, to different sources of variation and how the model depends upon the given input [92]. A Jackknife leave-one-out approach was used to assess the importance of environmental variables and their ranges [41,61,93,94]. The less contributing variables and ranges were eliminated, and the key ones were determined based on their contribution to the modeling process [24].

During the first stage of the analysis, we determined the variables and their ranges that had an influence on the distribution of G. thurberi. The models that had reached the highest values in the accuracy and concordance tests during validation were selected, and their ranges of AASW, FGDD, ALT, and Tminabs were analyzed to determine the sites and/or localities with extreme values. Images with the data falling within the range defined for each variable were generated, and a total of six models were built, making it possible to determine the variation in the predictive power of each model through the inclusion or exclusion of the extreme values. All 61 sites (sensu strictum and sensu lato) of the six models were quantified, following the methodology described in the validation phase, and the precision and concordance statistics of each model were obtained.

During the second stage of the analysis, the variables with the greatest influence (dominant variables) in the distribution of G. thurberi were determined as follows: Among the six models developed and tested during the first stage of sensitivity analysis, we selected the one with the highest values in the statistical tests. Its variables and their ranges were placed in a double-entry matrix to obtain all the possible combinations of variables considered in model development. To verify if the order in which the images were handled in the Spatial Analyst Tools did not affect the distribution model, we compared models with the same variables, but with a different order, i.e., if a model was developed considering AASW and FGDD, two models were generated from this combination: AASW-FGDD and FGDD-AASW. The resulting 10 models underwent the same analysis as in the validation phase. Their level of precision and concordance was determined by means of the seven non-parametric statistical tests. The model that obtained the highest values in the tests was selected as the distribution model of G. thurberi in northern Mexico under current climate. To determine the potential distribution area of the species throughout Mexico, a projection of the model was performed considering the whole Mexican territory.

2.7. Potential Impact of Climate Change

The objective of this part of the study was to estimate the potential impact of possible future climate change on the current potential distribution of the species in northern Mexico and to determine suitable areas in the country where the species can be introduced or re-introduced. For this, we generated a potential distribution map for G. thurberi with current climate data and then a potential distribution map for the species under the future climate.

With the ArcGIS Spatial Analyst Tools, future climate layers were generated assuming temperature increases of 1.5 °C and 2.0 °C for the near-future (2021–2040) and mid-future (2041–2060) climates [95,96], and a precipitation reduction of 100 mm [97,98,99] under both climates. The IPCC CMIP6 projections under SSP5.85 scenario show very likely temperature increase ranges of 1.3–1.9 °C in the near-future and 1.9 to 3.0 °C in the mid-future. The GWL (Global Warming Level) of 1.5 °C is predicted to be reached or exceeded within 2011–2030 to 2029–2048, and 2.0 °C is expected within 2023–2042 to 2044–2063 [95,96,100]. Hence, we considered these increases in our study. Almazroui et al. [101], in their assessment of the CMIP6 models’ performance for different regions, including Mexico in North Central America, concluded that over the 21st century, temperature increases will be large in the US and smaller but significant over Central America.

Mexico is predicted to be among the regions with exacerbated water deficits due to increased temperature and reduced rainfall [96,98,100]. Regional climate change scenarios from other studies show contrasts in the projected climate changes between the regions of Mexico, with precipitation changes ranging from −20.3% to 13.5%, depending on the scenario and analysis period [98].] Ureta et al. [102] mentioned that the GCMs vary greatly in their precipitation projections for Mexico (e.g., −67.6 to −18.2 mm for mid-future climate). Almazroui et al. [101] have likewise seen a wide spread of uncertainly among the CMIP6 models, particularly for the high-emission scenarios. In this study, we generated the future precipitation layers considering a 100 mm (approximately 20%) reduction in precipitation for northern Mexico [97,98,99], where existing populations of G. thurberi occur. Observed climatic changes in the Sonoran Desert region include warming temperatures and precipitation deficits [96], increased drought frequency, and a higher number of heat-stress days, signaling an arid, water-constrained climate [103] that is expected to remain, if not intensify, in the near and mid future. Severe drought has already been reported in Sonora and Chihuahua and the rest of the northern region [104,105].

The mean annual temperature and annual precipitation of the two future climates were generated using a baseline climate 1970–2000 with a 30” arc resolution [72]. The resulting climate layers were used as input in the generation of the thematic layers of future ETP (Equation (1)), AASW (Ppt-ETP), and FGDD (Equation (3)) variables. In the calculation of FGDD, in addition to tb, a threshold temperature or terminal stress of 30 °C was also considered [83,106,107]. The Map Algebra in the Spatial Analyst tools was used in this procedure.

The G. thurberi distribution model under current climate was projected onto the new future climate layers to obtain the potential distribution of the species under future climate. Three areas were identified: (1) high impact areas, where the species appears at present but will no longer be suitable in the future; (2) low impact areas, where the species can potentially exist in both present and future climates; (3) new suitable areas, where the species could potentially exist in the future, but are currently not suitable for natural occurrence [108].

The high and low impact areas were determined by projecting the current distribution model onto the digital layers of variables under the two future climate scenarios. To determine the areas that the species would gain in the future (new suitable areas), the current collection sites were projected onto the low-impact or stable areas, thereby obtaining, through Spatial Analyst Tools (Extract Multi Values to Points), the ranges of future AASW and FGDD. The ALT values of the populations projected to survive were also obtained so as to define the characteristics of the future suitable and unsuitable areas. With this information, distribution models of new suitable areas of the species were generated under each of the future climate scenarios. These models were then projected throughout the Mexican territory to generate maps.

2.8. Statistical Analysis

Seven non-parametric statistical tests were applied in the validation and sensitivity analysis processes: kappa index, overall accuracy, sensitivity, specificity, positive predictive power, negative predictive power and odds ratio [109,110,111,112]. The kappa index is a measure of agreement indicating the reliability of the data obtained, which ranges from 0, signifying no agreement between observed and predicted distributions, to 1, indicating perfect agreement. Overall accuracy is the measurement of correct classification made by the model. Sensitivity is the conditional probability that case X is correctly classified as true. Specificity is the conditional probability that a case that is not X is correctly classified as false. Positive predictive power assesses the probability that a case is X if the algorithm classifies the case as X. Negative predictive power assesses the probability that a case is not X if the algorithm does not classify the case as X. Odds ratio is the ratio of correctly assigned cases to incorrectly assigned cases. The following equations were used:

$$\mathrm{Kappa}\mathrm{Index}:\left(\left(\mathrm{a}+\mathrm{d}/\mathrm{n}\right)-\left(\left(\mathrm{a}+\mathrm{b}\right)\left(\mathrm{a}+\mathrm{c}\right)+\left(\mathrm{c}+\mathrm{d}\right)\left(\mathrm{d}+\mathrm{b}\right)\right)/{\mathrm{n}}^2\right)/\left(1-\left(\left(\mathrm{a}+\mathrm{b}\right)\left(\mathrm{a}+\mathrm{c}\right)+\left(\mathrm{c}+\mathrm{d}\right)\left(\mathrm{d}+\mathrm{b}\right)\right)/{\mathrm{n}}^2\right)$$

(4)

$$\mathrm{Overall}\mathrm{accuracy}=\mathrm{a}+\mathrm{d}/\mathrm{n}$$

(5)

$$\mathrm{Sensibility}:=\mathrm{a}/\mathrm{a}+\mathrm{c}$$

(6)

$$\mathrm{Specificity}=\mathrm{d}/\mathrm{b}+\mathrm{d}$$

(7)

$$\mathrm{Positive}\mathrm{Predictive}\mathrm{power}=\mathrm{a}/\mathrm{a}+\mathrm{b}$$

(8)

$$\mathrm{Negative}\mathrm{Predictive}\mathrm{power}=\mathrm{d}/\mathrm{d}+\mathrm{c}$$

(9)

$$\mathrm{Odds}\mathrm{Ratio}=\mathrm{ad}/\mathrm{cb}$$

(10)

where: “a” and “d” are presence and absence, respectively, correctly predicted by the model, and “b” and “c” indicate presence and absence, respectively, wrongly predicted by the model; “n” is the sum of “a, b, c, d”.

3. Results and Discussion

3.1. Model Construction and Calibration

The characterization of the occurrence sites of G. thurberi showed that despite being in a relatively compact area, territorially, in the states of Sonora and Chihuahua, G. thurberi presented wide ranges in the four variables (

Table 1. First approximation of characteristics of occurrence sites of wild cotton species Gossypium thurberi in the states of Sonora and Chihuahua, Mexico.

Variable	Range
Annual Available Soil Water (AASW)	0–90 mm year−1
Flowering Growing Degree Days (FGDD)	332–680
Altitude (ALT)	772–1486 amsl
Minimum Absolute Temperature (Tminabs)	−23.2–−15.2 °C

). For example, the areas it occupied ranged from those with zero water availability (i.e., annual precipitation equal to or less than actual evapotranspiration) to areas with 90 mm year⁻¹ of AASW, thus demonstrating the ability of the species to adapt to the arid and semi-arid zones of northern Mexico [11].

Of the three models initially built (Figure 3 and

Table 2. Gossypium thurberi models developed during the initial modeling construction stage and variables considered in their construction.

Model	Variable	Range
Mb_1	AAW †	0–90 mm year−1
	FGDD ‡	332–680
Mb_2	AASW	0–90 mm year−1
	FGDD	332–680
	ALT §	772–1486 amsl
Mb_3	AASW	0–90 mm year−1
	FGDD	332–680
	ALT	772–1486 amsl
	Tminabs ¥	−23.2–−15.2 °C

^† AASW: Annual available soil water. ^‡ FGDD, Flowering growing degree days. ^§ ALT: Altitude. ^¥ Tminabs: Absolute minimum temperature.

), the Mb_2 Model had the highest values in the kappa index, precision, sensitivity, and negative predictive power tests, with values of 0.72, 0.87, 0.96, and 0.98, respectively (Table S1). These values were similar to those obtained by the Mb_3 Model, which had the highest value (186) in the odds proportion test. Thus, the calibration process considered the variable ranges of both Models Mb_2 and Mb_3 and of the nine additional sites used for calibration. The extreme values of some sites caused some ranges to widen towards the upper end. Therefore, two types of models were built (

Table 3. Models generated during the calibration process to reduce errors of omission and commission in the distribution model of Gossypium thurberi in Mexico.

Model	Variable	Range
Cal_1	AAW †	0–90 mm year−1
	FGDD ‡	332–732
	ALT §	556–1474 amsl
Cal_2	AASW	0–90 mm year−1
	FGDD	332–732
	ALT	556–1474 amsl
	Tminabs ¥	−23.2–−15.2 °C
Cal_3	AASW	0–92 and 167–234 mm year−1
	FGDD	332–680 and 700–732
Cal_4	AASW	0- 92 and 167–234 mm year−1
	FGDD	332–680 and 700–732
	ALT	556–1486 and 1800–1810 amsl
Cal_5	AASW	0–92 and 167–234 mm year−1
	FGDD	332–680 and 700–732
	ALT	556–1486 and 1800–1810 amsl
	Tminabs	−23.2–−15.2 and −27.0–−24.5 °C

^† AASW: Annual available soil water. ^‡ FGDD, Flowering growing degree days. ^§ ALT: Altitude. ^¥ Tminabs: Absolute minimum temperature.

): (1) models with a single range containing the values of the sensu lato sites only and (2) models with two ranges of values from the sensu lato sites and the nine additional presence sites.

A visual analysis was made of the distributions predicted by the models with a single range of variables, i.e., Models Cal_1 and Cal_2 (Figure 4a,b), whose only difference between them was the inclusion of the Tminabs variable in Model Cal_2 (Table 3). The Tminabs represents the minimum temperature expected in the next 20 years [77]. In these two models, the projected area disappeared in the lower part of Baja California Sur as well as in the northern part of Chihuahua. The largest area adjustments (Figure 4c–e) occurred when the models included two ranges in each variable (Table 3).

3.2. Model Validation

All the models (

Table 4. Results of the concordance and accuracy tests of Gossypium thurberi models for northern Mexico (after calibration).

Model	Concordance and Accuracy Tests
	Kappa Index	Overall Accuracy	Sensitivity	Specificity	Positive Predictive Power	Negative Predictive Power	Odds Rate
Cal_1	0.67	0.84	0.84	0.84	0.77	0.89	28
Cal_2	0.70	0.86	0.92	0.83	0.72	0.95	52
Cal_3	0.65	0.83	0.77	0.88	0.84	0.82	24
Cal_4	0.64	0.83	0.86	0.81	0.70	0.92	27
Cal_5	0.71	0.86	0.96	0.82	0.70	0.98	99

) obtained kappa values in the substantial accuracy category (0.61–0.80), according to the Landis and Kock scale [113].

Model Cal_5 showed the smallest area, projecting the distribution of G. thurberi only in the states of Sonora, Chihuahua, Coahuila, and Nuevo León (Figure 4e). It obtained the highest values of kappa index (0.71) and sensitivity (0.96) (Table 4). The inclusion of the ALT and Tminabs variables in the model increased its reliability by 6% (0.71 vs. 0.65). The highest value in specificity (0.88) was obtained by Model Cal_3, whose kappa index was 0.65. As both Cal_5 and Cal_3 models best represented the current potential distribution of G. thurberi in northern Mexico, their ranges of variables were used in the sensitivity analysis.

3.3. Sensitivity Analysis

3.3.1. Significant Environmental Variables and Ranges in the Distribution of Gossypium thurberi

The variables and ranges of the six models built during the first stage of the sensitivity analysis are presented in

Table 5. [M1] Models generated during the sensitivity analysis process to determine significant variables and their ranges in the distribution of Gossypium thurberi in Mexico.

Model	AAW † ( mm year −1)	FGDD ‡	ALT§ (amsl)	Tminabs ¥ (°C)
I_SA II_SA	0–92, 167–234 0–92 0–92, 167–234 0–92 0–92, 167–234 0–92, 167–234	332–680 332–680, 700–732 332–680 332–680, 700–732 332–680, 700–732 332–680, 700–732	556–1486, 1800–1810 556–1486, 1800–1810 556–1486 556–1486, 1800–1810	−23.2–−15.2, −27.0–−24.5 −23.2–−15.2, −27.0–−24.5 −23.2–−15.2, −27.0–−24.5 −23.2–−15.2
III_SA IV_SA V_SA VI_SA

† AASW: Annual available soil water. ‡ FGDD, Flowering growing degree days. § ALT: Altitude. ¥ Tminabs: Absolute minimum temperature.

, and their predicted distribution is shown in Figure 5. The accuracy and concordance test results (

Table 6. Results of the concordance and accuracy tests applied to models during the first stage of the sensitivity analysis.

Model	Concordance and Accuracy Tests
	Kappa Index	Overall Accuracy	Sensitivity	Specificity	Positive Predictive Power	Negative Predictive Power	Odds Ratio
I_SA II_SA III_SA IV_SA V_SA VI_SA	0.68 0.65 0.75 0.71 0.71 0.66	0.84 0.83 0.88 0.86 0.86 0.84	0.82 0.77 0.96 0.96 0.96 0.95	0.86 0.88 0.85 0.82 0.82 0.80	0.80 0.84 0.75 0.70 0.70 0.66	0.87 0.82 0.98 0.98 0.98 0.98	27 24 127 99 99 99

) showed that Model III_SA (Figure 5c) had the highest values in five tests. The inclusion of ALT and Tminabs in this model decreased the commission error. Its kappa index of 0.75 puts it in the substantial accuracy category (0.60–0.80), according to the Landis and Kock [113] scale, making it the most reliable model in this stage of the study. As Model III_SA contains the significant environmental variables and ranges in the distribution of G. thurberi in northern Mexico, it was used in the second step of the sensitivity analysis to determine the dominant variable in the distribution of the species.

3.3.2. Dominant Variables in the Current Distribution of Gossypium thurberi

According to Gertseva et al. [114], the organisms in a population can change their distribution depending on the values of the dominant factor but not in response to the values of sub-dominant factors. To determine the dominant variable in this study, a double-entry matrix was developed with all four variables of Model III_SA in order to generate all possible combinations of variables (

Table 7. Results of the concordance and accuracy tests applied to ten models of Gossypium thurberi in northern Mexico generated during the second stage of the sensitivity analysis.

Model	Concordance and Accuracy Tests
	Kappa Index	Overall Accuracy	Sensitivity	Specificity	Positive Predictive Power	Negative Predictive Power	Odds Rate
FGDD ‡ -AASW †	0.71	0.86	0.82	0.88	0.84	0.87	3
FGDD-ALT §	0.68	0.84	0.82	0.86	0.80	0.87	27
FGDD-TAB ¥	0.72	0.86	0.85	0.87	0.82	0.89	38
AASW-ALT	0.74	0.88	0.89	0.87	0.80	0.93	54
AASW-TAB	0.79	0.90	0.91	0.89	0.84	0.94	82
ALT-TAB	0.69	0.85	0.84	0.85	0.79	0.89	31
AASW-ALT-TAB	0.70	0.86	0.88	0.84	0.75	0.93	40
FGDD-ALT-TAB	0.69	0.85	0.85	0.85	0.77	0.91	32
AASW-FGDD-TAB	0.78	0.90	0.94	0.87	0.80	0.96	112
AASW-FGDD-ALT	0.81	0.91	0.96	0.88	0.82	0.98	189

^‡ FGDD, Flowering growing degree days (332–680). ^† AASW: Annual available soil water (0–92, 167–234 mm year⁻¹). ^§ ALT: Altitude (556–1486,1800–1810 amsl). ^¥ Tminabs: Absolute minimum temperature (−23.2–−15.2, −27.0–−24.5 °C).

).

The results of the concordance and accuracy tests applied to the 10 models (Table 7) showed that the model combining the three variables AASW, FGDD, and ALT obtained the highest values in the kappa index (0.81), general precision (0.91), sensitivity (0.96), negative predictive power (0.98), and odds proportion (189) tests (Table 7). Its kappa index of 0.81 puts it in the almost-perfect accuracy category (0.81–1.0), according to the Landis and Kock scale [113]. We therefore conclude that the dominant variables and their ranges in the current potential distribution of G. thurberi are AASW (0–92, 167–234 mm year⁻¹), FGDD (332–680), and ALT (556–1486, 1800–1810 amsl). These variables and their ranges of values were used to generate the potential distribution map of the species in northern Mexico under current climate. In addition, they were used in a projection to determine the potential distribution areas of the species in the whole Mexican territory.

3.4. Current Potential Distribution

Maps of the current potential distribution of G. thurberi in northern Mexico and the whole country (Figure 6a,b) show that the current potential distribution of the species is wide in the country. Hence, we inferred that G. thurberi may be present and/or introduced in other areas with similar conditions as those of locations of current populations. To verify this, we considered the ex situ conservation garden of wild cotton species maintained by INIFAP (National Research Institute for Forestry, Agriculture, and Livestock Production) in Campo Experimental Iguala (red dot in Figure 6b) in Guerrero in southern Mexico [10]. The projected distribution areas of the species beyond the northern region covered the area of the said garden, where, to date, the species G. thurberi is growing and developing satisfactorily. The southern state of Guerrero has a warm subhumid climate Aw0, according to the Koppen classification modified by García [115]. The precipitation of the driest month is 0−60 mm. The altitude at which this ex situ garden is located is 721 amsl. The average annual temperature and the temperature of the coldest month are greater than 22 °C and 18 °C, respectively. This illustrates that G. thurberi can also develop and grow satisfactorily in regions with minimum temperatures above 0 °C in the coldest month of the year.

We compared our distribution map for the states of Sonora and Chihuahua (polygon, scale 1:1,000,000, GCS_WGS_1984) with those of previous studies on G. thurberi [9,49,116] by superimposing the shape files of the studies in ArcGIS. It was observed that in the model generated by CONABIO [116] (polygon, scale 1:1,000,000, GCS_WGS_1984) (Figure 7a) using 23 climatic variables from WorldClim and the GARP algorithm, the defined potential distribution covered the central part of the state of Sonora which fully coincides with what was modeled in our study. However, the CONABIO model [116] did not present any distribution in northern Sonora, such as what was shown in our study, possibly because only collections recorded in the period 2007–2008 were considered, and these did not include northern Sonora.

A later study was done by Cuervo-Robayo et al. 2019 [49], using the MaxEnt algorithm with seven bioclimatic variables, as well as land cover (bare soil and cultivated areas) and topographical information (altitude and slope). The modeled potential distribution of the species in the states of Sonora and Chihuahua (polygon, scale 1:1,000,000, GCS _WGS_1984) matched that of our study (Figure 7b). However, they projected the potential distribution of G. thurberi as a compact polygon whose total area is included in the distribution of the species (Figure 7b). A similar form of compact distribution polygon for G. thurberi was reported by Sanchez -Reyes et al. [9], who also used the MaxEnt algorithm and seven bioclimatic variables from WorldClim, plus elevation, soil type, and aridity index. In contrast, our model showed specific areas or zones of distribution, which might facilitate directed field exploration. Also, Cuervo-Robayo et al. 2019 [49] reported that in some cases they eliminated overestimated areas, whereas we considered those as potential distribution areas of the species. According to Raxworthy et al. [117], the regions that the model overpredicts show potential for identifying unknown distributional areas. The G. thurberi maps of Sanchez-Reyes et al. [9] and Cuervo-Robayo et al. [49] were part of larger studies involving multiple species. It would be interesting to compare the performance of these models, especially as Cuervo-Robayo et al. [49] also used CONABIO-SNIB data [30], although with a 70–30 data splitting ratio and MaxEnt. Unfortunately, no detailed statistical test results could be accessed. A comparative analysis could contribute to the evaluation of different approaches to data splitting and species distribution modeling [59].

3.5. Potential Distribution of Gossypium thurberi under Future Climates

3.5.1. Future Potential Distribution

Our study results estimate the current potential distribution area of the species in northern Mexico (Baja California, Baja California Sur, Sonora, Chihuahua, Nuevo Leon, Coahuila, and Durango) to be 112,727 square kilometers (Figure 8a,

Table 8. Potential suitable areas for Gossypium thurberi in northern Mexico under current and future climates.

Model	Area (Square km)	Percentage (%)
Current	112,727	100.0
Near-Future (2021–2040) ‡	25,307	−77.0 §
Mid-Future (2041–2060) ¥	16,017	−86.0 §

^‡ Near-future climate with 1.5 °C temperature increase and 100 mm precipitation reduction. ^¥ Mid-future climate with 2.0 °C temperature increase and 100 mm precipitation reduction. ^§ Percentage of reduction in relation to the current potential area.

). It is projected to be reduced to only 25,307 square kilometers under a near-future (2021–2040) climate presenting a 1.5 °C increase in mean annual temperature and a 100 mm reduction in precipitation (Figure 8b). It will be further reduced to 16,017 square kilometers under a mid-future (2041–2060) climate with a 2.0 °C increase in mean annual temperature and a 100 mm reduction in precipitation (Figure 8c). In short, a loss of 77 to 86% of the current potential distribution of the species is estimated under the two future climates.

This is alarming since, as shown in Figure 8b,c, the areas with the greatest loss (high impact areas) are in Sonora, where several populations of the species are currently located [30]. These areas border the state of Arizona, where a large area with a population of this species also exists [9,26] that might be similarly affected.

Based on our analysis, the high-impact areas that will not have the capacity to house populations in the future will have an average of −6.1- and −6.2-mm year ⁻¹ of AASW and 430 and 357 annual FGDD under the near-future and mid-future climates, respectively. The low impact distribution areas will present an average of 20.4- and 13.8-mm year ⁻¹ of AASW and 441 and 389 annual FGDD under near-future and mid-future climates, respectively (Figure 8b,c,

Table 9. Characteristics of suitable and unsuitable areas for the establishment of wild cotton species Gossypium thurberi in Mexico under near-future (2021–2040) and mid-future (2041–2060) climates with increases of 1.5 °C and 2.0 °C in mean annual temperature, respectively, and a 100 mm reduction in precipitation under both climates.

	Climate Change Scenario
	Near Future (2021–2040) ‡			Mid Future (2041–2060) ¥
	Suitable areas			Suitable areas
	AAW †	FGDD *	ALT §	AASW	FGDD	ALT
Mean	20.4	441	1034	13.8	389	1059
Maximum	63.8	638	1561	55.6	547	1561
Minimum	0.2	332	500	0.2	242	770
	Unsuitable areas			Unsuitable areas
	AASW	FGDD	ALT	AASW	FGDD	ALT
Mean	−6.1	430	927	−6.2	357	888
Maximum	−0.9	699	1438	−1.2	608	1438
Minimum	−8.5	271	562	−8.9	180	562

^‡ Near-future climate with 1.5 °C temperature increase and 100 mm precipitation reduction. ^¥ Mid-future climate with 2.0 °C temperature increase and 100 mm precipitation reduction. ^† AASW: Annual available soil water (mm year⁻¹). * FGDD: Flowering growing degree days. ^§ ALT: Altitude (amsl).

). Considering these results, it can be inferred that AASW will influence the species distribution most in the future. This is in agreement with IPPC CMIP6 regional projections indicating widespread soil moisture decline during summer across North America, with Mexico having the largest declines [100]. In recent decades, the demands of increased human population and agriculture have exceeded water availability in some areas in Sonora and other arid and semi-arid provinces in northern Mexico that are highly dependent on water availability for agriculture and ecosystem productivity [118].

3.5.2. Losing and Gaining Suitable Climatic Areas under Future Climatic Scenarios

Berry et al. [119] mentioned that modeled species could be categorized as those losing suitable climate space, those gaining it, and those showing little or no change. In the case of G. thurberi, it can lose suitable climate space but also gain it (Figure 9). For example, in northern Mexico, specifically in the states of Sonora and Chihuahua, the species will lose a large part of its current potential area (high impact areas) under future climates; however, a part remains on the border of both states (low impact areas). Moreover, this species may in the future have new suitable climatic spaces in the southern part of the state of Sonora that borders Chihuahua (Figure 9a,b).

In central and southern Mexico, the same behavior of loss and gain of suitable climate space is seen (Figure 9), especially in the states of Jalisco and Oaxaca, where the potential surface is currently very low but, under future climates, more extensive areas will be suitable for the species. This is projected to occur mainly in the central part of Jalisco and on the Jalisco–Nayarit border (Figure 9a–c).

On the other hand, it is observed that, as in the north, the impact of climate change will adversely affect large areas of central and southern Mexico, with greater impact on the high-temperature areas of Michoacan, where the suitable climate space will almost completely disappear. A similar substantial loss in potential distribution has been projected for avocado in these areas due to increased temperature and decreased precipitation [120]. It is possible that the high temperatures in these areas in the future will have a direct effect on the real annual evapotranspiration rate (ETR), which is a function of long-term mean precipitation and temperature, calculated using the Turc method in this study [75,76].

3.5.3. Refuge Areas for the Species under Future Climatic Scenarios

To avoid species extinction, it is important to locate refuge areas where the species can be conserved and/or introduced, considering habitat suitability. Figure 9b,c show at least three areas in Mexico where the species can be protected in the future. These refuge zones have wide areas that will continue to present the minimum climatic conditions necessary for the survival of the species even under more severe climate.

One possible refuge area is in the northern region, on the Sonora–Chihuahua border, where some populations currently exist (Figure 9b,c). Another is in the north-central state of San Luis Potosí. The third area is on the boundaries of the states of Guerrero, Puebla, and Morelos. This refuge zone is distinguished by being a large compact area in a section of the Trans-Mexican Volcanic Belt where the Sierra Madre Oriental and Occidental Sierras converge. Its location indicates that the species will have the opportunity for survival in high areas, although within a certain limit, since, as shown in Table 9, the future altitude range of the species is estimated to be 500–1561 amsl and 770–1561 amsl under near-future and mid-future climates, respectively. As our findings indicate, the joint interaction of the parameters of altitude, temperature, and precipitation plays an important role in the current and future distribution of G. thurberi.

3.6. Considerations

Predictions of future habitat suitability using SDMS and climate scenarios have their uncertainties and limitations and thus must be viewed as simply estimates of a species’ potential distribution [60,108]. It must be considered that our study results are based on a model that does not incorporate all the factors influencing distribution, some of which are more local and dynamic in nature. The maps resulting from our study show areas of potentially suitable climate space that were defined considering selected variables and climatic scenarios. Future distributions of G. thurberi could be affected by a different set of climate conditions or other factors influencing species survival that are not within the scope of this study.

Despite the limitations of the study, the findings can contribute to formulating preliminary guidelines on actions required for safeguarding the future of G. thurberi, particularly in matters of collection efforts and species establishment in suitable habitats. They can help conservation managers plan adequate preservation and expansion programs, considering the direction of distributional change under future climate scenarios. As mentioned by other authors [121], the greatest threat to biodiversity is habitat loss, and the best way of safeguarding biodiversity is habitat protection.

4. Conclusions

The current potential distribution of G. thurberi in northern Mexico, estimated at 112,727 square kilometers, will be greatly reduced by 77 and 86% under two future climate scenarios that consider increases in mean annual temperature of 1.5 °C and 2 °C in the near- and mid-future, respectively, and a 100 mm reduction in average annual precipitation under both scenarios. The greatest reduction is projected to occur in areas of Sonora (Mexico) adjoining Arizona (USA), where the largest populations of the species are currently concentrated. However, the species will also gain new suitable areas in the northern, central, and southern regions of Mexico that may serve as refuge areas for the species in the face of climate change.

Annual available water in the soil (mm year⁻¹), flowering growing degree days, and altitude (amsl) jointly influence the current potential distribution of G. thurberi in Mexico, with annual available soil water being possibly the dominant factor under future climate change. Based on our analysis, the areas that may continue to harbor populations of G. thurberi under future climate will present AASW values of 0.2–55.6 mm year⁻¹, FGDD of 242–547, and ALT of 550–1561 amsl.

Our study contributes an approach to the modeling of the potential distribution of G. thurberi under current and future climates. It shows the influence of water, temperature, and altitude on its distribution, adding support to the tenet that water–energy dynamics [122,123] may explain the spatial distribution of wild cotton species.

Considering the synergistic relationship between climate change and other environmental factors [118], a projected potential spatial distribution may seem an oversimplification of complex processes [120]. However, it is useful to evaluate the areas where a species will no longer experience optimal climatic conditions and to estimate suitable areas in the future [96,120]. The findings of this study can aid species-specific conservation programs, for planning the ex situ and in situ conservation of G. thurberi. They can also be of interest to researchers performing related studies on this valuable species.