Introduction

The Nelson’s spiny pocket mouse, Heteromys nelsoni Merriam, 1902, is a monotypic endemic species with an extent of occurrence of around 4 650 km² (^{Cuarón & Vázquez, 2018}). This species has been reported at elevations from 2 480 m to 2 850 m in discontinuous patches of cloud forests in Pinabeto (type locality) and Cerro Mozotal in Chiapas, Mexico and in Volcán Tajumulco in Guatemala (^{Rogers & Rogers, 1992}; ^{Patton, 2005}; ^{Reid, 2009}; ^{Cuarón & Vázquez, 2018}). Only two populations of H. nelsoni are known to occur in Chiapas both in mountain peaks, geographically isolated by a large valley of 25 km with a difference in altitude of 1 600 m, and known by a few specimens deposited in biological collections (^{Ríos, Lorenzo, & Álvarez-Castañeda, 2016}): 1) in the vicinity of Cerro Mozotal (municipalities El Porvenir and Siltepec), and 2) in the location of Pinabeto (municipality Motozintla).

Heteromys nelsoni is listed as endangered by the International Union for Conservation of Nature (IUCN, 2018) and is under special protection by the Mexican Government (SEMARNAT, 2010). Main threats for H. nelsoni are: 1) habitat loss and degradation in relation to human activities, which has evidenced some grade of demographic contraction in part of its distribution (e.g., Cerro Mozotal, Chiapas; ^{Ríos et al., 2016}); 2) stochastic events such as storms and flooding (^{Cuarón & Vázquez, 2018}); and 3) global warming, which can cause contraction in its restricted distribution range (^{Rogers & Rogers, 1992}).

Understanding how species respond to climate change is critical for implementing conservation strategies (^{Thomas et al., 2004}; ^{Rull & Vegas-Vilarrúbia, 2006}; ^{Cahill et al., 2013}). Climate change usually favor the presence of species that can expand their distribution range (^{Levinsky, Skov, Svenning, & Rahbek, 2007}; ^{Moritz et al., 2008}), but endemic and some narrowly distributed species are less likely to be in equilibrium with climate (^{Levinsky et al., 2007}) and exhibit some characteristics that are expected to increase a species’ sensitivity to climate change (^{Walther et al., 2002}; ^{Thomas et al., 2004}). Some examples include low fecundity, poor dispersal capability, high energetic requirements, narrow thermoneutral zone, and obligate relationship to a naturally patchy habitat (^{Smith & Weston, 1990}). ^{Moritz et al. (2008}) report the contraction of the distribution range of four species of mountain mammals and the collapse of two other species. (^{Mcdonald & Brown, 1992}) also report the disappearance of three species of mountain mammals for the Great Basin in Western North America.

In Mexico, climate change models show that there is a tendency of a 10-28 % precipitation decrease (2090) and a temperature increase (2050) of 1.5-2 °C (Banco Mundial, 2013), and while this may decrease the distribution range of some narrow distributed species, testing this hypothesis is very difficult since there are several factors that act at different spatial scales (^{Araújo, Thuiller, & Pearson, 2006}). One of the predicted biological responses to global warming is the upslope displacement of the species. There are often local endemic species in tropical mountain ranges that are distributed near the summits, which may be especially vulnerable to the total loss of habitat due to global warming (^{Raxworthy-Christopher et al., 2008}). However, the vulnerability of most tropical mountain species to extinction has not been well documented (^{Rull & Vegas-Vilarrúbia, 2006}).

In this paper, we attempt to explain the possible impacts of future climate change projections on the distribution of two small mammal species by modeling its potential distributions and altitudinal range variation. We expect the potential distribution of H. nelsoni to decrease due to its restricted distribution to highlands, while that of H. goldmani increases due to its wide distribution. These results can be used to re-assess the conservation status of the species and to propose conservation policies, considering specific areas of H. nelsoni occurrence, since having this kind of information it is fundamental to be able to establish conservation strategies to reduce the loss of biodiversity and protect important functions in the ecosystem (^{Mouillot et al., 2013}).

Materials and methods

Historical records of H. nelsoni and H. goldmani were obtained from the Global Biodiversity Information Facility (GBIF, 2018), as well as the data base of the Colección Mastozoológica of El Colegio de la Frontera Sur at San Cristóbal de Las Casas, Chiapas, México. Additionally, field monitoring was carried out in August and November, 2017 and February-March, 2018 in sampling sites of Pinabeto, Chiapas, 22.6 km away (in a straight line) from Cerro Mozotal, Chiapas (Fig. 1).

Fig. 1 Sampling sites of Heteromys nelsoni in southern Chiapas (stars). The upper right image shows the region of sampling sites in southeastern México and western Guatemala. The lower right image shows the altitudinal differences in the sampling sites. 

A total of 43 specimens were captured using Sherman traps. In Mozotal, 22 sampling days were carried out for a total of 2 720 traps-night, while in Pinabeto, 13 sampling days were carried out for a total of 1 500 traps-night. The external measures were registered for each individual (total lengths of body, tail, hind leg and ear), weight, sex, reproductive condition, habitat type and geographical position with a manual GPS (Garmin eTrex Vista, Kansas, USA). The two species can be distinguished by characteristics such as: the fur of H. nelsoni does not have hard bristles, H. goldmani has a lateral line that does not have H. nelsoni. At the skull, H. nelsoni has the upper third molar is equal to or wider than the premolar, while in H. goldmani the upper third molar is narrower than the premolar (^{Álvarez-Castañeda, Álvarez, & González-Ruiz, 2017}). All capture and handling methods followed the animal care and use guidelines of the American Society of Mammalogists (^{Sikes et al., 2016}).

Potential distribution models were generated using Maxent v.3.4.1 (^{Phillips & Dudík, 2008}). A Pearson correlation analysis was made between all layers, and all those with a significant correlation higher than 0.7 were discarded (^{Gormley et al., 2011}). All duplicate records, as well as those less than 1 km away from each other were discarded (^{Velasco et al., 2016}). We used remaining H. nelsoni (33 records) and H. goldmani (136 records), including altitudes for each one and the bioclimatic layers (version 1.4) generated (^{Hijmans, Cameron, Parra, Jones, & Jarvis, 2005}). A first model with ten replicates was executed to identify the variables with a contribution greater than 0.5 % (^{Warren, Wright, Seifert, & Shaffer, 2014}) and with greater importance for the model according to the Jackknife test of variable importance (which tests the highest gain environmental variable when used in isolation). Layers which did not meet these requirements were discarded. Finally, ten replicate models were run with a random seed to produce distinct model-to-model variation using the bootstrap type (^{Peterson et al., 2015}) and a maximum of 500 iterations, holding out 25 % of the samples for testing in each run of the model. Models were evaluated through the Area under the curve (AUC) and the ROC graphics (^{Hanley & McNeil, 1982}). To reduce uncertainty, the best five models were chosen (highest AUC and lowest omission rate, Table 1), each one was transformed to a binary map considering the logistic threshold of the 10th percentile of the training data (10 percentile training presence, ^{Proosdij, Sosef, Wieringa, & Raes, 2016}). Subsequently, we made a multiplication of these models using map algebra to obtain a consensus map of the current potential distribution for each species (consensus approach, ^{Aguirre-Gutiérrez et al., 2013}). Additionally, following the same procedure, we projected the potential distribution of H. nelsoni and H. goldmani for the years 2050 and 2070. The bioclimatic layers of (Hijmans et al., 2005) were used for the HadGEM2-ES (one of the major models used in the IPCC Third Assessment Report) and CCSM4 (the Community Climate System Model) models. These are two of the main models used in similar analyzes. We used the calculations of the RCP85 scenario (Representative Concentration Pathways) for both models. This scenario estimates that greenhouse gas emissions continue to rise throughout the 21^st century, with an estimated increase of 3.7 ºC for the year 2070 (IPCC, 2013).

Additionally, the Extent of Occurrence (EOO; area contained within the shortest continuous imaginary boundary which can be drawn to encompass all the known, inferred or projected sites of present occurrence of a taxon, excluding cases of vagrancy) and Area of Occupancy (AOO; area within its extent of occurrence which is occupied by a taxon, excluding cases of vagrancy) was determined according to the criteria of the IUCN through Criterion B1 (geographic range in the form of extent of occurrence estimated to be less than 5 000 km²) and B2 (geographic range in the form of area of occupancy estimated to be less than 500 km², IUCN, 2012; IUCN Standards and Petitions Subcommittee, 2017). This, in order to generate information about the risk status of H. nelsoni, maintains consistency with the Red List assessments. For EOO, we created a convex shell for all recent records of both species. For AOO, we created a 2x2 km grid for Chiapas (Mexico) and Guatemala. Then, records of both species were placed, and the number of grids of 4 km² within which there were records of the species was calculated. All analysis was developed using QGis software v.3.2 (QGIS Development Team, 2018).

Results

Bioclimatic variables such as Annual Mean Temperature (ºC) and Precipitation of Driest Month (mm) significantly determined the distribution of H. nelsoni. The same happened when modeling the distribution of H. goldmani, finding that bioclimatic variables such as Temperature Seasonality and Mean Diurnal Range can explain the distribution of this species (Table 2).

As expected, the present maximum suitability area of H. goldmani is greater than that of H. nelsoni, and there is no overlap in both species distributions. In projections to future scenarios, for H. nelsoni the size of its area of maximum suitability shows a slight decrease with CCM4 model (17.94 %) for 2070, while increase with HadGEM2-ES model (23.08 %) for 2070. For H. goldmani, the size of its area of maximum suitability shows a slight increase with CCM4 model (4.94 %) for 2070, while decrease with HadGEM2-ES model (28.29 %) for 2070 (Table 3; Fig. 2).

Fig. 2 Maximum suitability area for Heteromys nelsoni (black) and H. goldmani (light gray) for the two models tested. 

The EOO analysis show an area of 8 077.32 km² for H. goldmani and 487.47 km² for H. nelsoni. According to AOO analysis, H. nelsoni has a 32 km² AOO, and H. goldmani a 160 km² AOO (Fig. 3).

Fig. 3 Records of H. goldmani (light gray) and H. nelsoni (black), in study area. Extent of Occurrence (EOO) is shown for both species (A). AOO for H. goldmani (B) and H. nelsoni (C) is shown, whit 2x2 km grid for all records. 

Discussion

^{Cuarón and Vazquez (2018}) reported an extent of occurrence for H. nelsoni of about 4 650 km². In contrast, our field records and the potential distribution analysis show a distribution surface no greater than 33.44 km² according to one of the used models. Until now there is no evidence of the presence of H. nelsoni in the protected area of the Tacaná Volcano Biosphere Reserve in Chiapas, Mexico, which coincides with the result of its potential distribution, which does not predict suitable habitat sites in the Tacaná Volcano area. An ecological study of mammals in the Tacaná region at altitudes of less than 2 000 m does not mention the presence of this species (^{Mendoza-Sáenz & Horváth, 2013}). This result implies that the habitat available for the species is much smaller than what has been described.

Heteromys nelsoni exists between 2 480 and 2 850 m; ^{Lorenzo, Bolaños-Citalán, & Retana-Guiascón, 2019}), and H. goldmani between 1 250 m and 2 680 m (^{Ríos, Lorenzo, & Álvarez-Castañeda, 2016}). This difference in altitude was not important for the potential distribution models (its contribution value was low and highly correlated with other bioclimatic variables). Therefore, it seems like altitude is not a variable that determines a shift in these species distribution. In contrast, the bioclimatic variables, such as temperature and precipitation are important for the potential distribution of the Heteromys species under study. Our results indicate that future changes in climatic conditions will reduce the area of suitable habitat for H. nelsoni and will favor the presence of H. goldmani. This trend can be observed with the two models tested.

Heteromys goldmani, previously considered a subspecies of the widely distributed H. desmarestianus, but genetically differentiated and supported as a valid species (^{Rogers & González, 2010}; ^{Álvarez-Castañeda et al., 2012}), has been registered in coffee plantations, pine-oak forests and patches of cloud forest in Cerro Mozotal and also patches of cloud forest in Pinabeto (^{Ríos et al., 2016}). It is a habitat generalist species, and it is expected to benefit from future climate change (^{Levinsky et al., 2007}; ^{Moritz et al., 2008}), which is supported by our models that indicate an increase in its distributional range (its distribution area could be extended up to 4.74 % in by 2070, according to the CCSM4 model). However, in case the climate conditions predicted by the HadGEM2-ES model are registered in 2070, this species could also be affected by climate change, despite being a generalist species.

On other hand, H. nelsoni only occurs in the mountain peaks of the Sierra Madre de Chiapas, which present unique bioclimatic conditions that do not allow it to occur in other areas. This mouse is restricted to cloud forests (^{Rogers & Rogers, 1992}), which is an endangered habitat due to unrestricted logging, extensive livestock farming and fire, and presents a low genetic diversity with a very recent demographic contraction in the population of Cerro Mozotal (^{Ríos et al., 2016}). This specialist habitat species is expected to be negatively affected from future climate change (^{Walther et al., 2002}; ^{Thomas et al., 2004}; ^{Levinsky et al., 2007}), which is supported by our models. Although models show an increase of available area for H. nelsoni in 2050, these models also show a similar trend with a reduction of the available area by 2070. By geographically projecting the areas of maximum suitability for H. nelsoni, the increase in temperature seems to create new available sites in a higher altitudinal range. Probably, the availability of habitat for this species is favored by promoting the presence of certain species that form its habitat (^{Liu, Kumar, Katul, & Porporato, 2018}; ^{Mitchell et al., 2018}). However, the prediction for 2070 is a reduction in the size of areas of maximum suitability, probably due to limitations in the altitudinal range, limiting the areas where the species could take refuge. This demonstrates the need to make longer-term projections for these species.

For both species, temperature was the most important variable, measured through the effect of bioclimatic variables. For H. nelsoni the Annual Mean Temperature (ºC) contributed 59.1 % to explain its distribution, and for H. goldmani the Temperature Seasonality contributed with 44.1 %. This is something reported in previous studies, where the temperature and its derivatives in the bioclimatic variables had a significant effect on the size of the future potential distribution of the species. This seems to be explained mainly by the direct effect of temperature on the distribution of habitat, refuge or food for the species of interest (^{Fernández & Hamilton, 2015}; ^{Bradie & Leung, 2016}; ^{Sántiz, Lorenzo, Carrillo-Reyes, Navarrete, & Islebe, 2016}; ^{King, Karoly, & Henley, 2017}).

It is evident that the EOO and AOO values are closer to the IUCN criterion for critically endangered species than to the criteria for endangered species. Considering the bias of this criterion, by not considering non-habitable areas for species of punctual distribution, as is the case with H. nelsoni (IUCN Standards and Petitions Subcommittee, 2017), we suggest this species status be reevaluated mainly in its B criterion (geographic range in the extent of occurrence and/or area of occupancy), and to take specific actions for the conservation of the species and its habitat. We also suggest that the creation of a terrestrial protected area for this species is urgent.

Ethical statement: authors declare that they all agree with this publication and made significant contributions; that there is no conflict of interest of any kind; and that we followed all pertinent ethical and legal procedures and requirements. All financial sources are fully and clearly stated in the acknowledgements section. A signed document has been filed in the journal archives.