nueva página del texto (beta)
Inglés (pdf)
Artículo en XML
Referencias del artículo
Traducción automática
Enviar artículo por email
Citado por SciELO 
Similares en
SciELO 
Permalink
Introduction
Romerolagus diazi (Ferrari-Pérez, 1893) is an endemic species of the central-eastern Trans-Mexican Volcanic Belt (TVB) biogeographic province (Barrera 1966; Velázquez et al. 1996), and it has been listed as endangered species with a decreasing estimated population of 7,000 individuals (Velázquez and Guerrero 2019). This rabbit is also known as the zacatuche, teporingo, or volcano rabbit; it is a monospecific and ancient genus with taxonomic, anatomic, and biogeographic features similar to Pentalagus and Pronolagus (Hoffman et al. 1994). The volcano rabbit is a gregarious species that form groups of two to five individuals, although the age and sex composition of these groups are unknown. To date, knowledge about home range or dispersal capacity is low (Rizo-Aguilar et al. 2014). However, Galindo-Leal and Velázquez (1996) suggested low dispersal capacity than other lagomorphs, and Cervantes and Martínez-Vázquez (1996) proposed a home range of 2,500 m2 based on their field research. In terms of reproduction, the gestation period is longer, and the litter size is smaller in these ancient rabbits compared to other rabbits (Cervantes 1982).
Romerolagus diazi faces intense human pressure due to agriculture expansion, poaching, and development for tourism and because its restricted distribution is surrounded by large cities, including Puebla, Toluca, and Mexico City. These situations modify, fragment, or destroy the specific habitat where R. diazi lives. All of these factors make this lagomorph species a priority conservation target. Also, its geographic distribution rarity increases its vulnerability (Lawler et al. 2003). If the species is extirpated from its known geographic distribution, there is no other place in the world to find it.
The Iztaccíhuatl-Popocatépetl National Park (IPNP) is located in the Mexican states of Puebla, México and Morelos (19.2362° N, -98.6634° W), and it has an area of 39,819 ha and a human population of 244 people (Instituto Nacional de Estadística y Geografía (INEGI 2010); Comisión Nacional de Áreas Naturales Protegidas (CONANP 2013); Secretaría de Medio Ambiente y Recursos Naturales (SEMARNAT 2013). A large part of the geographic distribution of R. diazi (Martínez-Meyer 2005; Farías et al. 2015) is located within the natural protected area of IPNP (Figure 1), which includes pine forest, alpine meadows, rocks without vegetation and the Popocatépetl volcano. In the IPNP, Osuna et al. (2020) identified two of the four linages of the teporingo, one at northern and the other at southern part.
The geographic distribution of a taxon depends on its ecological niche and dispersal ability (Soberón and Peterson 2005). Correlative methods link presence records of a taxon and its associated environmental variables, and therefore, they can be used to produce maps of potential distribution based on the theory of ecological niche (Soberón et al. 2017). In particular, species distribution models are representations in geographical space of the suitability of a place for a species’ presence, where suitability is the mathematical or statistical relationship between the actual distribution and a set of predictor variables (Mateo et al. 2011).
We aimed to prioritize zones with high suitability for the volcano rabbit within Iztaccíhuatl-Popocatépetl National Park (IPNP), using data derived from satellite and fieldwork, relating the utility of correlative methods with helpful conservation strategies.
Materials and Methods
Fieldwork. The fieldwork aimed to obtain georeferenced presence records of R. diazi. Our study area for fieldwork was located on the southern part of the IPNP between the Popocatépetl and Iztaccíhuatl volcanos, a 74 km2 site. Trap cameras were placed in the study area with a separation of 1 km between them. The cameras were installed 50 cm from the ground and were not impeded by vegetation. The cameras were operational from April 2018 through October 2018, with 2,138 camera days. Image processing was performed using Wild.ID 0.9.28 software (TEAM Network 2017). The following information was captured in the databases: project name, camera ID, geographic coordinates, date and time, type of photo, file name, taxonomic identification, number of animals, the person that identified, start date, end date, person who placed and removed the camera, camera model, and institution responsible. In addition to the camera images of volcano rabbits, their presence was recorded by direct and indirect observations. Throughout the study, transects were made to identify volcano rabbits and their latrines visually. The identification of latrines was based on a latrine surface area of approximately 20 cm x 20 cm, at least 20% fresh scat, a uniform discoid shape of feces, and a uniform 5 to 9 mm size of feces (Cervantes and Martínez-Vázquez 1996), with at least ten feces per latrine.
Data analysis. We performed different suitability models at a scale of 30 m to identify places with ideal conditions. To do this, we compiled the data obtained from our fieldwork into one database (model a) and in a separate database (model b) compiled data from fieldwork, GBIF data (biological collections or authors in Appendix 1), and data from the mammal database of the JM055 project funded by CONABIO (Escalante 2014). In both data sets, the locations were reviewed, and the points were filtered with a spatial resolution of 30 m. We used satellite images from November 7, 2017, January 10, 2018, and January 29, 2019, from Landsat 8 OLI/TIRS C1. United States Geological Survey (USGS 2019). First, we performed a radiometric calibration (conversion to reflectance with angular correction) of the Near Infrared (NIR), Red (R), and Shortwave Infrared 1 (SWIR1) bands (Ariza 2013; USGS 2013). Then, we used these calibrated bands to calculate the Normalized Difference Vegetation Index (NDVI; Rouse et al. 1973). To relate the NDVI values to specific vegetation types, six 50-m transects were performed in different areas and elevations of the IPNP. A 50-m rope was placed on the ground, and the vegetation under the rope was collected at one-m intervals. Plant species were identified ex-situ by César Miguel-Talonia (unpublished data). The NDVI values used for this vegetation characterization were from the satellite images from January 10, 2018, with the EPSG projection: 4326 - WGS84. The humidity was obtained using the Normalized Difference Moisture Index (NDMI; USGS 2013).
Figure 1 Iztaccíhuatl-Popocatépetl National Park is located at the center east of the Transmexican Volcanic Belt. The potential distribution for R. diazi (Farías et al. 2015) is shown in the lower-left corner into the IPNP.
We estimated the surface temperature of the area using algorithms, as proposed by Wang et al. (2015) and Avdan and Jovanovska (2016). To do this, we first calculated the TOA (Top of Atmospheric Spectral Radiance; Barsi et al. 2014) from band 10 (TIR1) due to high uncertainty in the values of band 11 (TIR2; Wang et al. 2015). The value of O10 reported by Wang et al. (2015) is 0.29 (W·m−2·sr−1·μm−1), based on USGS files for dates before February 3, 2014. For later dates, as in our study, this value should not be included because the downloaded data is already processed, including this value (Wang et al. 2015). Secondly, we converted the reflectance to the brightness temperature (BT; USGS 2013), derived from an approximation of the Planck radiance function using the constants that appear in the product metadata. Third, we calculated the proportion of vegetation (Pv) from the range of NDVI (maximum and minimum) depending on the area (Carlson and Ripley 1997; Sobrino et al. 2004; Dash et al. 2005). Sobrino and Raissouni (2000) proposed values of NDVIS = 0.2 and NDVIV = 0.5 for global conditions. For particular areas, NDVIS and NDVIV can be extracted from the NDVI histogram (Sobrino et al. 2008). Next, we calculated the emissivity of the Earth’s surface (LSE, Land Surface Emissivity). In the wavelength range of Landsat 8 band 10, the emissivity can be calculated using the simplified threshold method for the NDVI (SNDVITHM; Sobrino et al. 2008). To calculate the emissivity, we used the values for general conditions proposed by Sobrino and Raissouni (2000). For NDVI values less than 0, the emissivity was assigned a value of 0.991; for NDVI values between 0 and 0.2, we considered that the surface is covered by soil and corresponds to an emissivity of 0.966; values between 0.2 and 0.5 corresponded with a mixture of soil and vegetation, while values above 0.5 were considered fully covered by vegetation, with an emissivity of 0.973 (Wang et al. 2015). We obtained the land surface temperature using the method of Stathopoulou and Cartalis (2007).
Finally, we used topography because it could affect species distributions. Therefore, we calculated the slope using the Continuous Mexican Elevation 3.0 (INEGI 2013). All layers obtained were projected to SCR: EPSG 4326 - WGS84.
For model calibration, and because our study is focused on the Iztaccíhuatl-Popocatépetl National Park, whose boundaries are based on natural factors (CONANP 2013), we considered the polygon of the park with a 2 km buffer as the accessible area. The environmental suitability was modeled using MaxEnt version 3.4.1 (Phillips et al. 2006). The model’s parametrization strongly influences the final results, and the output format has implications for its interpretation. Therefore, we used the cloglog transformation as the output format, derived from the interpretation of MaxEnt as a non-homogeneous Poisson process (IPP), producing a more robust theoretical justification than the logistic transformation (Phillips et al. 2017). Ecological theory suggests that response curves are unimodal for fundamental niches (Austin 2007), so quadratic features may be more appropriate.
On the other hand, linear features may be sufficient when the species’ niche is cut on one side of the unimodal curve (Merow et al. 2013). For some authors, such as Radosavljevic and Anderson (2014), the regularization multiplier value must be greater than the default value (1) to achieve the model’s optimal complexity. It is also advisable to eliminate highly correlated environmental layers since the features are already strongly correlated (Merow et al. 2013). In this work, to (a) have a unique criterion, (b) facilitate replication and (c) maximize the robustness of the resulting models, we used the ‘kuenm’ package (Cobos et al. 2019) in R Studio version 1.2.1k (RStudio Team 2018). We used the kuenm_cal function to produce candidate models combining all possible features (29 combinations) and the following values for the regularization multiplier: 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5 and 6. To develop candidate models, we used different combinations of predictive variables: the first set with all the predictive layers, the second set only with the layers that did not show collinearity (based on a collinearity matrix and estimation of the VIF; R routine in Appendix 2), and the third group without topographic layers. The data partition for training and testing was done with the ‘ENMeval’ package (Muscarella et al. 2014) in RStudio version 1.2.1k (RStudio Team 2018) using the random k-fold method where 75 % of the data were selected for training and the remaining 25 % for testing. We evaluated each candidate model using the kuenm_ceval function based on the statistical significance given by the partial ROC curve, the omission rate, and the complexity of the model, calculated with the Akaike information criterion corrected for small sample sizes (AICc; Akaike 1974; Burnham and Anderson 2002; Warren and Seifert 2011). The models that met the requirements were selected, and we generated the final models with the kuenm_mod function using ten bootstrap replicates with a cloglog output format. Finally, we evaluated those final models with the kuenm_feval function, based on an independent data set excluded before model calibration, and we selected the model with the best performance.
The conversion of the suitability model (with continuous values from 0 to 1) to a binary model (only absence and presence) was performed using the ‘dismo’ package (Hijmans et al. 2017) with an acceptable threshold of E = 10 (Peterson et al. 2008). To establish a priority area for the conservation of R. diazi within the IPNP, the suitability model was converted to a distribution model reduced to an area of greater suitability with a threshold of E = 50, because some work such as Waltari and Guralnick (2008) suggests that a more stringent threshold provide high consensus presence and better refuge area. To obtain one unique distribution model for both thresholds, we summed all of the final binary models for each date and data set, and we selected the pixels where at least half of the models predicted presence.
Polygons for conservation. Finally, to obtain priority conservation polygons, we used the distribution model reduced to an area of greater suitability (threshold E = 50) and a road map (Google Inc. 2019) from Google Maps Layer TMS (Tile Map Service) using the XYZ Tiles plug-in. We selected this distribution model’s continuous areas that did not reach crossroads in QGIS 3.4.4 (QGIS Development Team 2019). Subsequently, we merged the chosen areas and produced four minimum convex polygons using the QGIS convex envelope function. We calculated the area and perimeter for each polygon, and to characterize them, their external vertices were extracted with their respective latitude-longitude coordinates. Furthermore, we obtained the elevational range for each polygon using the Continuous Mexican Elevation 3.0 (INEGI 2013) with a 90 m-resolution.
Results
Fieldwork. In the camera traps, we obtained photos of 33 individuals of R. diazi. Adding the visual identification of latrines and direct observations of individuals, we obtained 82 occurrence records for the species. There were nine main plants recorded along the vegetation transects: Pinus hartwegii, Lupinus montanus, Eryngium proteiflorum, Draba jorullensis, Senecio sp., Festuca tolucensis, Muhlenbergia sp., Trisetum sp. and Calamagrostis tolucensis. We compared the NDVI value with the georeferenced transects and obtained a vegetation characterization, which we mapped. The results were the following: NDVI < 0, without vegetation; NDVI between 0 and 0.1, alpine meadow with C. tolucensis and bare soil and rock; NDVI between 0.1 and 0.2, meadow with Muhlenbergia sp. / C. tolucensis; NDVI between 0.2 and 0.3, meadow with Muhlenbergia sp.; NDVI between 0.3 and 0.4, meadow-forest ecotone; and NDVI > 0.4, pine forest.
Suitability models. For the model using only data from fieldwork (model a), we had 51 points for calibration (13 points used as testing data and 38 points as training data) and 17 points for independent evaluation. For the model with mixed data (model b; fieldwork, GBIF, and JM055 project), we had 66 points for calibration (50 points as training data and 16 points as testing data) and 22 independent points for final evaluation. Overall, 4,466 candidate models were evaluated (Table 1).
Table 1 Results of the evaluation of the final models (model a with fieldwork data only and model b with mixture data) for each date. P-value-pROC is the value of ‘p’ given in the Partial ROC analysis. OR is the omission rate. Features types: T = threshold, L = linear, Q = quadratic, P = product and H = hinge. AICc is the Akaike Information Criterion with a correction for small sample sizes. RM is the regularization multiplier.
| Models date | pvalue-pROC | OR | AICc | Feature types | Variables | RM |
|---|---|---|---|---|---|---|
| Model ‘a’ dates | ||||||
| Nov 7, 2017 | 0 | 0.08 | 1154.32 | T | all | 2 |
| Jan 10, 2018 | 0 | 0.23 | 1147.74 | L, Q, P, T | all | 1.5 |
| Jan 29, 2019 | 0 | 0.08 | 1161.14 | T | no collinearity | 1 |
| Model ‘b’ dates | ||||||
| Nov 7, 2017 | 0 | 0.06 | 1608.13 | Q, T | all | 4 |
| Jan 10, 2018 | 0 | 0.13 | 1551.30 | P, T | all | 2.5 |
| Jan 29, 2019 | 0 | 0.06 | 1636.80 | L, Q, T, P, H | all | 5 |
Potential distribution model. The area potentially occupied by the volcano rabbit, based on this final distribution model, reached 132.5 km2 within the IPNP (Figure 2). Compared to the potential distribution (Martínez-Meyer 2005; Farías et al. 2015), the volcano rabbit’s distribution was reduced to 33 % of the IPNP polygon.
Potential distribution model reduced to a greater suitability area. The more suitable area covered 7 km2 (Figure 2). We overlapped this area with the vegetation characterization. We found that bare soil with sparse vegetation occupied 0.29 km2 of this most suitable area, alpine meadow dominated by Muhlenbergia sp. and C. tolucensis covered 1 km2. The meadow dominated primarily by Muhlenbergia sp. occupied 4.6 km2 of greatest suitability, being the most suitable vegetation for the volcano rabbit.
Proposal of polygons for conservation. We obtained four polygons as refuges for the conservation of R. diazi inside the IPNP (A, B, C & D; Figure 2). Polygon A had an area of 0.69 km2 and a perimeter of 3.33 km, with an average altitude of 3,912 m, a minimum altitude of 3,848 m, and a maximum altitude of 4,029 m. Polygon B represented an area of 1.13 km2 with a perimeter length of 4.39 km, and the mean altitude is 3,931 m with a minimum and maximum altitude of 3,835 m and 4,054 m, respectively. For Polygon C, the selected area had 1.61 km2 and a perimeter of 4.92 km, where the mean altitude is 3,874 m. The minimum altitude is 3,780 m, and the highest value is 3,987 m. Finally, Polygon D had an area of 3.14 km2 and a perimeter of 9.1 km, and its average altitude is 3,946 m with a minimum value of 3,823 m and a maximum elevation of 4,066 m. The coordinates for each external vertex are shown in Appendix 3. The total proposed area to prioritize for the conservation of R. diazi within the IPNP occupies 6.57 km2.
Figure 2 The potential distribution model for R. diazi in the Iztaccíhuatl-Popocatépetl National Park and potential distribution model with the area reduced to those with greater suitability. All the models obtained with the threshold E = 10 or E = 50 were added together, and the pixels where at least 50 % of the models predict presence were selected. On the right side, four polygons (A, B, C, and D) are shown based on the potential distribution model E = 50 and Google Roads map.
Discussion
Ecological niche modelling and its use as a hypothesis of potential or actual geographical distribution are helpful in prioritizing conservation areas (Sánchez-Cordero et al. 2004). The IPNP is a decreed conservation area; however, administration and management’s decision-making should coincide under a theoretical-practical framework.
The information derived from satellite products is a powerful tool for generating models with better performance on detailed scales (e. g.,Rödder et al. 2016; Vila-Viçosa et al. 2020). In this study, we included the NDVI as an abiotic predictor, although it can also be interpreted as a biotic predictor (He et al. 2015). Different plant species have different leaf anatomy that results in variations in the reflectance captured by remote sensing (He et al. 2009). Although it is possible to find areas with diverse vegetation where the NDVI calculated could be similar, the comparison between studies situated in other geographical regions should be taken cautiously. Identifying georeferenced vegetation and its subsequent relationship with the NDVI values allows the replication of our research. For the future, the addition of NDVI values to flora catalogs, discriminating seasonally, altitudinally, and latitudinally, will facilitate sampling and knowledge of the vegetation of an area and its subsequent relationship with the geographical distribution of animals.
This study included different sets of variables from various data sources to obtain a broader spectrum of information. Because statistical methods controlled the evaluations, and we cannot be ensured that the results of the selected final models guarantee accurate geographic distribution, we decided to combine all the models and select only those pixels where at least half of them coincided with increasing the certainty of the results. Distribution models should be generated with high precision and strictly interpreted not to incur unhelpful or unproductive conservation practices that may put the target species at risk and other taxa. However, a potential distribution model from a suitability model is strongly influenced by the different existing configurations and settings. Thus, it is advisable to use methods that allow standardizing the calibration process of a model.
The distribution of R. diazi showed that the volcano rabbit’s occupation could reach most of the IPNP surface (Martínez-Meyer 2005; Farías et al. 2015). However, this is uncertain, and the IPNP administration does not have the information required to determine the appropriate area for implementing a conservation strategy. For this reason, the search for a more realistic distribution area is necessary to provide that basic knowledge. At a 30-meter m scale, the potential distribution model significantly reduced volcano rabbit occupancy within the IPNP polygon to up to 33 % park’s total area. This difference may be caused by the use of a different scale, but also by the use of different variables or the parameterization of the model, so comparisons should be made with caution.
On the other hand, the environmental predictors used corresponded only to the dry season in México. Therefore, the resulting models can be considered informative for the dry season. However, since R. diazi has a very small home range (Cervantes and Martínez-Vázquez 1996), the distribution of this rabbit will likely be very similar in the rainy season.
Several studies have been published about the volcano rabbit over the past few decades. Velázquez and Bocco (1994) considered agriculture one of the factors that most threatened this species on the Tláloc and Pelado volcanoes. They established areas with different degrees of suitability based on this risk factor. Rizo-Aguilar et al. (2015) showed that vegetation structure and altitudinal range are directly related to the abundance of R. diazi. Moreover, the percentage of meadow cover, short grass, and the scrub cover have a positive relationship with the abundance of the volcano rabbit.
In contrast, the closed forest, tall grass, cattle pasture, hunting, bare terrain, and slope have a negative relationship with its abundance (Hunter and Cresswell 2015). We demonstrated that the suitability and potential distribution are not uniform in the IPNP through modelling techniques and the cross-linking of information with data derived from satellite products. The alpine meadow dominated by Muhlenbergia sp. is the most suitable area for R. diazi in the IPNP, followed by the alpine meadow composed of Muhlenbergia sp. and Calamagrostis tolucensis, both belonging to the Poaceae family. The relationships between the volcano rabbit and plant communities were studied by Velázquez and Heil (1996) in the Pelado and Tláloc volcanoes, where the associations of Festuca tolucensis and Trisetum spicatum-Festuca tolucensis had the greatest abundance of R. diazi. While the floristic study of our work was not the main objective, the inclusion of the plants allowed us to corroborate the importance of alpine meadow without trees in the distribution of R. diazi. Conservation strategies should focus on preserving the alpine meadows in the prioritized polygons, avoiding tourist and unskilled personnel’s access to those specific areas. Human activities in high-mountains ecosystems adversely affect animals (Pęksa and Ciach 2015). The speed of road traffic and the emission of noise in the surrounding areas should be controlled (Garriga et al. 2012), along with other practices to mitigate human impact, such as avoiding the discharge of inorganic and organic waste. In the coming years, scientific studies within the IPNP should be carried out jointly with the park administration and human communities to provide updated tools for decision makers.
Our conservation polygons for R. diazi are based on the most suitable areas and existing access roads in the IPNP because of the significant influx of tourists. Polygons were selected using a remote sensor from three specific dates, so we expected that the vegetation, humidity, and the land surface temperature would change at other times. Besides, our polygons have an average altitude upon 3,800 masl, so they could be helpful in a climate change scenario according to the ascent of the lowest altitudinal limit of the volcano rabbit’s distribution (Anderson et al. 2009). Furthermore, the relative abundance of this rabbit is significantly higher in altitudes above 3,600 masl with abundant bunchgrass cover (Osuna et al. 2021). According to the significant evolutionary units proposed by Osuna et al. (2020), these polygons are located in the area of the Nevada south unit and can also be useful as a natural refugee for genetic management.
Global warming in the medium-long term could lead to changes in the distribution of the volcano rabbit since in the mountains there is an altitudinal temperature gradient that will likely lead species to move uphill (Gottfried et al. 2012). This phenomenon may eventually result in an inability to reach optimal conditions, and the species may become extinct (Colwell et al. 2008). Therefore, we consider that the next step for conserving the teporingo could be selecting the most suitable areas within the IPNP under alternative climate change scenarios.
In addition to the proposed polygons, the most suitable areas identified by niche models should be used as a theoretical and practical basis to propose and execute any conservation strategy in situ. In particular, the vegetation within the reduced distribution of volcano rabbit should be conserved. The meadow with Muhlenbergia sp. is a priority area, and reforestation with pine in them is not advisable.
The boundaries of the IPNP were well designed to protect R. diazi. Although this park has recovered some of its forest covers and has a low amount of transformed area (Aguilar-Tomasini et al. 2020), the conservation of an area does not end with its declaration as a Natural Protected Area. Continuous updating is needed to provide precise conservation tools within protected areas. The abiotic and biotic conditions of the park polygon will vary. The methodological tools in the future will provide various techniques that will enable the rise of information about one of the most threatened and emblematic lagomorphs of México: the volcano rabbit.
