Rangewide habitat suitability analysis for the Mexican wolf (Canis lupus baileyi) to identify recovery areas in its historical distribution

To develop an updated distribution model and habitat suitability analysis for the Mexican wolf, to inform the recovery efforts in Mexico and the United States.


| INTRODUC TI ON
The Mexican wolf, Canis lupus baileyi (Nelson & Goldman, 1929), the smallest and genetically most distinct grey wolf subspecies in North America (Nowak, 1995), formerly inhabited the temperate forests of Mexico and the southwestern United States before being nearly extirpated during the 20th century (Gish, 1977;McBride, 1980). In 1976, the U.S. Fish and Wildlife Service (USFWS) listed the Mexican wolf as an endangered species (USFWS, 1976). At that time, there were fewer than 50 individuals in the wild, all in the Sierra Madre Occidental, Mexico (Brown, 1983).
In 1998, after a successful captive breeding programme (Siminski, 2016), the USFWS and its partners released the first wolves from captivity into central parts of Arizona and New Mexico (currently designated as the Mexican Wolf Experimental Population Area or MWEPA), to establish a self-sustaining wild population (USFWS, 1982). Currently, this area has a wild population of at least 163 individuals (USFWS, 2020). In the early 1980s, recovery efforts began in Mexico with initiatives to identify appropriate sites to establish a Mexican wolf population (CONANP, 2009) In the last 20 years, several efforts, using information available at the time, have been undertaken to identify suitable areas for the recovery of the Mexican wolf. Most of these efforts have focused on specific regions of the United States and Mexico (Araiza, 2001;Araiza et al., 2012;Carnes, 2011;Carroll et al., 2003Carroll et al., , 2006Carroll et al., , 2014Carroll et al., 2004;Martínez-Gutiérrez, 2007). One of the challenges to complete past rangewide habitat analyses has been the availability of reliable information spanning both countries, particularly on the availability of wild ungulate prey (Carroll et al., 2004).
However, in recent years, regional and global databases have been available to help overcome this limitation. To our knowledge, only one study (Hendricks et al., 2016) has attempted a rangewide analysis to redefine the historical distribution of the Mexican wolf and to identify suitable areas in its expanded range to guide recovery.
Unfortunately, their analysis contains critical methodological flaws for delineating suitable areas that render their results and conclusions questionable for informing recovery planning (Heffelfinger et al., 2017a).
Here, we present the modelled distribution and habitat suitability analysis for the Mexican wolf that was developed for, and integrated into, the Revised Mexican Wolf Recovery Plan (USFWS, 2017). This analysis used a niche-based species distribution modelling approach and incorporated natural and anthropogenic factors, derived from global databases, that previously were identified as influencing wolf population establishment and persistence, including land cover and vegetation, human population density and road density (Carroll et al., 2014;Jedrzejewski et al., 2004;Oakleaf et al., 2006). Working with the agencies responsible for the management and monitoring of wild ungulate populations, we included information on relative wild prey density for the first time, one of the key factors influencing wolf population success (Fuller et al., 1992(Fuller et al., , 2003.

| ME THODS
We carried out the analysis in four steps: (1) modelled the distribution of the grey wolf in the southwestern United States and Mexico using an ensemble niche-based modelling approach; (2) compiled and standardized ecological and anthropogenic habitat variables for Mexico and the United States; (3) modelled the habitat suitability for the Mexican wolf throughout the modelled distribution; and (4) quantified the largest, continuous high-quality habitat patches across its modelled distribution. We describe each step below and include details in the Appendix S1, Detailed methodology.

| Modelling the distribution of the grey wolf in Mexico and the southwestern United States
The delineation of the historical distribution of Mexican wolves has been controversial during the last two decades of recovery planning (Heffelfinger et al., 2017a;Heffelfinger et al., 2017b;Hendricks et al., 2016Hendricks et al., , 2019Leonard et al., 2005). For our analysis, we designated the extent of the modelling area ("M" sensu Barve et al., 2011) from central Utah and Colorado in the United States to southern Mexico ( Figure 1). In this way, we included an area beyond the historical distribution of the subspecies and the expanded distribution recognized by Parsons (1996) adopted by the USFWS.
To model the distribution of the grey wolf in the southwestern United States and Mexico, we followed a niche-based species distribution modelling (SDM) approximation. Niche modelling algorithms characterize the ecological niche (sensu Hutchinson, 1957) of a taxon by looking for non-random associations between a collection of occurrence records and environmental conditions of the region where this taxon occurs. Then, similar conditions across the study region are identified and potential distribution is mapped (Guisan & Thuiller, 2005;Guisan et al., 2017;Pearson & Dawson, 2003;Peterson et al., 2011).

| Input data
We established an area ("operational area"; Figure 1) that not only included the geographic delineation of the Mexican wolf historical range proposed by Goldman (1944), Hall (1981) and Nowak (1995), but also the expanded delineation used by the USFWS (Parsons, 1996), based on Bogan and Mehlhop (1983), who incorporated C. l. mogollonensis and C. l. monstrabilis into C. l. baileyi. Therefore, we included the totality of records of specimens recognized as Canis lupus baileyi by all the authors, but also the records within the operational area of specimens classified as C. l. baileyi, C. l. mogollensis, C. l. nubilus, C. l. youngi or C. l. monstrabilis, depending on the author. We acknowledge that our operational area includes areas of intermixing between the Mexican wolf and its neighbouring subspecies, but because range limits are not rigid and our knowledge of them is incomplete, any delineation will be somewhat imprecise.
We gathered grey wolf occurrences within this operational area from the literature (Araiza et al., 2012;Brown, 1983;Goldman, 1944;Hall, 1981;Martínez-Meyer et al., 2006;Nowak, 1995), electronic databases (i.e., GBIF, VertNet) and oral records from local trappers (Brown, 1983;Servín et al., 2007), from 1848 to 1980. We excluded existing records for the Mexican wolf after 1980 because these correspond exclusively to individuals reintroduced or born in the areas designated for its recovery in the United States and Mexico.
We reviewed each record and discarded those for which location description was too ambiguous for accurate georeferencing. Next, we divided the records according to their reliability into primary (specimens preserved in natural history collections) and secondary (observations or interviews with no physical evidence). We used only the primary records to calibrate ecological niche models and secondary records for model validation. Due to the strong clustering observed in the primary records in some areas, we filtered occurrences at 25 km separation distance, that is, we eliminated records that were <25 km apart, aiming to reduce clustering and thus over-representing in the model the environmental conditions of such areas (Boria et al., 2014). To do so, we used the thin function in the spThin R package ( Aiello-Lammens et al., 2015). Our final data set consisted of 41 primary and 296 secondary spatially unique occurrences ( Figure 1).
We used 19 bioclimatic variables obtained from the WorldClim database (Hijmans, Cameron, et al., 2005; Table S1) and three topographic variables: elevation, slope and topographic heterogeneity from the HYDRO1k database (USGS, 2001). To avoid model overfitting, we used only the most informative variables (>1% contribution) selected via the permutation method implemented in the Maxent program Searcy & Shaffer, 2016). The final set included 15 variables at 30 arc seconds (approx. 1 km 2 ) spatial resolution (Table S1).

| Ecological niche and distribution modelling
Considering that no single niche modelling algorithm performs better under any condition (Qiao et al., 2015), we tested eight algorithms and support vector machine (SVM). We ran all algorithms in the R packages sdm (Naimi & Araújo, 2016) and dismo (Hijmans, Phillips, et al., 2005), except Maxent, for which we used its stand-alone interface ; see Detailed methodology in Appendix S1).
All models had a final spatial resolution of 30 arc seconds. F I G U R E 1 Occurrence records used for the construction of niche models. Primary records (yellow dots) were used for calibration and secondary records (black dots) for validation. The dashed orange area was established to obtain the gray wolf occurrences and was delineated based on the taxonomic classification for Canis lupus baileyi by Goldman (1944), Hall (1981) and Nowak (1995), and the extended range proposed by Parsons (1996), based on Bogan and Mehlhop (1983). The dashed square represents the extent of analysis ("M" sensu Barve et al., 2011) We converted the resulting continuous maps into binary (presence-absence, Figure S1) based on a 10-percentile threshold value to account for undetected uncertainty/error in occurrences (Liu et al., 2005), and validated each model using a combination of four metrics: omission and commission errors (i.e., the number of presences predicted as absences and vice versa), true skill statistic (TSS) and a chi-square test (Allouche et al., 2006;Fielding & Bell, 1997).
Then, we generated a consensus map with the four algorithms that performed best by summing each binary map, selected the areas where two or more models coincided and converted that into a binary map, resulting in the modelled distribution map (i.e., areas where suitable climatic conditions for the Mexican wolf exist, Figure S2). We approximated the historical distribution of Mexican wolves from the modelled distribution map ( Figure S2) by clipping off climatically suitable areas in biogeographic provinces with no historical occurrence records of the subspecies, that is Mediterranean California (Comer et al., 2003) and Baja California (Morrone, 2005), assuming that such biogeographic boundaries have represented barriers to the dispersal or establishment.
Finally, we characterized climatic suitability across our modelled distribution of Mexican wolves based on the notion that optimal conditions for a species prevail towards the ecological centroid of its multidimensional niche (Hutchinson, 1957;Maguire, 1973). We estimated the distance to the ecological niche centroid as a biologically meaningful measure of environmental suitability (Manthey et al., 2015;Martínez-Meyer et al., 2013;Osorio-Olvera et al., 2019Yañez-Arenas et al., 2012). To do so, we calculated the multivariate Euclidean distance of each pixel to the multidimensional mean in environmental space (i.e., the niche centroid). Then, we

| Environmental and anthropogenic habitat variables
We considered the following natural variables for the habitat analysis: (1) the climatic suitability score described above; (2) land cover and vegetation types; and (3) ungulate biomass. The anthropogenic variables considered were as follows: (1) human population density and (2) road density. All variables were clipped to our modelled distribution map of the Mexican wolf ( Figure 2a) and resampled from their native spatial resolution to 30 arc seconds (approx. 1 km 2 ) to have the same extent and spatial resolution for further analysis.

| Land cover and vegetation types
We used land cover information for the entire study region generated by the European Space Agency ( CCI-LC, 2015). This map represents the major land cover and vegetation types of the world produced for 2015 at a spatial resolution of 300 m. We performed a use/availability analysis for the Mexican wolf via a chi-square test (Table S2). Finally, the land cover layer was standardized to values from −1 to 1 based on the proportional occurrence of wolf records in the different land cover classes in the raster calculator of ArcGIS 10.2 (ESRI, 2014; see Appendix S1).

| Human population density
We used the Gridded Population of the World, ver. 4 (GPWv4) raster map (CIESIN- FAO-CIAT, 2005) Mladenoff et al. (1995), we established a threshold value of 1.52 humans/km 2 for the pessimistic scenario, so we rescaled pixel values below this density from 0 to 1 (with 1 being a human population density = 0), and we rescaled values above this threshold from 0 to −1 (where −1 was the maximum population density in the region). We calculated 1 and 2 SE above the pessimistic threshold resulting in a human population density of 3.13 humans/km 2 and 4.74 humans/ km 2 , which corresponded to the threshold values of the intermediate and optimistic scenarios, respectively. Then, these layers were also rescaled from −1 to 1.

| Road density
We used two data sources for roads: OpenStreetMap (http://www. opens treet map.org/), downloaded from Geofabrik (http://downl oad.geofa brik.de/), which is a vector map of the roads of the world at a maximum scale of 1:1,000 in urban areas, and we complemented this information with a road map for Mexico at a scale of 1:250,000 (INEGI, 2015). From these two maps, we selected paved roads and dirt roads suitable for two-wheel-drive vehicles and combined them into one consistent layer. From the unified, binational map we calculated road density (linear km/km 2 ) using the Line Density function in ArcGIS 10.2. Previous studies have found that wolves can persist in human-dominated landscapes with varying road density thresholds, ranging from 0.15 to 0.74 km/km 2 before wolf populations are affected (Fuller et al., 1992;Mladenoff et al., 1995Mladenoff et al., , 2009Sazatornil et al., 2016;Thiel, 1985;Vickery et al., 2001). Road density values were rescaled from −1 to 1 in the same way as we did with the human density map to construct pessimistic, optimistic, and intermediate scenarios ( Figure S4) using the following threshold values: 0.74 km/ km 2 for the optimistic, 0.15 km/km 2 for the pessimistic and the average of these two values, 0.445 km/km 2 , for the intermediate.

| Ungulate density
We used ungulate density estimations in the United States and Mexico to calculate an Ungulate Biomass Index (UBI) across wolf modelled distribution (Table S3). The UBI is a standardized value that uses a weighting factor based on mean animal mass to compare relative biomass available across different predator-prey systems (Fuller et al., 2003). Ungulate density estimates in the  García-Chávez, et al., 2014;López-González et al., 2012 and from 193 Unidades de Manejo para la Conservación de la Vida Silvestre (UMAs) in four states of Mexico: Sonora, Chihuahua, Durango, and Sinaloa, from 1999 to 2010 (Servín et al., 2008(Servín et al., , 2018. We also used cameratrap surveys for mule deer density estimates in northern Mexico . After preliminary analyses to model the UBI across the study area, we made several adjustments to improve the accuracy of our predictions for each species (see Appendix S1). We analysed rangewide density estimations for the three ungulate species using a generalized linear model (GLM) and random forest (RF) modelling to establish the best parameters to estimate UBI. For calibrating the model, we used 15 climatic, topographic and ecological variables (selected from an initial set of 27 based on their levels of significance; Table S4). We measured the reliability of individual species' models with a Pearson correlation analysis (r 2 ) and Akaike's information criterion (AIC). We calculated the UBI within the modelled distribution considering the ungulate abundance/density information and a weighting factor (β) for each species (see Appendix S1). We built the UBI distribution maps of each species across the whole study area in QGIS (QGIS Development Team, 2016) using the best fit GLM/RF models. Then, the UBI map of each species was clipped to its known distribution using the IUCN polygon maps (IUCN, 2016; Figure S5).
The three individual UBI maps were summed together in QGIS to produce a combined UBI map ( Figure S6), which was then clipped to our modelled distribution and finally rescaled from 0 to 1 to match the other layers for the habitat suitability model. This map represented the first spatial layer ever developed to estimate relative ungulate biomass available to Mexican wolves across their range.

| Habitat suitability models
We implemented an additive model with the rescaled variables to produce two sets of habitat models, with and without the UBI map, to assess the sensitivity to this variable. For models without UBI, We present results for the intermediate scenario (optimistic and pessimistic scenarios can be found in Appendix S1, Figures S7 and   S8). Finally, we identified and quantified the areas of high-quality habitat across the modelled distribution of the Mexican wolf in each scenario using the aggregate polygons tool of ArcGIS 10.2. To identify large areas of high-quality habitat, we selected groups of pixels in the upper quartile for each scenario with and without UBI (≥3.57 and ≥2.513, respectively) that represented areas ≥2,500 km 2 in which high-quality pixels were not separated by more than 10 km.

| Distribution modelling of the grey wolf in the southwestern United States and Mexico
According to validation metrics, the niche modelling algorithms that performed best were Maxent, RF, CART, and GAM (Table 1)

| Environmental and anthropogenic habitat variables
The climatic suitability map showed that the highest suitability scores were found in the western portion of the modelled distribution TA B L E 1 Model performance metrics for binary maps generated by ecological niche modelling algorithms  Ungulate density estimations were affected by methodological differences between sources of data. In general, the variance explained with the RF regression models was relatively high for elk but low for the mule deer and white-tailed deer (Table 2)

| Habitat suitability without the Ungulate Biomass Index (UBI)
Results of the additive habitat suitability models excluding the UBI for   Figure 3a).

| Habitat suitability including the Ungulate Biomass Index (UBI)
When the UBI layer was included in the model, the general geographic patterns remained unchanged, but the estimated amount of high-quality habitat changed in most areas (Figure 3b) TA B L E 3 Area estimates (km 2 ) of high-quality patches for the intermediate scenario without and with the UBI layer (see Figure 3 for reference of the spatial distribution of the four regions)

| D ISCUSS I ON
The recovery of Mexican wolves has been a tremendous endeavour undertaken by governments, scientists and conservation agencies in the United States and Mexico over the last 40 years (CONANP, 2009; Moctezuma-Orozco et al., 2011;Servín et al., 2008;Siminski, 2016;USFWS, 1982). Today, more than 190 Mexican wolves roam free in two areas designated for their recovery, with at least 163 individuals in the United States (USFWS, 2020) and about 30 in Mexico (CONANP, 2020). However, the long-term persistence of this subspecies will be enhanced with additional recovery areas, highlighting the need for a comprehensive analysis of habitat suitability for the Mexican wolf throughout its historical range to identify potential areas to conduct specific investigations of local conditions, both ecological and sociological, to select suitable release sites.
We first determined the geographic extent of ecological conditions that potentially could have supported the Mexican wolf.
Different niche-based distribution models concordantly identified a region that included central Arizona and New Mexico southwards into Mexico. This area was then the logical focus of the subsequent evaluation of habitat suitability to refine the areas that meet the ecological conditions that could maximize the probability of successful recovery. Therein, we identified large enough areas (>81,000 km 2 ) of highly suitable habitat to support the recovery of the Mexican wolf, both in Mexico and in the United States ( Figure 3, Table 3).

| Habitat suitability for the Mexican wolf
One of the limitations to evaluate habitat suitability for the Mexican wolf has been a lack of information on prey availability (Araiza et al., 2012;Carroll et al., 2003Carroll et al., , 2014. In this analysis, we addressed this limitation for the first time by modelling relative ungulate density across the entire range of the Mexican wolf using empirical data to calculate an Ungulate Biomass Index (UBI, Fuller et al., 2003). Our analyses partially support hypothesis (5), but further data and analyses are needed to clarify this intriguing question.
The NSMOcc and SSMOcc also hold extensive areas of highly suitable habitat interspersed in a mosaic of varying quality habitat ( Figure 3). The smaller Mexican wolf evolved preying on Coues' white-tailed deer (Odocoileus virginianus couesi) and smaller mammals and birds (Brown, 1983), but seems to have benefited from the multispecies ungulate guild in the northern periphery of its historical range. Finally, the SMOr also holds patches of high-quality habitat and may also play a role in the recovery of wolf populations, but the mountain ranges in that region are naturally smaller and more sparsely distributed than in the SMOcc (Figure 3). This analysis shows that there is ample high-quality habitat in Mexico for the Mexican wolf to be recovered in its historical range.

| Comparison with other habitat suitability analyses
Most previous habitat suitability analyses have focused only on select geographic areas in the United States or Mexico, in part because of the lack of comparable source data for the two relevant countries (Araiza et al., 2012;Carroll et al., 2004Carroll et al., , 2006Carroll et al., , 2014. Carroll et al. (2004) used what geospatial data and tools were available at the time to evaluate habitat throughout the Mexican wolf's historical range and extralimital areas to the north, prior to wolves being released in Mexico. However, their analysis was influenced by the application of a differential spatial adjustment to base mortality risk that decreased the suitability of most habitat in Mexico and boosted suitability of extralimital areas, such as the privately owned Vermejo Park Ranch and Grand Canyon National Park. Despite these adjustments, they concluded that Mexico had suitable quality habitat to sustain 2,600 Mexican wolves (Carroll et al., 2004). On the other hand, the work of Araiza et al. (2012) was not intended to be a comprehensive analysis of Mexican wolf historical range, but rather an evaluation of six areas identified by knowledgeable biologists as holding potential for wolf recovery in Mexico. Unfortunately, this analysis has been portrayed inaccurately (Carroll et al., 2014;Hendricks et al., 2017) as evidence that Mexico does not have enough suitable habitat for wolf recovery. However, their results differed substantially from ours and warrant further discussion.
We found problems in several steps taken by Hendricks et al. (2016) to produce their niche models that we detail in Appendix S1, Extended discussion. First, the authors built two different niche models with Maxent , one "typological," using occurrence records from individuals morphologically recognized as Mexican wolves, and another "genealogical," in which they included additional occurrences from individuals belonging to a genetic cluster referred to as the "southern clade" (Hendricks et al., 2016(Hendricks et al., , 2019Leonard et al., 2005), but morphologically recognized as a different subspecies. The records used for the typological model were highly clustered and did not include available records associated with museum specimens in the southern portions of the historical distribution of the Mexican wolf. Because the authors did not reduce the sampling bias, the resulting model included areas that the Mexican wolves rarely occupied, such as the lowlands of the Chihuahuan Desert (Brown, 1983;Gish, 1977;McBride, 1980). Furthermore, the genealogical model was more problematic. Besides including seven additional northern wolf samples-including one from Nebraska and two from Utah-Hendricks et al. (2016) changed the parameterization of Maxent without explanation to produce a more widespread distribution (Merow et al., 2013), resulting in a misleading model extending far to the north of the known historical and extended distribution, even as far north as Oregon (Goldman, 1944;Hall, 1981;Heffelfinger et al., 2017a;Nowak, 1995;Parsons, 1996).
Nonetheless, even with these changes, the resulting model still did not capture the locations of four of the seven "southern clade" northern wolf records. In our analysis, we used similar sources to obtain the occurrence records, but we avoided these pitfalls and obtained a distribution model that better resembles the historical and extended distribution of this subspecies (Figure 1).
Another substantial difference between our results and those of Hendricks et al. (2016) was that they determined habitat in Mexico, changing abruptly at the international border with the United States, to be mostly unsuitable (Heffelfinger et al., 2017a).

| Informing recovery actions
Our geographic description of ecological conditions associated with the Mexican wolf's historical occurrences coincides with all historical range delineations, as summarized by Heffelfinger et al. (2017a), but made before the more recent advocacy of an expanded range to facilitate recovery only in the United States (Carroll et al., 2014;Hendricks et al., 2016Hendricks et al., , 2017Hendricks et al., , 2019Leonard et al., 2005). The latter position seems especially incongruous because the greater part of the historical range of the Mexican wolf was in Mexico (Heffelfinger et al., 2017a;Parsons, 1996). Some have suggested that a significant portion of the subspecies' range must be occupied before delisting (Vucetich et al., 2006); if such a position were adopted, the Mexican wolf could not be recovered without extensive efforts in Mexico. Furthermore, the recovery of a subspecies outside of the ecological conditions in which it evolved is fraught with legal and ecological concerns (Odell et al., 2018).
The Mexican wolf was extirpated both in Mexico and in the United States mainly because of conflicts with livestock operations and not because of habitat destruction. Nonetheless, in past analyses Mexico has been marginalized or parameterized in such a way as to make it appear unsuitable (Carroll et al., 2014;Hendricks et al., 2016Hendricks et al., , 2019. In contrast, we found that abundant high-quality habitat still exists across the historical range of the Mexican wolf to ensure recovery and allow Mexico to play a vital role in this endeavour. Until now, all efforts to recover the Mexican wolf in Mexico have focused on the northern portion of the Sierra Madre Occidental (SMOcc). Here, we demonstrate that the alternatives are much broader as large tracts of highly suitable habitat remain in the southern SMOcc and to a lesser extent in the Sierra Madre Oriental. In fact, there is a continuum of suitable habitat among the four regions that we identified in our analyses (Table 3) Recovery efforts for the Mexican wolf must favour a collaborative process of common conservation interests to strengthen the currently successful binational effort to return Mexican wolves to the landscape where they can play their historical ecological role.
Therefore, we suggest building on the present analysis to identify additional candidate areas to implement detailed local ecological and social analyses, including prey availability assessments, potential conflicts with livestock, social tolerance, safety issues for the field teams and the ability to maintain a long-term monitoring programme, all with the inclusion of stakeholders in the process.