Bighorn sheep (Ovis canadensis) have a wide geographic range, spanning western Canada, the western United States, and northern Mexico (Figure 1) (Brewer et al., 2014). Given the diversity of landscapes and environmental pressures across this range, management strategies must be adapted to address region-specific conservation challenges.

ENMs need two types of data to estimate potential distribution: a set of georeferenced locations where the species has been detected and a set of quantitative raster maps describing the environmental conditions under which the target species lives. One of the most-used algorithms for ENMs is MaxEnt (Maximum Entropy) (Phillips et al., 2006), which estimates the potential distribution index of a species by finding the most spread-out, or “maximum entropy,” distribution given environmental constraints. MaxEnt is robust with small sample sizes and performs well when only presence data are available, making it a widely recognized tool for species distribution modeling. It is used for estimating the potential distribution index and potential range shifts under changing environmental conditions (Elith and Leathwick, 2009; Phillips et al., 2006; Merow et al., 2013). This study uses climatic variables that are among the most important factors driving species’ distribution (Grinnell, 1917; Guisan et al., 2013), especially at large spatial extents, as they have a direct influence on organisms’ behavior and physiology. They are particularly important for plants, which cannot evade adverse weather by sheltering or migrating.

Historical occurrence data for the bighorn sheep subspecies were obtained from an open-access repository, the Global Biodiversity Information Facility (GBIF.org) web database. GBIF collects biodiversity data by aggregating species occurrence records from global institutions, researchers, and citizen scientists, using standardized formats. Occurrence data were filtered from 1994 to 2020, and only the subspecies O. c. mexicana, O. c. cremnobatres, and O. c. weemsi were included. The 1,158 records underwent a spatial verification process, where data points carefully were reviewed to ensure their accuracy and consistency with the species’ historical known range. This involved identifying and removing any erroneous or mislocated entries from the database. After this quality control, a final set of 280 georeferenced records remained, providing a reliable foundation for analysis. The calibration area, also referred to as the M region (Barve et al., 2011), defines the accessible area for species, crucial for accurate niche modeling by setting boundaries for potential distribution predictions and minimizing bias. For these subspecies in Mexico, the calibration area was based on the ecoregions in which they are distributed (Escobar-Flores et al., 2015), as these units represent geographic barriers for several species’ dispersal. The selected ecoregions for O. canadensis in Mexico included the Californian, Baja Californian, Sonoran, Altiplano Norte (also known as the Chihuahuan Desert), Tamaulipas, and the Mexican transition zone Sierra Madre Occidental (Figure 2).

Nineteen bioclimatic layers representing annual, seasonal, and extreme climatic patterns were downloaded from the Chelsa database (https://chelsa-climate.org) (Brun et al., 2022), providing high-resolution climate data that offers a more accurate fit compared to simple monthly or yearly averages. The datasets in Karger et al. (2018) cover both the current period (1979-2013) and future projections (2014-2060). All layers were downloaded at a 30 arc-second resolution and clipped to the calibration area. Future GCMs were incorporated, based on scenarios of GHG emissions or representative concentration pathways (RCPs), to generate climate change projections. The layers were cropped to the M region, which defines the species’ accessible area, enhancing the models’ precision. Uncertainty in the climate projections was addressed by using four different RCPs (2.6, 4.5, 6.0, 8.5) alongside six GCMs (Supplementary Appendix 1). We developed a GCM consensus map (Figure 3) by averaging the outputs of multiple GCMs over 20-year horizon (2041–2060) to capture a robust range of climate projections. The RCPs represent a range of potential GHG emission trajectories, from the low-emissions pathway (RCP 2.6), which assumes emissions peak and decline shortly after 2020, to RCP 8.5, which assumes continued growth emissions through 2,100 (Weyant et al., 1996). Although RCP 8.5 has been widely used in ecological modeling to represent worst-case climate outcomes, recent literature suggests it no longer reflects a plausible “business-as-usual” scenario. Its underlying assumptions —such as high fossil fuel dependency and negligible mitigation— are increasingly inconsistent with global policy and energy trends (Peters and Hausfather, 2020). The IPCC Sixth Assessment Report (2023) similarly notes that recent mitigation efforts and the declining cost of low-emissions technologies are driving reductions in global energy and carbon intensity, making extreme high-emission pathways less likely. Nevertheless, we retain RCP 8.5 in this analysis to explore upper-bound climate risks and to apply the precautionary principle in conservation planning (Cooney, 2004; Rye et al., 2021), ensuring that even low-probability but high-impact outcomes are considered.

After inputting the data, we followed a structured approach to process and refine the information. The first step involved preselecting the environmental variables to ensure model accuracy. To avoid redundancy, we assessed multicollinearity among environmental variables using the Pearson correlation coefficient in R (version 4.3.3, R Core Team, 2024). The correlation analysis was conducted within the R base environment, allowing us to identify and manage highly correlated variables.

When two variables had a Pearson correlation coefficient |r| > 0.85 (indicating strong correlation) (Pradhan, 2016; Schober and Schwarte, 2018), only one was retained for model development. Selection criteria were based on each variable’s predictive strength and expert knowledge to ensure the most relevant covariates were used. This analysis determined that 10 environmental variables best explained bighorn sheep presence (Supplementary Appendix 2).

The predicted model was initially calibrated for the period 1979–2013 and then projected for 2041–2060, incorporating various Representative Concentration Pathways (RCPs) and General Circulation Models (GCMs) to account for different climate scenarios. The final dataset was analyzed using Tableau (version 2018.3), an interactive data visualization software, to better interpret spatial patterns and model outcomes.

Uncertainty in climate change projections was addressed by averaging outputs from six GCMs and four RCPs. This approach reduces model-specific biases and integrates a broad spectrum of greenhouse gas emission scenarios, from low to high. By combining multiple models and scenarios, the analysis captures both structural and scenario-related uncertainties, resulting in more reliable climate impact estimates.

This improved reliability strengthens conservation planning and habitat management by providing robust, data-driven insights. To apply these findings at a local scale, we incorporated a raster layer of ejido polygons in Mexico. Using ejidos in Baja California, Baja California Sur, and Sonora as a reference, we identified highly suitable areas for bighorn sheep. Identifying these areas is crucial for policy intervention, as ejidos play a key role in habitat management, conservation efforts, and species protection under changing environmental conditions.

The MaxEnt model was run with five replicates, using 25% of the data for testing, and a maximum of 8,000 background points. The model was set to run for up to 500 iterations with a regularization parameter of 0.5. To estimate potential distribution, we used a probability transformation method that accounts for species presence likelihood in a nonlinear way (log–log output format). The random test percentage was 25 percent, which means that 75 percent of the total database was used as the random sample to train the model, and the other 25 percent was used to test the model predictions. To avoid overfitting the test data, we set the regularization multiplier value to 0.5. The output format was a default cloglog transform, which has a stronger theoretical justification than the logistic transform (Phillips et al., 2017). This study used linear, quadratic, product, threshold, and hinge feature classes (auto features) based on sample sizes (Phillips and Dudík, 2008). Phillips and Dudík (2008) classified the features for sample sizes: an auto features setting for more than 80 records, quadratic and hinge setting for 15–79 records, linear and quadratic for 10–14 records, and linear setting for sample sizes fewer than 10 records.

The ENM established a relationship between the bighorn sheep occurrence and predictor variables and estimates the present and future climatic niche for Mexico’s bighorn sheep. A potential distribution map for Ovis canadensis was produced by utilizing the area under the threshold-independent receiver operating characteristic curve (AUC) weight averages of the five log–log output format maps produced by five-fold cross-validation. The goodness-of-fit test of the MaxEnt model was evaluated by AUC in which the relative suitability ranged from 0 to 1, among which that of 0.5 suggests that the models show no predicting capability, while that of >0.7 represents that the models are acceptable. The model also was evaluated by the true skill statistic (TSS) that measures model performance, comparing predicted and actual species occurrences (Allouche et al., 2006) based on five-fold cross-validation. The TSS ranges between −1 and +1, where +1 indicates perfect agreement and values of zero or less indicate a performance no better than random (Zhang et al., 2020). The following ranges were used to interpret TSS statistics: values < 0.4 were poor, 0.4–0.8 useful, and > 0.8 were good to excellent (Zhang et al., 2015).

The average AUC value of niche models from five-fold cross-validation was 0.9887, indicating high discrimination accuracy between presence and background locations (Allouche et al., 2006). The True Skill Statistic (TSS) value was 0.4955, suggesting a useful and reliable model for predicting bighorn sheep distribution.

Precipitation and temperature are well-documented key climate variables influencing species distributions (Zhong et al., 2010). For bighorn sheep, three bioclimatic variables had the highest positive contributions to their probability of occurrence:

Extreme temperatures can constrain bighorn sheep distribution limits, affecting their ability to persist in certain regions and influencing distribution shifts over time.

Figure 4 illustrates that the current potential distribution of bighorn sheep in Mexico is primarily restricted to the northwest, including Baja California, Baja California Sur, and Sonora, as well as Tiburón Island. The map displays only areas with high potential distribution index (0.4 to 0.99) under current climate conditions. Baja California (Figure 4, number 1) has the highest number of ejidos with high potential distribution index (n = 35). Baja California Sur (Figure 4, number 2) has the highest number of ejidos (12) with high potential distribution index (0.6 and even close to 0.9). Sonora (Figure 4, number 3) is the third state with a high potential distribution index in five ejidos, including Tiburon Island.

We compared the potential distribution in all polygons. Figure 5 shows high potential distribution (30–99 percent) areas in the ejido polygons. The potential distribution estimated in each GCM differed. In general, the bighorn sheep potential distribution decreases in the future, especially in the extreme RCP 8.5. Notably, climatic suitability declines in key conservation areas, including Baja California, ejido Bonfil, and Tiburón Island, not only under the high-emission scenario (RCP 8.5) but also under RCP 6.0 and RCP 4.5 across nearly all General Circulation Models (GCMs). The only exception is the Earth System Model (ESM) from the College of Global Change and Earth System Science, Beijing Normal University, which projects more stable climatic conditions in these regions (Figure 5).

Figure 3, column 1, shows that ejido polygons in Baja California (BC), Baja California Sur (BCS), and Sonora have an average potential distribution index of 0.85 under current conditions. However, when considering future climate scenarios, the potential distribution index changes, with a reduction in areas offering favorable precipitation and temperature conditions, especially under RCP 8.5.

Figure 3, column (BC, BCS, Sonora) shows that under RCP 8.5, the potential distribution for bighorn sheep populations decreases, with favorable areas predominantly shifting to the northern part of Baja California. The study shows significant changes in the potential distribution for bighorn sheep across future climate projections. Figure 3, column 2 (Baja California Sur), reveals a decline in potential distribution in the southern polygons of Baja California Sur (with an average potential distribution of 0.57, so a 32 percent decrease) and Sonora (with an average potential distribution index of 0.37, so a 56 percent decrease). However, under the more extreme RCP 8.5 scenario, some polygons in Baja California show an increase in potential distribution. It also is observed that the potential distribution index within the polygon boundaries of ejido Bonfil decreases by 38 percent under climate change scenarios from 0.87 to 0.55 (Figure 3, column Ejido Bonfil potential distribution).

MaxEnt outcomes give estimates of relative contributions of the environmental variables to the model (Table 1). Accordingly, Bios17—Precipitation of Driest Quarter and Bios12—Annual Precipitation are the variables contributing the most to potential distribution predictions of Ovis canadensis in Mexico.

Figure 6 suggests that negative changes to climatic suitability are more drastic at the Sonoran and Baja Californian ecoregion, where more than half of the ejidos present a negative percentage change with respect to the current potential distribution index, while ejidos in the Californian ecoregion seem to have a negative percentage change in at least half of its ejidos. Thus, the Californian ecoregion is where ejidos have a lower decrease in the bighorn sheep potential distribution index.

Climate change poses a significant threat to bighorn sheep habitat in Mexico, with projections indicating a reduction in suitable areas under future climate scenarios. Under the high-emission RCP 8.5 scenario, the species’ potential distribution is expected to decline, highlighting the substantial impact that climate change may have on habitat suitability. However, model results also reveal a geographic shift in favorable conditions, with some polygons in northern Baja California showing increased potential distribution under this more extreme scenario.

This seemingly paradoxical outcome may be explained by a northward shift in climatic conditions that more closely align with the species’ bioclimatic niche—particularly in areas experiencing a relative decrease in seasonal water stress—. These localized increases in potential distribution highlight the complex, non-linear, and spatially heterogeneous nature of climate impacts. They also reinforce the need for geographically targeted and adaptive conservation strategies. One of the most concerning projections is the 38% decline in potential distribution within ejido Bonfil, where the potential distribution index drops from 0.87 to 0.55 under future climate change scenarios. Given that this area has been a key conservation stronghold for bighorn sheep, the projected reduction suggests that current management efforts may be insufficient to ensure the species’ long-term viability in this location. This underscores the urgency of integrating climate adaptation measures into conservation planning, particularly in areas projected to undergo the most severe declines in the species’ occurrence.

In this context, mapping potential future distributions becomes critical for guiding species management. Identifying and prioritizing high-potential distribution areas enables conservation managers to focus resources on the landscapes most likely to support bighorn sheep in the future. Moreover, the use of multiple GCMs highlights uncertainty in climate projections, yet consensus maps, such as the one developed in this study, serve as essential tools for policy-oriented wildlife management decisions.

Results from the MaxEnt model indicate that Bio17 (Precipitation of the Driest Quarter) and Bio12 (Annual Precipitation) are the most influential environmental variables (Supplementary Appendix 3). These findings reinforce the critical role of precipitation. In arid and semi-arid ecosystems, such stress is a key limiting factor, and for bighorn sheep, the driest part of the year imposes critical physiological and ecological challenges (Turner and Weaver, 1980). These variables capture the minimum precipitation during this period, influencing forage availability, hydration, and movement patterns. Incorporating these variables into the model allows the model to assess potential distribution under conditions of seasonal drought, which is projected to intensify in much of the species’ range.

Given the variability in potential distribution across different GCMs, it is essential to integrate multiple computational approaches to address uncertainties in wildlife management planning. RCP 8.5 shows significant negative changes in temperature and precipitation, which could adversely affect species distributions. However, some northern Baja California polygons remain highly suitable, suggesting that shifts in conservation priorities may be necessary. While RCP 8.5 is increasingly considered less likely under current emissions trajectories (IPCC, 2023; Peters and Hausfather, 2020), its inclusion helps identify areas at high risk and stress-test management strategies under extreme climate futures. However, some northern Baja California polygons remain highly suitable, suggesting that shifts in conservation priorities may be necessary. Decision-makers must evaluate all possible futures by considering multiple GCMs and RCPs (Lempert et al., 2006), ensuring that management strategies remain flexible and responsive to new climate data.

A key conservation strategy involves focusing on areas with a high potential distribution index, particularly in northern Baja California under RCP 8.5. Prioritizing these regions for habitat protection, restoration, and monitoring is essential to maintaining suitable environments for bighorn sheep as climate change impacts their current range. Understanding which populations are under the most climate-related stress is crucial for anticipating future conservation and management actions (Zamora-Maldonado et al., 2021).

The Bonfil ejido, with 28 years of experience (1996–2024) in bighorn sheep conservation and sustainable hunting programs, has developed extensive expertise in species management. This accumulated knowledge could serve as a foundation for establishing a training and collaboration program to assist other ejidos in adapting to future habitat shifts. As conservation efforts may need to shift toward northern Baja California, ejido-based management strategies should be expanded and restructured to align with these geographic changes. Importantly, ejido owners are already engaged in ongoing dialog to promote collaborative and sustainable decision-making regarding bighorn sheep management.

Effective wildlife management must balance conservation priorities with economic incentives. Large herbivores such as bighorn sheep play a key role in ecosystem dynamics, influencing primary production, nutrient cycling, and habitat structure (Danell, 2006). However, they also provide economic benefits through regulated hunting programs, which generate revenue for local communities and conservation initiatives.

Assessing the impacts of climate change on economically valuable species is critical for biodiversity stewardship and local livelihoods (Saba et al., 2012). While much research has focused on fisheries and climate change (Islam et al., 2014; Sumaila et al., 2011), there is a need to further quantify how climate change will affect terrestrial wildlife species with economic significance, such as bighorn sheep (Advani, 2014).

While this study provides valuable insights into the potential effects of climate change on bighorn sheep potential distribution, we recognize its limitations. Landscape features such as slope, terrain ruggedness, and distance to water are undeniably critical for bighorn sheep at finer spatial scales. However, this study is not intended to replace localized assessments but rather to serve as a complementary tool for understanding broad-scale climate-driven habitat shifts.

Our decision to focus exclusively on bioclimatic variables is based on the study’s primary objective: to assess how climate change influences potential bighorn sheep distributions. Bioclimatic factors are particularly suited for this purpose because they:

Directly and indirectly shape species distributions – Climate affects physiological tolerances (Epps et al., 2004) and influences vegetation patterns and water availability (Cruz et al., 2024), both of which are crucial for bighorn sheep survival.

Enable long-term projections – Unlike terrain variables, which remain largely unchanged over time, bioclimatic variables allow us to model habitat shifts under future climate scenarios (e.g., RCPs 2.6, 4.5, 6.0, 8.5) using multiple General Circulation Models (GCMs).

We acknowledge that integrating additional ecological and landscape variables could further estimate habitat suitability. However, the results of this study offer a valuable foundation for broad-scale conservation planning, helping to identify priority areas for further, more localized research. Future studies can build upon this analysis by incorporating landscape features and species movement patterns, ensuring a more comprehensive approach to bighorn sheep conservation under changing climate conditions.

This study underscores the importance of using ENM as a strategic tool for conservation planning. Future research should aim to deepen our understanding of how species interactions, land-use changes, and socioeconomic factors influence potential distribution. In this regard, we propose a series of proactive strategies to improve bighorn sheep conservation and management under future climate conditions: