Does size matter? An analysis of the niche width and vulnerability to climate change of fourteen species of the genus Crotalus from North America

The niche comprises the set of abiotic and biotic environmental conditions in which a species can live. Consequently, those species that present broader niches are expected to be more tolerant to changes in climatic variations than those species that present reduced niches. In this study, we estimate the amplitude of the climatic niche of fourteen species of rattlesnakes of the genus Crotalus to evaluate whether those species that present broader niches are less susceptible to the loss of climatically suitable zones due to the projected climate change for the time period 2021–2040. Our results suggest that for the species under study, the breadth of the niche is not a factor that determines their vulnerability to climatic variations. However, 71.4% of the species will experience increasingly inadequate habitat conditions, mainly due to the increase in temperature and the contribution that this variable has in the creation of climatically suitable zones for most of these species.


INTRODUCTION
Global climate change is one of the main factors that impact biodiversity and the distribution of species (Barnosky et al., 2011). Each species has a tolerance to various environmental factors, and when this tolerance is exceeded, the species cannot optimally carry out their life cycle (Peters, 1990;Walther et al., 2002;Hardy, 2003;Dawson & Spannagle, 2009). When this occurs, the distribution and abundance of the species is altered (Hughes, 2000;Peterson et al., 2005;Root et al., 2005;Parmesan, 2006), and in some cases, it can result in the direct disappearance of some species and populations (Walther et al., 2002;Thomas et al., 2004). This in turn creates conditions that could modify the structure in the composition of species in the ecosystem and, consequently, disturb the ecological balance of a landscape (Gray, 2005;Walther, Beißner & Burga, 2005).
Niche modeling provides a predictive measure about how the climatic suitability of a species may change under different climate change scenarios (Morin & Lechowicz, 2008;Thuiller, Lavorel & Araújo, 2005b;Lawler et al., 2009). Currently, most niche models have been developed from a correlative approach, particularly when more than one species is involved (Hijmans & Graham, 2006). In this approach, the environmental variables that characterize the places where a species occurs (or is absent) are used to develop correlative models that analyze the effect of climate change on the species' climatic suitability (Wiens et al., 2009).
Rattlesnakes of the genus Crotalus are widely distributed across the New World from southern Canada to Argentina (Campbell & Lamar, 2004). There are currently 53 recognized species, with the greatest number found in Mexico (Sánchez et al., 2020). Various authors point out that temperature and precipitation are important factors in the ecology of the species of this genus (Paredes-García, Ramírez-Bautista & Martínez-Morales, 2011;Sunny et al., 2019;Yañez-Arenas et al., 2020). As such, Crotalus represent a good model to predict the response of snake species to climate change. However, there are few studies that evaluate the effects that these environmental variations will have on the future distributions of species of this genus (Greene & Campbell, 1993;Gibbons et al., 2000). In this regard, and under the criterion that the niche comprises a set of environmental conditions in which a species may exist (Gaston, Blackburn & Lawton, 1997), it has been suggested that those species with broader niches could be less vulnerable to abrupt environmental variation under anthropogenic climate change. By contrast, those species with narrow niches would be particularly threatened by climatic disturbances (Brown, 1984;Johnson, 1998;Boyles & Storm, 2007;Botts, Erasmus & Alexander, 2013;Ozinga et al., 2013).
From this perspective, the question arises: can the breadth of niche, by itself, be considered as a determining factor that helps to predict the vulnerability of Crotalus species to climate change? Few studies have provided sufficient evidence to answer this question and thus the effects that climate change will have on each of the species of this genus, remain unknown (Greene & Campbell, 1993;Gibbons et al., 2000). The present study aims to analyze whether there is a relationship between niche width and vulnerability to climate change, projected for the period 2021-2040, in a sample of fourteen species of the genus Crotalus distributed in North America. This information is of great relevance for the establishment and development of conservation strategies for species of the genus Crotalus.

Presence data
We obtained geographical data of occurrences of 14 species of Crotalus, including C. atrox, C. basiliscus, C. cerastes, C. enyo, C. intermedius, C. lepidus, C. molossus, C. pricei, C. ravus, C. ruber, C. scutulatus, C. tigris, C. viridis, and C. willardi (following the taxonomy of Campbell & Lamar, 2004). We obtained geographical data from published scientific information (scientific articles, scientific notes, books), scientific collections from Mexico and other countries (Table S1), information generated by the National Commission of Protected Natural Areas (CONANP), as well as from the database of the Global Biodiversity Information Facility (GBIF; http://www.data.gbif.org). We selected these 14 species of the genus Crotalus because, after the geographic data purification process, they were the species that had the most complete base of geographic records with the best distributed geographic records in the known range of these species, reflecting with greater precision the total range of the species under study (Campbell & Lamar, 2004). As has been previously demonstrated, the precision of geographic records is of great relevance in the performance of species distribution models (Hefley et al., 2014;Fei & Yu, 2015;Velásquez-Tibatá, Graham & Munch, 2015). Data 'cleanliness' is particularly important for data coming from species distribution data warehouses such as GBIF (Hijmans & Elith, 2013). Using the ''dismo '' library (Hijmans et al., 2017) in the statistical software R (version 3.1.3, R Core Team, 2015), we checked the geographic projections of each record and eliminated duplicate records. We further cross-checked coordinates through visual inspection (Hijmans et al., 1999) and assessed sampling bias by subsampling the geographic records (Hijmans & Spooner, 2001;Phillips et al., 2009). Records with unreliable coordinates (according to the known distribution of the species; Campbell & Lamar, 2004) were removed from the database. In total, we generated a data set with 4,813 presence points (C. atrox = 1,241, C. basiliscus = 125, C. cerastes = 676, C. enyo = 135, C. intermedius = 41, C. lepidus = 239, C. molossus = 516, C. pricei = 76, C. ravus = 52, C. ruber = 568, C. scutulatus = 610, C. tigris = 72, C. viridis = 429, and C. willardi = 33; Fig. 1).

Climatic variables
Current weather data for North America was obtained with a resolution of 2.5 min (∼5 km) from the WorldClim database (version 2). This is an online database with 19 bioclimatic variables derived from monthly averages  of temperature and precipitation (Fick & Hijmans, 2017). We selected a subset of these variables on the basis of ecological theory and subsequently reduced, when necessary, through statistical analysis (Austin, 2007). In the preselection of the variables related to temperature, we considered those proposed by Rodder & Lotters (2009), who suggested that this set of variables were of great ecological relevance, particularly for those taxa limited by thermoregulation, such as squamates. The variables related to precipitation included descriptors that have been mentioned as key factors for the species of the genus Crotalus, which may become more relevant when thermal conditions are not optimal, for example in periods of time with extreme temperatures (Glaudas, 2009;Phadnis et al., 2019). Subsequently, to eliminate variables that provide similar information, we developed a Pearson correlation matrix (r < 0.7) to reduce the collinearity error. After this process, the retained variables were Annual Mean Temperature (bio1), Mean Diurnal Range (bio2), Mean Temperature of Wettest Quarter (bio8), Annual Precipitation (bio12), Precipitation of Wettest Month (bio13), Precipitation of Driest Month (bio14), Precipitation Seasonality (bio15), Precipitation of Warmest Quarter (bio18) and Precipitation of Coldest Quarter (bio19). In general, the bivariate correlation analysis was carried out by providing information on the 19 climatic variables to the presence records of the species under study. In our case, the climatic information was provided to 10,000 randomly distributed geographic points in the distribution area of the species under study to avoid discarding areas with relevant climate information (non-repetitive) (Becerra-López et al., 2016).

Climate profile and niche range
With the selected variables of the current climate, a Principal Component Analysis (PCA) was carried out in the software R using the ''ecospat'' library (Broennimann et al., 2014) to identify the climatic profile within the distribution area of the species under study.
We also evaluated the climate profile for the climate change models BCC-CSM2-MR, CNRM-CM6-1, and IPSL-CM6A-LR, considering the shared socio-economic pathway 5 8.5 W/m 2 (SSP5 8.5) proposed for the period 2021-2040. These climate models were randomly selected from a total of eight models.
For each selected variable, we then performed an Analysis of Variance and Tukey's post hoc tests to evaluate if there were statistical differences between the current climate data and the climate change scenarios. Subsequently, in the software R, the distribution of the species under study in the climatic space (niche range) were identified through a PCA using the nine current climate variables used in this study, following the methodology proposed by Becerra-López et al. (2020). This representation of the records of the species in a climatic context is based on the Hutchinson duality that indicates that there are two spaces, the geographic one and the multidimensional abstract space, denoted by climatic variables that establish the conditions in which a species can simply exist (Colwell & Rangel, 2009).
For the selection of SSP5 8.5 W/m 2 , we took into account that the narrative of this route considers a socioeconomic development driven by fossil fuels, which implies a scenario with increasing CO 2 emissions (Riahi et al., 2016;Kriegler et al., 2017). Considering that fossil fuels meet current energy demand, and it is estimated that they will supply at least 80% of the energy demand required in 2040 (Beltrán-Telles et al., 2017), we decided to use only SSP5 8.5 W/m 2 to model the availability of suitable climatic environments for the presence of the species under study. Likewise, we considered that SSP5 8.5 W/m 2 is the climate scenario that provides the best test of our hypothesis as this will result in the largest difference between current and future environmental conditions.

Vulnerability of climatic suitability in the face of environmental variations
The Maximum Entropy (MaxEnt) approach was used to model the climatic suitability of the 14 species of Crotalus. MaxEnt uses the principle of maximum entropy on presenceonly data to estimate a set of functions that relate environmental variables and climatic suitability to approximate the species' niche and potential geographic distribution (Phillips et al., 2017). Therefore, the species distribution model considered in this study represents a correlative species distribution model (Phillips, Anderson & Schapire, 2006), subject to the challenge of balancing goodness of fit with model complexity, as models that are inappropriately complex or simple show reduced ability to infer habitat quality, reduced ability to infer relative importance of the variables in the restriction of the distribution of the species, and a reduced transferability to other time periods (Warren & Seifert, 2011). In our case study, using the ''ENMeval'' library (Muscarella et al., 2014) in the software R, the calibration of the model for each species considered the choice of (a) accessible area (background or M area), (b) the type of variables that MaxEnt constructs (features), (c) regularization multiplier, and (d) the type of model output (raw, cumulative, logistic), as these considerations affect the inferences to be made (Fourcade et al., 2014).
Using the MaxEnt software, the information obtained from the calibrated models was projected within the known distribution area of the species under study. We used the layers of the current climate and those of the future climate mentioned above. All climatic layers were obtained from the WorldClim database v2.1 (https://www.worldclim.org/). The models were generated with a climatic suitability gradient from 0 (low suitability) to 1 (high suitability), which were then converted to binary models (presence/absence). For each species, the threshold Maximum training sensitivity plus specificity (MaxSS) provided by MaxEnt in each model was chosen. The threshold MaxSS has been reported to show good performance for models that work only with presence data (Liu, White & Newell, 2013). The importance of each bioclimatic variable in the observed distribution of the species under study was evaluated according to the relative importance of each variable, which was obtained by adding the percentage of contribution (PC) and the importance of permutation (IP), evaluated by MaxEnt, and the result was divided by two et al., 2015). As a last step, the climatic suitability of the realized niche of each species was measured under current and future climate conditions. The vulnerability of the climatic suitability of each species to climate change was also identified, using the following change rate analysis: % of change = (S1−S0) S0 * 100, where S0 is the total surface of the study area, according to the base scenario, and S1 is the total surface occupied in the study area under change conditions.

Climate profile
The PCA suggested that, for our study area, the climate profile could be explained by considering the first two components. In all cases between components one and two, they explained at least 95% of the variation in the data. Under current weather conditions, for example, component one explained 96.2% of this variation, while component two only explained 2.8%. Considering the climate change scenarios, the scenario that presented the value with the lowest percentage in the sum of the two components was the BCC-CSM2-MR scenario with a value of 95.1%. The highest value was presented in the CNRM-CM6-1, and IPSL-CM6A-LR scenarios with 96%.
Regarding the contribution of the variables for each component, for both the current climate conditions and the climate change scenarios, the variable Annual Precipitation was the one that presented the greatest contribution in component one. For component two, considering the current climate conditions and the climate change scenarios, the variables Precipitation of Warmest Quarter and Precipitation of Coldest Quarter were the ones that presented the greatest contribution (Table 1); however, the Analysis of Variance and Tukey post hoc tests suggested that only the climatic variables Annual Mean Temperature and Mean Temperature of Wettest Quarter presented significant statistical differences in their means with respect to the three climate change scenarios used in this study. While the variable Mean Diurnal Range only presented significant differences in its means when contrasted with the climatic information proposed for the scenarios BCC-CSM2-MR and CNRM-CM6-1, the rest of the variables did not present significant differences (Table 2).
Regarding the size of the niches, our results showed that these amplitudes varied among components. For example, C. ravus presented the greatest niche width considering the     Table 3).

Vulnerability of climatic suitability in the face of environmental variations
The models obtained for the species under study showed an area under the curve ranging from 0.80 to 0.95, indicating low levels of commission (predicts the presence of the species where it does not exist, false positive) and omission (predicts the non-presence of the species where it really exists, false negative) ( Table 4). The relative importance of each variable in the generation of climatically suitable zones for the presence of the species under study indicated that the variable Annual Mean Temperature presented a greater contribution for 42.8% of these species. The variables Annual Precipitation, Mean Temperature of Wettest Quarter, and Precipitation of Coldest Quarter presented a higher contribution for the 28.5%, 14.2%, and 7.14% of species under study, respectively. The rest of the variables did not present a marked influence on the generation of climatically suitable zones for the species under study (Table 4). The models allowed the identification of three groups of species according to the percentage of loss of climatic suitability between current climatic conditions and the three climate change scenarios considered in this work (Fig. S1). In the first group (high vulnerability), the species C. viridis, C. scutulatus, C. molossus, and C. ravus showed a loss of climatic suitability of between 40% and 66% in at least two climate change scenarios used in this study. In the second group (medium vulnerability), C. pricei, C. ruber, C. lepidus, C. basiliscus, C. tigris, and C. cerastes showed a loss of climatic suitability of between 1% and 34%. In group three (low vulnerability), the species C. willardi, C. intermedius, C. enyo, Table 4 The relative importance values of each variable in the generation of habitat suitability models for the Crotalus species under study. Area under the curve (AUC) values also provided that allow the evaluation of habitat suitability models. Variables are bio1 = Annual Mean Temperature, bio2 = Mean Diurnal Range, bio8 = Mean Temperature of Wettest Quarter, bio12 = Annual Precipitation, bio13 = Precipitation of Wettest Month, bio14 = Precipitation of Driest Month, bio15 = Precipitation Seasonality, bio18 = Precipitation of Warmest Quarter, bio19 = Precipitation of Coldest Quarter. and C. atrox showed an increase in climatic suitability for the climate change scenarios considered in this study (Table 5). Hutchinson (1957) defines the niche of a species as a n-dimensional space, where each dimension represents the response of a species to the variation of a certain variable. In this way, each site on Earth is characterized by a set of environmental conditions that define a specific habitat inhabited or uninhabited by a community of species (Kearney, 2006). In this sense, our results indicate that for current climate conditions, according to the PCA, the climatic profile of the distribution area of the species under study can be viewed from two approaches. The first is approach one (PC1), where the climate profile is determined to a greater extent by the Annual Precipitation. With approach two (PC2), the greatest contribution is provided by the variables Precipitation of Warmest Quarter and Precipitation of Coldest Quarter. For the climate change scenarios used in this study, the variables Annual Precipitation, Precipitation of Warmest Quarter, and Precipitation of Coldest Quarter will continue to make the greatest contribution to the climate profile. Climate change in the last 30 years has produced numerous changes in the distribution and abundance of species (Parmesan & Yohe, 2003;Root et al., 2003) and has been implicated in the extinction of several species (Pounds, Fogden & Campbell, 1999). For the period 2021-2040, our results of climatic suitability loss identify three levels of vulnerability (high, medium, and low) for the species under study. For the group with high vulnerability, we identified C. viridis, C. scutulatus, C. molossus, and C. ravus, which represents 28.5% of the analyzed species. In the group of medium vulnerability, we identified C. pricei, C. ruber, C. lepidus, C. basiliscus, C. tigris, and C. cerastes, which represent 42.8% of our studied species. The species with low vulnerability includes C. willardi, C. intermedius, C. enyo, and C. atrox, representing 28.6% of our considered species. Various authors have pointed out that the breadth of the niche can have an important effect on the risk of extinction of a species because species with broader niches could be less vulnerable to abrupt environmental variation under anthropogenic climate change. At the opposite extreme, species with narrow niches would be particularly threatened by climatic changes (Brown, 1984;Johnson, 1998;Kotiaho et al., 2005;Pearson et al., 2014;Saupe et al., 2015). There is substantial evidence from a variety of taxa that supports the theory that narrowed niches drive the risk of extinction of species in the face of climate change variations (e.g., fishes (Munday, 2004), bats (Boyles & Storm, 2007, birds (Seoane & Carrascal, 2008), frogs (Botts, Erasmus & Alexander, 2013), and plants (Ozinga et al., 2013). In relation to this, for the period 2021-2040, within the high-vulnerability and medium-vulnerability groups, C. viridis, C. molossus, C. tigris, C. scutulatus, C. ruber, and C. cerastes showed reduced niches for the variables related to temperature. This coincides with the aforementioned predictions since it would be expected that the species under study with reduced niches related to temperature present a greater disturbance in their habitat with respect to the increase in temperature projected for the period 2021-2040. However, other species in these same two groups (C. ravus, C. basiliscus, and C. lepidus) show a greater niche width compared to several species classified in the low-vulnerabilty group (C. atrox, C. pricei, and C. intermedius). This finding contrasts what is proposed above. In this context, Carrillo-Angeles et al. (2016) suggest that although various studies reinforce the hypothesis that species with narrow niches are more susceptible to climate change, there is no single trend in the fate of species with narrow niches and their vulnerability to environmental variations. For example, projections for an increase in greenhouse gases and, consequently, in temperature, for the year 2050 in Europe suggest that some of the most affected species will be those that inhabit colder northern regions, species with low densities, and species with less tolerance to aridity (Huntley et al., 1995;Thuiller et al., 2005a).

DISCUSSION
Related to this last point, evidence suggests an increase in temperature and low rainfall for the period 2021-2040. For example, the comparison of means indicates that the variables Annual Precipitation, Precipitation of Warmest Quarter, and Precipitation of Coldest Quarter will present a relative stability for the period 2021-2040, with respect to what is shown in the current climate. However, for the variables Annual Mean Temperature, Mean Diurnal Range, Mean Temperature of Wettest Quarter, an increase in the averages of between 1.74 • C and 1.99 • C is expected; 0.11 • C and 0.49 • C, and 1.1 • C and 1.8 • C, respectively. In this regard, various studies have mentioned that the significant increase in temperature and the low availability of water will lead to a reduction in humidity of the air and substrate (Seager et al., 2007;Ye & Grimm, 2013;Kunkel et al., 2013). This is a condition that may have significant detrimental effects on reptiles that are less tolerant to aridity (Inman et al., 2014;Hatten et al., 2016).
Our results show that for C. ravus, C. basiliscus, and C. lepidus, despite presenting wide climatic niches for the variables related to precipitation and temperature, their ideal habitat is influenced to a greater extent the Annual Mean Temperature and Mean Temperature of Wettest Quarter, respectively. Like the rest of the species classified as high and medium vulnerability, they are also influenced to a greater extent by the variables Annual Mean Temperature and Mean Temperature of Wettest Quarter. In contrast, for C. atrox, C. enyo, C. willardi, and C. intermedius, four species identified with low vulnerability to climate change, the variables related to temperature show little contribution to the generation of suitable climatic environments for their distribution. In this way, the evidence suggests that for our species identified with high vulnerability to climate change, they can be considered as less tolerant to the increase in aridity projected for the period 2021-2040.
In conclusion, the increase in the variables Annual Mean Temperature and Mean Temperature of Wettest Quarter may compromise the climatic suitability of at least 71.4% of the species considered in our study. In this sense, for the species under study, the niche width, by itself, cannot be considered as a determining factor that helps to predict the vulnerability of their climatic suitability under rapid environmental change. However, evidence from our study shows how the relative importance of climatic variables in the construction of niche modeling can help us understand the vulnerability of the climatic suitability of the species under study to global climate change.
In this study, we used correlative methods to model the climatic suitability of the species under study and estimate niche width. Soberón (2007) pointed out that the realized niche is determined by biotic restrictions in the fundamental ecophysiological niche, population dynamics (e.g., source-sink dynamics) and dispersion limitations (that is, accessibility). Therefore, in our study, we are not considering the physiological limits of the species and, although Cuervo-Robayo et al. (2017) comment that correlative ecological niche models are a good technique to capture exposure to climate change, we cannot rule out that we could be underestimating or overestimating our results. However, mechanistic (physiological) methods can also be subject to overestimation or underestimation of the niche (Peterson & Holt, 2003;Strubbe et al., 2015).