Introduction
Updating the distribution pattern of all endemic and endangered species proves relevant (Smith et al., 2020), and most critically, those that have been controversial, as has been for the volcano rabbit, Romerolagus diazi. Hoth et al. (1987) conducted the most thorough study in this area 35 years ago. Controversial new evidence has contested original findings (Gonzalez et al., 2014; Monroy-Vilchis et al., 2020). The volcano rabbit, locally known as zacatuche, an endangered species (Velázquez & Guerrero, 2019), is the smallest lagomorph and endemic to the central mountains of the Trans-Mexican Volcanic Belt, specifically in the Sierra Chichinautzin and Sierra Nevada which comprises the Popocatépetl and Iztaccíhuatl volcanoes. Its range covers 386 km2 (Velázquez, 1994), though recent studies suggest the area might be larger (Rizo-Aguilar et al., 2015). It is restricted to bunchgrasses (Muhlenbergia spp., Festuca spp.) within forests at elevations of 2,800-4,200 m (Osuna et al., 2021).
Species with a high level of habitat specificity are often associated with ecological factors as proximal drivers to explain their distribution (Ottaviani et al., 2020). Long-term underlying geo-ecological factors (García & Di Marco, 2020), as well as short-term anthropic factors, have also proven to be relevant to explain distribution patterns of endemic and endangered species (López et al., 1996; Uriostegui-Velarde et al., 2018; Velázquez, 1993). This is even more relevant when there are many sympatric endemic species (Fa et al., 1992), so that habitat, rather than one species on its own, must be considered endemic and endangered (Velázquez & Heil, 1996).
The volcano rabbit and its habitat have experienced human-caused and climate change threats (Anderson et al., 2009; Velázquez et al., 2011). Current research on the species has confirmed that dense bunch grassland habitats favor its presence (Hunter & Cresswell, 2015; Rizo-Aguilar et al., 2015; Uriostegui-Velarde et al., 2018). Monroy-Vilchis et al. (2020) recently found that similar dense bunchgrass land habitats are unsuitable. Hence, ecological conditions seem to be only part of the drivers explaining the volcano rabbit distribution pattern. Local surveys of limited scientific outreach have shown that many other areas have been overlooked. These areas may be potentially suitable habitats (Osuna et al., 2021; Velázquez & Guerrero, 2019). To date, the factors that determine the occurrence or absence of the zacatuche at the local and regional level have been a poorly documented aspect.
This research aimed at comparing 2 areas (Monte Tláloc and Nevado de Toluca) with similar ecological habitat characteristics but with different geological histories, where fieldwork in both areas was extensive and the presence of the volcano rabbit has been controversial. The results are discussed, considering their implications for biogeographical conservation contexts.
Materials and methods
Our research took place in areas in a dispute concerning the presence of the volcano rabbit, namely, Sierra Nevada and Nevado de Toluca. Sierra Nevada comprises the Iztaccíhuatl, Popocatépetl, Telapon, and Tláloc volcanoes (Fig. 1). This area was formed around 1.4 Ma to recent (Arce et al., 2003; Espinasa-Pereña and Martin-Del Pozzo, 2006). Nevado de Toluca is one massive structure, locally known as Xinantecatl, that was formed in the Late Pliocene - Holocene around 2.6 Ma to recent (Arce et al., 2003; Astudillo-Sánchez et al., 2017; Table 1). These 2 areas are detached from the Sierra Chichinautzin, where the volcano rabbit has been systematically reported as abundant (García et al., 2018; Rizo-Aguilar et al., 2015; Velázquez, 1993). This is also the case with the Iztaccíhuatl and Popocatépetl volcanoes, where there is well-documented evidence of the volcano rabbit presence. The Telapon Volcano was also thoroughly surveyed recently with no evidence of the volcano rabbit, although Osuna et al. (2020) reported its presence. Our current research focuses on the last 2 controversial areas, namely the Tláloc Volcano (locally known as Monte Tláloc) and the Nevado de Toluca Volcano.

Figure 1 Area of study: Tláloc and Telapon are in the Sierra Nevada and Nevado de Toluca (both in green color). Present protected areas in Sierra Nevada (Izta-Popo-Zoquiapan National Park) and Nevado de Toluca (both are delineated by dashed lines). The new volcano rabbit distribution area in Tláloc and Telapon is not fully embraced by the Protected Area in Izta-Popo National Park. Map by Luis Antonio García Almaraz.
Table 1 Geological history of the volcanoes Nevado de Toluca and Tláloc. Source: Macías et al. (1997); Montero (2002); Macías (2005); D´Antonio (2008); García-Tovar (2011); García-Palomo (2015); Weber et al. (2019).
| Tláloc | Nevado de Toluca | |
|---|---|---|
| Type of volcano | Stratovolcano | Stratovolcano |
| Age | 2.6 million years | 1.8 million years |
| Eruptive activities | 1. 2.6 - 1.15 million years | 1. 1.82 and 1.58 million years |
| 2. 42,000 and 10,500 years | 2. 14,000 to 12,500 years |
The 2,441 sampling sites were surveyed from April 2017 to November 2020. Out of these places, 634 were from the Tláloc Volcano and 1,807 from the Nevado de Toluca. These sampling sites were located above 3,400 m asl, along contour lines (BOLFOR et al., 2000), with an elevational separation of 100 m (Mayer & Ott, 1991). The sampling sites were circles of 0.1 ha with a 17.86 m radius. Sites were systematically distributed on each curve at a 300 m equidistance. We followed Dauber (1995) to obtain the minimum sampling intensity (0.89% recommended and 1.36% achieved). All sampling sites were located by UTM coordinates and elevation, and data on dominant plant species were recorded. Plant species were identified in situ to genus and species, using taxonomic guides and local knowledge. Each sampling unit was characterized by slope steepness, slope exposure (°), percentage of occupation of the dominant herbaceous and shrub species (calculated in m2), evidence of recent fire (< 1 year), rocky areas (% coverage), if there was any type of road or trail., and other reference data for the sampling site (e.g., reforestation, extraction, ravine, associated fauna). In addition, all trees were inventoried (≥ 7.5 cm of normal diameter), recording their normal diameter and total height. Evidence of the volcano rabbit was recorded through droppings, direct sightings, and carcasses.
Following Velázquez (1994), the abundance was estimated by latrine counts (group of 30 or more pellets) in a 9.78 m (300 m2) radius within the sampling site. Interpolation was carried out with the Natural Neighbor method (Childs, 2004; Etherington, 2020; Sibson, 1981) to calculate the area occupied by the species in ArcGis Desktop software v. 10.8 (ESRI, 2019).
Following Velázquez and Heil (1996), we conducted Canonical Correlation Analysis (CCA) habitat analyses (CANOCO v. 4.5; ter Braak, 2002) to test habitat affinities among the study areas, where the largest part of the variation could be explained by the environmental and floristic variables. In addition, data on the presence or absence of the volcano rabbit at all the sampling sites were subjected to factorial analysis with the extraction method of principal components with varimax rotation (Kaiser, 1974). We ran this in SPSS Statistics v. 26.0 (IBM Corp., 2019), considering the variables elevation (m asl), exposure and slope (°), top of the trees, soil cover percentages (rock, herbaceous, shrub), fire, habitat, records of Sylvilagus sp. and R. diazi (through latrines), reforestation (management practices), and road proximity. To calculate the elevation range with the highest presence of zacatuche, the Kaiser-Meyer-Olkin suitability measure and Bartlett´s test of sphericity were performed (Bartlett, 1950).
To evaluate habitat preferences, each habitat type was categorized considering the dominant species of each soil cover (herbaceous, shrub, trees); the observed frequency of the latrine number in each habitat was recorded. The Pearson’s chi-square goodness-of-fit test was applied to obtain frequencies. The result of this analysis was represented following Monroy-Vilchis and Velázquez (2002) .
Results
For the Tláloc Volcano, 4 types of habitats were surveyed: pine forest-bunchgrass land (65 sampling sites with volcano rabbit latrines), alder forest (2 sampling sites with latrines), cypress forest (no evidence of the volcano rabbit), and other habitats (no evidence of the volcano rabbit). In contrast, in the Nevado de Toluca, 3 habitats were surveyed: pine forest-bunchgrass land, alder forest, and other habitats (Fig. 2). No evidence of the volcano rabbit was found in the 1,807 sampling sites in the Nevado de Toluca.

Figure 2 Ordination diagrams showing habitat affinities among the study areas. The triangle symbols represent plant species, whereas arrows indicate variable locations within the ordination diagram. The top diagram (denoted as A) shows the Tláloc Volcano where the 4 plant communities depicted by their dominant species (here listed) occurred. The bottom diagram (denoted as B) shows the Nevado de Toluca Volcano where 3 out of the 4 plant communities depicted by their dominant species (here listed) occurred. (A) The Tláloc Volcano: 1, pine forest-bunchgrass land: Pinus hartwegii-Senecio cinerarioides-Festuca- Barkleyanthus salicifolius-Lupinus montanus- Agrostis-Calamagrostis. 2, Alder forest: Alnus jorullensis-Roldana platanifolia-Pinus pseudostrobus-Senecio-Salix cana-Acaena elongata-Gnaphalium-Ageratina pazcuarensis-Castilleja pectinata-Trisetum-Abies religiose. 3, Cypress forest: Cupressus lusitanica-Arbutus xalapensis-Cirsium jorullense-Ribes ciliatum-Quercus laurina-Symphoricarpos microphyllus-Rumex acetosella-Ribes microphyllum-Baccharis conferta. 4, Other habitats: Juniperus monticola-Robinsonecio gerberifolius, Cirsium nivale-Roldana angulifolia-Alchemilla procumbens-Senecio toluccanus. Axis eigenvalues l: 1: 0.368, 2: 0.065, 3: 0.044 and 4: 0.035. (B) The Nevado de Toluca Volcano: 1, pine forest-bunchgrass land: Pinus hartwegii-Senecio cinerarioides-Festuca-Barkleyanthus salicifolius-Lupinus montanus-Agrostis-Calamagrostis-Eryngium proteaflorum-Penstemon gentianoides-Senecio tolucanus- Ribes microphyllum-Muhlenbergia. 2, Alder forest: Alnus jorullensis-Roldana platanifolia-Pinus patula-Senecio-Acaena elongata-Castilleja toluccensis- Symphoricarpos microphyllus-Baccharis conferta-Roldana angulifolia-Stipa-Quercus laurina-Trisetum-Abies religiosa. 3, Other habitats: Salix cana-Cupressus lusitanica, Pinus montezumae-Arbutus xalapensis-Buddleja cordata. Axis eigenvalues l: 1: 0.182, 2: 0.060, 3: 0.023 and 4: 0.015.
On the Tláloc, the presence of R. diazi was recorded in 67 of the 634 sites (Fig. 3). Most were on the southwest slope, which covers 1,537 ha of the volcano rabbit habitat in the sampled area. According to the latrine number per surface, the abundance of R. diazi on Tláloc was 0.047 latrines / m2. The elevational distribution ranges between 3,400 and 3,900 m asl, with a higher abundance between 3,700 and 3,800 m asl (p < 0.05, 95% confidence) (Table 2), as well as in sites with evidence of recent burning (25 sampling sites) and reforestation (12 sampling sites).

Figure 3 Abundance and distribution of Romerolagus diazi on Tláloc Volcano (1,537 ha). Colors contrast different vegetation types and areas comprising different volcano rabbit abundances. Low: 0.0026-0.0279 l/m2; medium: 0.0280-0.0532 l/m2; high: 0.0533-0.1921 l/m2. Sampling sites surveyed are denoted by white spots. Map by Luis Antonio García Almaraz.
Table 2 Contrasts between Kaiser-Meyer-Olkin suitability measurement and Bartlett´s sphericity test. Both measurements are consistent with the 3,700-3,800 elevation range as the most suitable one for the presence of the volcano rabbit on Tláloc.
| Kaiser-Meyer-Olkin suitability measurement | 00.626 | |
|---|---|---|
| Bartlett´s sphericity test | Chi-squared | 948.054 |
| Degrees of freedom | 78 | |
| Level of significance | p < 0.01 | |
There was significant variation in the frequency of volcano rabbit latrine among habitats (Fig. 4). The pine forest-bunchgrass land (10% of the total area) and the alder forest (0.16% of the total area) habitats had higher frequency values than expected.

Figure 4 Observed and expected frequencies among habitat types in the Tláloc Volcano. Positive values represent volcano rabbit habitat preference greater than expected, while negative values represent volcano rabbit habitat preference less than expected (CC = 274.87, df = 3, p < 0.05).
The Principal Component Analysis was the relationships between variables and the influence on each component (Fig. 5). According to this, fire and reforestation variables were positively correlated with each other. This means that the presence of any of these variables in the highland pine forest and the forest bunchgrass land habitat increases the probability of finding R. diazi.
Discussion
Our results demonstrate that ecological conditions are not the only driving factor to explain the present distribution pattern of the volcano rabbit. Here it is documented that the Tláloc and Nevado de Toluca volcanoes share similar ecological characteristics. They also share these with those reported in the Sierra Chichinautzin, Iztaccíhuatl, and Popocatépetl volcanoes. These are places where the volcano rabbit’s presence has been proven indisputably (Velázquez & Guerrero, 2019). In Figure 2, we documented the structural and species composition similarities among habitats on Tláloc and Nevado de Toluca. Velázquez and Heil (1996), as well as Hunter and Cresswell (2015), strongly state that ecological factors were key drivers of the presence of the volcano rabbit. The intensive sampling conducted in this study (as shown in Figure 3) leaves no doubt that high-elevation habitats from these 2 volcanoes are alike ecologically.
The presence of the volcano rabbit on Nevado de Toluca was reported by local farmers. The most academically outstanding evidence of this was given by González et al. (2014) in 1998. However, we assume that this evidence was either erroneous or derived from an introductory exercise that was done in at least 2 attempts (pers. com.), therefore, there were never native populations of R. diazi on Nevado de Toluca. No trace of the current presence of the volcano rabbit was found on Nevado de Toluca despite all the ecological affinities. This result supports the contribution of Hoth et al. (1987) and, more recently, of Murga-Cortés et al. (2020) and Monroy-Vilchis et al. (2020), who conducted photo-trapping and reached the same conclusion.
Our findings let us infer that geological and biogeographical attributes play a role in explaining the absence of the volcano rabbit on Nevado de Toluca. The Tláloc Volcano arose 1.8 million years ago (Osuna et al., 2021) and the Nevado de Toluca arose 2.6 million years ago (García-Palomo et al., 2002). These 2 sites have gone through many volcanic events. Nonetheless, the most recent volcanic activity in the area has only been experienced in the Nevado de Toluca and the Popocatépetl (this volcano is still in a period of activity).
Furthermore, recent research using ultraconserved genetic elements among lagomorphs (Cano et al., 2021) demonstrated that the volcano rabbit diversified from its ancestor during the Pliocene/Miocene transition (time scale: 5.33 Ma), while Osuna et al. (2020) estimate that it began its diversification ca. 1.4 Ma (Sierra Nevada and Sierra Chichinautzin). As stated by Montero (2002) and Siebe and Macías (2006), the Sierra Nevada and Nevado de Toluca volcanoes developed during the Pleistocene (time scale: 2.5 ~ 0.1 Ma). During the Late Pleistocene and the Upper Holocene (around 0.01 million years ago), many drastic climatic changes took place. These changes impacted species distribution patterns.
Based upon the present results and those of Cano et al. (2021), we postulate that the populations of R. diazi found refuges in high volcanoes during the ice retreat of the Early Holocene. The volcano rabbit populations were partially depleted on Popocatépetl and totally depleted on Nevado de Toluca because of recurrent eruptions during the transition from the Late Pleistocene to the Upper Holocene (Siebe & Macías, 2006). This is without discarding the urban expansion and overgrazing that occurs in the Nevado de Toluca, as there are human settlements up to 3,500 m asl; human disturbance of habitats advances from the bottom up, reducing and isolating them more. Some of the consequences that can come with rising temperature, as well as changes in precipitation, are the extinction of species and the decline of their populations (Domínguez-Pérez, 2007); areas potentially habitable by the zacatuche tend to be confined to the higher elevation zones.
Romerolagus survives from the late Pleistocene, as its presence was recorded from a tooth belonging to a zacatuche in Valsequillo, Puebla (Cruz-Muñ̃oz et al., 2009), although it remained at the Iztaccíhuatl and Tláloc volcanoes of Sierra Nevada during the Late Holocene. The ecological effects of climate change during the Pleistocene led to the loss or fragmentation of habitats (Koch & Barnosky, 2006), which probably completely extinguished the habitable areas for R. diazi in Valsequillo. Later, during the Northgrippian and Meghalayan Holocene periods, it expanded its present distribution to the Sierra Chichinautzin. This hypothesis is coherent with the theory of island biogeography (MacArthur & Wilson, 1967), which is based on the principle that large, connected islands support greater resilience compared to small, isolated islands. This hypothesis is similar to Luna-Vega (2018), who sustained that Central Mexico has been subject to paleoclimatic, tectonic, and glacier advance and retreat events that have caused contraction, isolation, differentiation, speciation, and range expansion of local species. The Popocatépetl and Iztaccíhuatl volcanoes function as biogeographic islands in the midst of warmer climates and diverse types of vegetation, limiting the migration of the zacatuche. In addition, the Pleistocene-Holocene boundary extinction of megafauna was important in reducing predation or vegetation change associated with the loss of disperser species as it altered the distribution of smaller species such as the zacatuche (Ferrusquía-Villafranca et al., 2010).
The present distribution range of the volcano rabbit includes the Tláloc Volcano in the Sierra Nevada and excludes Nevado de Toluca. Although Tláloc is adjacent to Iztaccíhuatl, one of the larger and potentially better-protected areas of habitat (Hunter & Cresswell, 2015), 35 years ago, periodic visits were made in this area without finding evidence of the volcano rabbit (Hoth et al., 1987). Based on the above, it is possible to deduce that disturbances such as geological events and human activities have occurred in the same way the habitats of the entire range of distribution and the populations only translocate but regionally remain, namely, the populations undergo local distributional shifts but rarely go extinct from an entire region. Geological events, biogeographical processes, ecological conditions, and human activities are all connected to explain the present distribution pattern of this endemic and endangered species. Our results are expected to have positive implications for conservation in the Izta-Popo National Park and especially for the zacatuche populations on Tláloc.
Currently, Romerolagus diazi conservation on the Tláloc Volcano in Sierra Nevada is mainly the result of local actors who have engaged in managing their land favoring the conservation of this emblematic species. Ongoing research on the potential for participatory landscape conservation to engage local actors as allies in conservation tasks is still to be documented (sensuVelázquez et al., 2003). Further research to document if these connected driving forces also explain the distribution of species that are sympatric with the volcano rabbit is yet to be conducted.










nueva página del texto (beta)




