Factors affecting the distribution and abundance of the Endangered volcano rabbit Romerolagus diazi on the Iztaccihuatl volcano, Mexico

Abstract

The volcano rabbit Romerolagus diazi is an Endangered species endemic to Mexico, with a range of < 400 km². We investigated threats from destruction, fragmentation and degradation of habitat, hunting, and cattle grazing intensity in relation to the distribution and abundance of the volcano rabbit on the Iztaccihuatl volcano. Faecal pellet counts were taken as a proxy for rabbit abundance in 1,718 random 0.2 m² quadrats at 859 sampling points along 25 transects, covering an area of c. 100 km² at altitudes of 3,400–4,000 m. Presence of the species was significantly associated with absence of closed forest, absence of long grass types (not bunchgrass), shallow inclines, absence of cattle grazing, lower altitude, low hunting pressure (measured by proximity to ranger station), absence of bare ground and, contrary to previous findings, increased frequency of fire. The species was significantly more abundant in habitats with a greater percentage cover of zacaton (bunchgrass) and short grass types. It was significantly less abundant in areas with more hunting (measured by proximity to ranger station) and cattle grazing. Key conservation priorities are therefore the protection of the subalpine zacaton bunchgrass-dominated habitat type, strict enforcement of hunting laws and the removal of livestock from relevant national park boundary areas. However, the results suggest that frequent fires have a significant positive effect on the occurrence of the volcano rabbit as a result of habitat improvement and this is often a consequence of anthropogenic management of land for cattle grazing.

Introduction

The volcano rabbit Romerolagus diazi is categorized as Endangered on the IUCN Red List (AMCELA et al., 2008). It is endemic to the central Trans-Mexican Neovolcanic Belt and has one of the most restricted ranges of any mammal (Granados, 1979; Thornback & Jenkins, 1984; Bell et al., 1985). There are a number of distinct populations of the species in the uplands surrounding Mexico City, distributed in discontinuous patches on the volcanoes of Popocatepetl, Iztaccihuatl, El Pelado and Tlaloc (including Sierra Chichinautzin; Fa & Bell, 1990). Its range covers c. 386 km² (Velázquez, 1994) and the species has disappeared from parts of its historical range in the central Trans-Mexican Neovolcanic Belt, such as the Nevado de Toluca (Cervantes et al., 1990; Fa & Bell, 1990). The most recent population estimate was c. 7,085 (Portales et al., 1997). R. diazi is a subalpine habitat specialist (Velázquez & Heil, 1996) and occurs at its highest density in zacaton bunchgrass communities dominated by the genera Festuca, Calamagrostis, Muhlenbergia, Agrostis and Stipa (Velázquez, 1993). The greatest area of suitable habitat remaining is within the boundaries of the Iztaccihuatl-Popocatepetl National Park (henceforth Izta–Popo), at 3,400–4,200 m.

The volcano rabbit's habitat is shrinking, mainly as a result of destruction and degradation (Velázquez, 1993) caused by cattle grazing, encroachment of agriculture and human settlements, logging, harvesting of zacaton and anthropogenic fire (Fa & Bell, 1990; Portales et al., 1997). Its range is becoming increasingly fragmented as a result of habitat loss and highway construction, causing the fragmented populations to become genetically isolated and increasing their risk of local extinction from random processes (Velázquez, 1993; Frankham, 1997). However, genetic variation in R. diazi is similar to that of more widespread and less fragmented genera of Mexican rabbits, such as Sylvilagus and Lepus (Cervantes et al., 2002). The species occurs in a number of protected areas, principally the Izta–Popo National Park, but hunting, cattle grazing and grass burning persist within the park boundaries (Cervantes et al., 1990; Fa & Bell, 1990; Velázquez, 1993). Although the volcano rabbit is protected under Mexican law there is a lack of local knowledge of its protected status and a lack of enforcement of the law; consequently hunting could cause further declines of the species (Fa & Bell, 1990; Portales et al., 1997).

Here we focus on the slopes of the volcano Iztaccihuatl because it is one of the larger and potentially better protected areas of habitat. We consider the factors affecting populations of the volcano rabbit at two spatial scales. At the larger scale of the entire volcano we aim to determine why the species is present in some habitats and not in others. On the smaller scale we aim to determine what factors predict the abundance of the species within the habitats where it is present. Our objective is to derive conservation recommendations for the volcano rabbit (Portales et al., 1997) by assessing the importance of vegetation, hunting, cattle grazing intensity, fire frequency and altitude in determining its presence or absence and its abundance in areas where it occurs.

Study area

The study was conducted during July–September 2011 on the volcano Iztaccihuatl, located between the states of Mexico and Puebla, covering c. 60 km² on the western slopes and 40 km² on the eastern slopes (Fig. 1). It was not possible to include areas on the volcano Popocatepetl because it was active at the time of the study, with an exclusion zone of 12 km from the crater. We chose the eastern and western slopes of Iztaccihuatl because they allowed us to study the effects of a rapid change in altitude. The study was carried out during the rainy season, which is the breeding season for volcano rabbits (Cervantes et al., 1990).

Fig. 1 (a) The location of the study site in Mexico, where the white rectangle shows the location of (b) within the central Transverse Neovolcanic Belt. (b) The four sites where the volcano rabbit Romerolagus diazi occurs near Mexico City: 1, El Pelado; 2, Tlaloc; 3, Iztaccihuatl; 4, Popocatepetl. The grey rectangle indicates the location of (c), the study area on Iztaccihuatl. (c) The white area indicates the extent of the forest/grassland habitat where the volcano rabbit occurs on Iztaccihuatl and includes the transect locations. The central grey areas mark the largely unvegetated summit of the volcano; the peripheral grey areas are farmland and urban areas.

Methods

We used faecal pellets as a proxy for both species presence and abundance, as validated by Velázquez (1994). Faecal pellet counts were used because R. diazi is crepuscular and shy (Cervantes et al., 1990), time of day does not affect pellet count data, and counts are temporally independent. We considered degradation effects on pellets to be of minor significance. Pellet counts were used as a proxy for abundance only in habitats where the species occurred over a small scale (kilometres), where climatic factors affecting degradation rates were likely to be similar. Furthermore, pellets were not used to ascertain population size or density but only to investigate variations in abundance. Martinez-García et al. (2011) analysed the change in defecation rate with vegetation for R. diazi and concluded that faecal pellets can be used as a reliable index for this species.

The total sampling effort was 16 days (150 hours) of data collection, with an additional 3 days (20 hours) of experimental study design in the field. The survey comprised 15 transects on the western slopes and 10 on the eastern slopes of Iztaccihuatl, parallel to the altitudinal gradient, which is approximately east–west (Fig. 1). The survey area was c. 15 × 15 km. On each transect we took point samples at every 20-m increment of altitude from 3,400 m to 4,000 m. The number of point samples per transect varied with the topography but was generally 30–35 (mean 37.2 ± SE 2.3 points per transect, n = 25 transects). Transects covered c. 4.1 ± SE 0.2 km (straight line distance from start to finish) and the mean distance between sample points was 111.2 m. The mean distance between point samples containing pellets on a transect was 223.0 ± SE 13.4 m (n = 266), with 99 (37.3%) quadrats separated by ⩾ 223 m and 35 (13.1%) separated by ⩾ 334 m. Most distances between samples probably exceed the mean range size (the pygmy rabbit Brachylagus idahoensis, which is slightly smaller than the volcano rabbit, has a range of < 0.3 ha; Katzner & Parker, 1997), and therefore point samples can be considered to be reasonably independent. The arrangement of sampling points along transects was for logistical reasons; random sampling would have yielded the same results because of the short duration of the study and the fact we counted pellets, which do not move. Nevertheless, we explore the effects of transect sampling on the model results.

At each sampling point we took two 0.2 × 0.2 m quadrat samples, 5 m apart and perpendicular to the transect line. Faecal pellets were counted in each quadrat and pooled into one point sample. In total 1,718 quadrat samples were taken in pairs at 859 points. At each sampling point we recorded vegetation abundance in terms of percentage cover within a 10 m radius, for zacaton bunchgrass, Salvia spp., Lupinus spp., long grass types (not zacaton bunchgrass), short grass types, scrub bushes and small trees, and closed pine forest. Other variables recorded at each point included the number of mature and juvenile pines, latitude & longitude (position), altitude, incline, cattle grazing intensity, fire frequency and transect number (Table 1).

Table 1 The variables recorded at each of 859 sampling points on Iztaccihuatl volcano, Mexico (Fig. 1).

Data selection was carried out using PASW Statistics v. 18 (SPSS, Chicago, USA). We used northing and easting coordinates to calculate the relative distance of all point samples from Alsomoni station, the only large ranger station on Iztaccihuatl. We used this variable (distance to station Alsomoni) to attempt to account for illegal hunting pressure, the hypothesis being that less hunting occurs close to the station. Hunters were seen on several occasions during fieldwork and none of these sightings were proximal to the ranger station but nevertheless this spatial variable was not explicitly linked to hunting pressure. We used R v. 2.13.0 (R Development Core Team, 2012) for data analysis. Six point samples had missing data as a result of equipment malfunction and these were removed, leaving 853 point samples. The first analysis examined the factors determining the occurrence or distribution of the species on a large spatial scale and the second analysis examined the factors determining the species’ abundance.

For both analyses we constructed a full general linear model of all possible variables (Table 1) affecting the distribution of the species.

For the first analysis a binomial regression model was used to analyse presence/absence: pellets were present in 266 point samples and absent in 587. For the second analysis only point samples where the species was present (i.e. samples with faecal pellets, n = 266) were included, to determine factors affecting the abundance of the species. Pellet count was natural log transformed, facilitating a reasonable fit of a quasi-Poisson distribution.

The full models were then simplified using the dredge command in R and parameter estimates were averaged for all models within two AIC units of the best model to represent the best, statistically similar models according to the information theoretic approach (Akaike, 1974). The top models were then evaluated for both biological and statistical significance. Residuals of the top models met model assumptions reasonably well, as demonstrated by the plot command in R and according to criteria in Crawley (2007). However, any one of the top models obtained using the dredge command, or the mean parameter estimates, could have been used, giving similar results in relation to statistical significance and biological effects.

In both analyses the variables mature pines, juvenile pines, position and transect were excluded from the full models to remove autocorrelation and to simplify analysis within biological reason. Mature pines and juvenile pines were excluded on the basis that they measured the same component of the vegetation as the percentage of closed forest. Position was excluded because spatial considerations were incorporated in the variables distance to station Alsomoni and altitude. Transect was included in top models as a random effect (using generalized linear mixed models in R, lmer), even though the clustering of samples along transects was unlikely to have biased pellet counts in any systematic way (transects were evenly distributed around the volcano). The results for the random effects models indicated that there was some significant variance accounted for by some transects, and therefore there was additional spatial variation in rabbit distribution that we were not accounting for that was captured by the location of points along a transect. All other variables retained in the model and their biological and statistical significance remained largely the same regardless of whether transect was included (AIC mixed model = 878.7, AIC without transect = 879.2). We therefore removed transect from all subsequent models. Where biologically sensible (e.g. altitude, because there should be an optimal range of altitudes for the rabbits rather than a simple linear function), quadratic functions were also considered by including the square of variables in both full and top models. This did not improve any of the models, and therefore quadratics were not retained in any models.

Results

The presence/absence analysis resulted in 15 models within 2 AIC of the top model (Table 2) and the best model retained seven significant variables (Table 3). The probability of occurrence of the species decreased with increasing closed pine forest cover, cattle grazing intensity, long grass cover, incline, hunting pressure, and altitude; the probability of occurrence increased with increasing fire frequency (Fig. 2; Table 3). As the percentage cover of closed pine forest increased from 0 to 100% the probability of the species’ occurrence decreased by 97.1%. As the cattle grazing intensity increased from low to high the probability of the species’ occurrence decreased by 57.9%. As the percentage cover of non-zacaton long grass increased from 0 to 80% the probability of the species’ occurrence decreased by 82.4%. As the incline of the site increased from shallow to steep the probability of the species’ occurrence decreased by 78.5%. As the distance from the ranger station Almosoni increased from 0 to 12 km the probability of the species’ occurrence decreased by 44.4%. As altitude increased from 3,400 to 4,000 m the probability of the species’ occurrence decreased by 48.3%. As fire frequency increased from low to high the probability of the species’ occurrence also increased by 39.3%. The percentage cover of scrub was included in some top models but had no significant effect on the occurrence of the volcano rabbit in the top model.

Fig. 2 Sunflower plots showing the significant biological effects of the variables investigated on the distribution of the volcano rabbit on the volcano Iztaccihuatl, using data from 853 random point samples and using faecal pellets as a proxy for the species’ presence. Plots (a) to (i) are in order of biological significance (Table 3). The raw data are plotted as ‘petals’, which represent the number of sample points, with presence points at the top of the graph and absence points at the bottom. Lines represent the predicted values from an optimal binomial regression (general linear model), using median values for all other variables, solid lines represent mean predicted values and dashed lines 95% confidence limits.

Table 2 Model variables for presence/absence analysis for the volcano rabbit Romerolagus diazi on the Iztaccihuatl volcano (Fig. 1), with mean estimates and relative weighting for all 15 models within 2 AIC of the best model.

* Significant predictor variables

Table 3 The results of the top binomial regression general linear model describing the presence/absence of the volcano rabbit, with model variables, estimates, P values and mean % change in occurrence ([maximum predicted value – minimum predicted value]/minimum value, i.e. the percentage change of a variable relative to its minimum value).

* Shown as a measure of biological significance, using median values for all other parameters

The abundance analysis resulted in 10 models within 2 AIC of the top model (Table 4) and the best model retained four significant variables (Table 5). In general the species was more common as the percentage of zacaton and short grass increased and less common as grazing intensity and hunting pressure increased, assuming that hunting increases with distance from ranger station Alsomoni (Fig. 3; Table 5). As the distance from ranger station Almosoni increased from 0 to 12 km the species’ abundance decreased by 72.2%. As the percentage cover of zacaton bunchgrass increased from 0 to 95% the species’ abundance increased by 124.2%. As cattle grazing intensity increased from low to high the species’ abundance decreased by 46.6%. As the percentage cover of short grass increased from 0 to 70% the species’ abundance increased by 220.6%. The percentage cover of scrub and percentage cover of Lupinus spp. were also retained in the top model but they had no significant effects on the abundance of the species.

Fig. 3 Scatterplots showing the significant biological effects of the variables investigated on the abundance of the volcano rabbit on Iztaccihuatl, using data from the 266 random point samples where faecal pellets were found, with faecal pellet counts as a proxy for the species' presence. Plots (a) to (f) are in order of biological significance (Table 5). Each point represents a single point sample, with associated values of log transformed count data. Lines represent the predicted values from an optimal quasi-Poisson regression (general linear model) using median values for all other variables, solid lines represent mean predicted values and dashed lines 95% confidence limits.

Table 4 Model parameters for the abundance analysis for the volcano rabbit, with their mean estimates and relative weighting for all 10 models within 2 AIC of the best model.

* Significant predictor variables

Table 5 The results of the top model for the abundance analysis for the volcano rabbit, with model variables, estimates, P values and mean % change in abundance ([maximum predicted value – minimum predicted value]/minimum value, i.e. the percentage change of a variable relative to its minimum value).

* Shown as a measure of biological significance, using median values for all other parameters

Discussion

Our results indicate that volcano rabbits are present where closed forest, long grass, cattle grazing, hunting, bare ground and incline occur at low levels and where there is a higher frequency of fires, and their abundance is higher in areas with higher percentage cover of zacaton bunchgrass and short grass types and lower levels of hunting and cattle grazing. It is important to note, however, that our measure of presence and abundance is an index based on faecal pellet counts and may have been subject to false negatives as a result of degradation, only measuring the effects of recent fire, and basing our estimate of hunting pressure on an untested assumption that there was proportionally less hunting closer to where park rangers were based. Nevertheless we believe that these potential problems do not introduce a systematic bias and are likely to render our hypothesis tests more conservative, thus strengthening the significance of any effects we identified.

Closed forest dominated by Pinus hartwegii, also known as highland temperate forest (Galicia & García-Romero, 2007), has a strong biological effect because as the habitat changes from open to completely closed forest cover the likelihood of occurrence of the volcano rabbit decreases to almost zero. This supports the current understanding of the dependence of the species on the subalpine zacaton grassland habitat found above the forest edge (Velázquez & Heil, 1996). In many cases the transition zone between highland temperate forest and subalpine grassland habitats can occur as low as 3,400 m and this may be dependent on anthropogenic effects on the forest. Although the abiotic altitudinal upper limit to the species’ range is intrinsic, its occurrence at lower levels is constrained by forest occurrence, which can be managed.

Another vegetation variable that has a significant effect on the distribution of the species is the percentage cover of non-zacaton long grass. The probability of the species’ occurrence decreased by 82.4% as the cover of long grass types increased from 0 to 80% (Fig. 2). Although this may seem like a large biological effect, the majority (91%) of sample points had < 20% cover of long grass types and the mean cover in the subalpine habitat above the forest edge was only 4.5 ± SE 0.41% (n = 554). However, in the transition zone between this and subalpine habitats (3,400–3,600 m) percentage cover increased dramatically and then decreased with increasing altitude to < 1% at 4,000 m (Fig. 4). Therefore, the effect of this vegetation component will follow this pattern; it will have a greater negative effect on species occurrence in the transition zone just above the forest edge, especially below 3,600 m, and the effect will decrease with increasing altitude.

Fig. 4 Mean percentage cover of long grasses at altitude 3,200–4,000 m. Data are separated into 100 m intervals to show change in abundance with altitude. There is a transition zone between highland temperate forest and subalpine grassland at 3,400–3,600 m.

The scientific consensus is that fire frequency has a negative effect on R. diazi (Hoth et al., 1987; Velazquez, 1989; Fa et al., 1992; Velázquez, 1994). However, we found the opposite, at least for recent fires. In support of our result, it is recognized that generally bunchgrasses are adapted to and encourage a frequent fire regime to remove competition and permit new growth (Ellsworth & Kauffman, 2010). This suggests that fire could be important in the conservation of the subalpine habitat through the promotion of the dominant zacaton bunchgrass, and at the very least a policy of active fire suppression is unwarranted (Keating, 2007). Zacaton bunchgrass is the main food source for R. diazi (Cervantes et al., 1990; Cervantes & Martinez, 1992), and therefore maintaining its dominance through a frequent fire regime while promoting new growth should benefit R. diazi in turn. However, the effects of fire frequency need to be considered in conjunction with cattle grazing because fire is used specifically to promote new pasture growth for cattle grazing (Fa & Bell, 1990; Portales et al., 1997).

An increase in the intensity of cattle grazing had strong negative effects on both the distribution and abundance of R. diazi. These results agree with the findings of other researchers (Hoth et al., 1987; Velazquez, 1989; Fa & Bell, 1990; Fa et al., 1992), and therefore we reiterate the recommendation that cattle should be removed from the subalpine habitat, especially within park boundaries. However, it is possible that the negative effect of cattle grazing has been overestimated in studies that have not considered the potential positive effect of increased fire frequency as a result of pastoralism. If cattle are removed from the National Park the positive effects of fire frequency may also be reduced because 98% of fires are deliberately started by people to promote grazing. However, this does not mean that pastoralism should continue because the negative effects of cattle grazing probably outweigh the positive effects of the associated fire frequency: fire resulted in a 39.3% increase in the occurrence of R. diazi over the range of the study whereas cattle grazing resulted in a 57.9% decrease over the same range. The removal of cattle and the restoration of native bunchgrass species could increase the frequency of natural fires without anthropogenic intervention, with a build-up of plant debris  as a result of reduced grazing pressure.

An increase in hunting pressure also had negative effects on both the distribution and abundance of R. diazi. We estimated hunting pressure using the distance from the ranger station Alsomoni on the basis that there would be less hunting of the species, which is illegal under Mexican law, in the proximity of enforcers. Consistent with this assumption was the observation of a marked change in abundance at c. 3 km from the ranger station (Fig. 3), coinciding with the limit of visibility and reduced audibility from the ranger station caused by a large east–west ridge. Previous studies have also identified the negative effect of hunting pressure on volcano rabbits (Hoth et al., 1987; Velazquez, 1989; Fa & Bell, 1990; Fa et al., 1992).

Abundance was also affected by the percentage cover of zacaton bunch grass and short grass, which is consistent with zacaton bunchgrass being the rabbits’ main dietary resource as well as providing cover (Cervantes et al., 1990; Cervantes & Martinez, 1992; Rizo-Aguilar et al., 2014). From the model the effect of short grasses appears to be more important biologically but in the subalpine habitat of R. diazi the short grass type has a mean cover of 3.7 ± SE 0.45% (n = 554) and zacaton bunchgrass has a mean cover of 57.3 ± SE 0.97% (n = 554). Furthermore, there are no short grass types present in 80% of the habitat (Fig. 3). Therefore although zacaton bunchgrass is less biologically significant over the full range of the study its effect on the abundance of R. diazi is greater because it is the dominant grass type. Thus conservation effort in the subalpine habitat should be focused on maximizing zacaton abundance.

Recommendations for conservation

Habitat destruction is perhaps the greatest threat facing the volcano rabbit, although it is less so in the Izta–Popo National Park, which has been a protected area since 1935. There has been a low rate of habitat loss (0.35% per decade) in the Park compared with unprotected areas in Mexico (3% per decade; Ochoa-Gaona & González-Espinosa, 2000; Galicia & García-Romero, 2007). The Park boundary is well placed, at 3,600 m, to protect the species in its current range but if the subalpine habitat were put at risk this would have a greater effect on the abundance of R. diazi than either hunting pressure or cattle grazing. Therefore protecting the largest area of habitat, the Izta–Popo National Park, is key to this species’ survival.

Habitat creation could also be considered. Rabbits may become more abundant if the availability of short grass increases and at higher altitude; clearings in upper-level forest could therefore increase the species’ abundance. We recommend further research on the species’ abundance and distribution in highland temperate forest clearings to ascertain how creating such areas could potentially increase populations of the species.

Another priority should be the enforcement of the existing laws prohibiting hunting of the species (Portales et al., 1997). Education and awareness programmes are necessary (Portales et al., 1997), as local people are often not aware of the species’ existence or its protected status. The establishment of clear signs of vigilance by the Park authorities should also be considered. For example, a station could be established on an outcrop near the flat plateau above San Rafael, relatively far from the existing station at Alsomoni, to maximize the area of ground that would be protected.

R. diazi is endemic to a limited area of < 400 km², which is one of the smallest mammalian ranges in Mexico. Its largest contiguous area of habitat is in the Izta–Popo National Park, on the volcanoes Iztaccihuatl and Popocatepetl, which together account for just 0.02% of Mexico's land area but contain > 10% of all Mexican mammalian species (Velázquez et al., 2001). The Transverse Neovolcanic Belt, within which the Izta–Popo National Park lies, also has a high level of endemism, with 52% of Mexico's endemic mammalian species and all of its endemic mammalian genera represented (Koleff et al., 2004). Therefore it is likely that protecting the habitat of one indicator species, in this case the volcano rabbit, would also benefit a number of other species.