This study was conducted in Don Carmelo Multiple Use Private Reserve, a protected area of about 44,000 ha, located in La Invernada valley (31° 10’S, 69° 46’ W), in San Juan province, Argentina (Figure 1). The study area is located in the southernmost extent of the Puna ecoregion, a desert environment at an elevation of > 3,100 m (Ellenberg, 1979; Matteucci, 2012). Vegetation is composed of low xerophytic shrubs and grasses (i.e., a shrub-steppe environment), with a large amount of bare soil (Márquez, 1999; Lara et al., 2007; Ripoll & Martínez Carretero, 2019). The climate of Puna at this latitude is dry, cold, with a wide daily temperature range, and primarily summer rainfall. The mean annual temperature is 8.15 °C, maximum absolute temperature is 26 °C, and minimum absolute temperature is -22 °C. Precipitation is about 100 mm/year, and snowfall occurs mainly between May and October with thickness up to 50 cm. Snow usually remains on the ground for no longer than 15 days (Andino & Borghi, 2017).

The most abundant herbivores occurring in the reserve are guanacos (Lama guanicoe) and the subterranean rodent Ctenomys mendocinus, locally called tuco-tuco. C. mendocinus are subterranean hystricognath rodents with a mass of 108-253 g, that dwell in areas ranging from 400 m to 3,200 m.a.s.l. in arid and semi-arid environments of western Argentina (Rosi et al., 2005). C. mendocinus are strict herbivores, primarily consuming woody plants, but also feeding on grasses when available (Rosi et al., 2003), creating discrete patches through their intense herbivory, characterized by an increase in bare soil cover and a decrease in grass and shrub cover (Andino & Borghi, 2017, Borghi et al. 2020) (Figure 2). The area of these disturbed patches ranges from around 700 to 20,000 m². C. mendocinus home range depends on sex and locality, ranging from 12 to 43 m² (Rosi et al., 2005). In Don Carmelo Multiple Use Private Reserve, density of tuco-tucos ranges from 3.3 to 11.7 / ha (Borruel, 2013), with this burrowing activity also affecting soil nutrient concentration (Lara et al., 2007).

The study site was located on a high-altitude plateau, situated between the north-south oriented Sierra del Tigre and Sierra de la Invernada mountain ranges. Spanning approximately 4 km in width (west to east) and 15 km in length (north to south), this plateau maintains a consistent elevation ranging from 3,000 to 3,100 m.a.s.l. Such uniformity ensures comparable climatic and environmental conditions across all sampled transects, providing a stable context for our study. Within this expanse, patches are found in various states of disturbance ranging from very old patches, where vegetation is in the process of recovering, to patches where disturbance activity by tuco-tucos has been very recent. For this study, we selected only those patches that were either highly disturbed (but not in a recovery process) or relatively undisturbed (undisturbed) by C. mendocinus. In the disturbed patches, the abundance of burrow holes was approximately 9.6 times greater than in the undisturbed ones, soil mounds were about 7.4 times more prevalent, and plant cover was 55 % less than in undisturbed patches (Table 1).

Fieldwork was conducted from February 2001 to October 2004. In February 2001, sampling was carried out to characterize the plant cover, as well as to quantify the abundance of mounds and burrow holes created by Ctenomys mendocinus. This sampling was performed across two distinct types of areas: those heavily disturbed (disturbed) by the activities of C. mendocinus and those that remained relatively undisturbed (undisturbed). The degree of disturbance by C. mendocinus on the vegetation was assessed visually, ensuring consistency in our categorization of each patch. To accurately describe this influence, we randomly established thirty transects, each measuring 30 m in length, across both disturbed and undisturbed sites. Specifically, fifteen transects were placed in areas heavily disturbed by tuco-tucos, and fifteen were situated in relatively undisturbed patches (Figure 2). For each transect, we collected samples from ten plots, each covering an area of 2 m² and spaced 1 m apart, resulting in a total sampled area of 20 m² per transect. The sampling unit was the transect.

To assess the influence of Ctenomys mendocinus on bird populations, fieldwork was conducted from February 2003 through October 2004. We selected two types of patches: (1) patches highly disturbed by tuco-tucos (with high density of holes [1.73 ± 0.11 /m², mean ± SE] and mounds [0.91 ± 0.03 /m²];= ‘‘disturbed patches’’), and (2) patches relatively undisturbed by Ctenomys (with low density of holes [0.18 ± 0.05 /m²] and mounds [0.12 / ± 0.02 / m²]; = ‘‘undisturbed patches’’) Lara et al. (2007). We randomly selected 109 patches (with a diameter of 50 meters or larger) for this purpose, divided between disturbed (56) and undisturbed (53), ensuring each patch was represented by a single transect. These transects, 50 m in length and 40 m in width, were placed at least 200 m apart to minimize overlap in bird territories. Surveys were carried out during times of peak bird activity - in the early morning (7:30 am to 10:30 am) and late afternoon (4:30 pm to 7:30 pm) - across three seasons: summer (February 2003), autumn (May 2003), and spring (October 2004). Each transect was walked at a steady pace for 30 minutes by the same observer, who recorded sightings of bird species, identified with the help of the field guide by Narosky and Yzurieta (2003). The scientific and common names of the birds observed are listed in Table 2.

To describe and assess the ecological effect of Ctenomys activity on plant cover, burrow openings, and soil mounds, we employed Generalized Linear Models (GLMs) with appropriate error distributions for each type of data. Gaussian error distribution was used for the continuous variable of plant cover, whereas a negative binomial distribution was employed for the count data of burrow holes and soil mounds to accommodate overdispersion. We incorporated disturbance level as a fixed effect in all models to directly evaluate its ecological effect. The statistical analyses were performed in R (R Core Team, 2021), utilizing the functions ‘glm’ for Gaussian models and ‘glm.nb’ for negative binomial models, provided by the ‘stats’ and ‘MASS’ packages (Venables & Ripley, 2002), respectively. Overdispersion was evaluated using the ‘dispersiontest’ function from the ‘AER’ package (Kleiber & Zeileis, 2008). Additionally, we employed the ‘emmeans’ package to estimate and interpret the fixed effects with their confidence intervals, thereby quantifying the magnitude of disturbance’s influence.

We applied Generalized Linear Mixed Models (GLMMs; Zuur et al., 2009) with Poisson error distribution and Negative Binomial error distribution if variables exhibited overdispersion (Ĉ > 1; Crawley, 2007) to analyze if C. mendocinus disturbance affects the abundance of all bird species occurring in the area. The response variable was abundance of all bird species, the explanatory variable was disturbance by Ctenomys (two levels: disturbed and undisturbed patches), and the random effect used was the month of survey (three levels: February, May and October). GLMMs were fitted using the glmer function from the lme4 package (Bates et al., 2015). The significance of fixed and random effects was assessed using the likelihood ratio test (LRT; Bolker et al., 2008).

To determine the impact of Ctenomys disturbances on avian populations, we concentrated on species with total abundances greater than 20 individuals, analyzing data from months with significant recorded abundances. This approach aimed to enhance the statistical robustness of our intra-specific abundance comparisons, focusing on the effects of disturbances on four bird species: Oreopholus ruficollis and Thinocorus rumicivorus in February, and Muscisaxicola cinereus along with Geositta cunicularia in May. We encountered significant numbers of zero counts in our bird species abundance data, exceeding what standard Poisson or Negative Binomial distributions would predict, indicative of excess zeros. The data also exhibited over-dispersion, where variance surpassed the mean, suggesting that these standard models were insufficient for our analysis (Zuur et al., 2012). To accurately assess the effect of Ctenomys disturbances, we employed Zero-Inflated models, which are adept at managing datasets characterized by excess zeros and over-dispersion. These models include two components: a logistic regression model that predicts the likelihood of zero counts, and a count model (Poisson or Negative Binomial) for the positive counts. Depending on the dispersion index (Ĉ > 1; Crawley, 2007), we applied Zero-Inflated Poisson Models (ZIP) and Zero-Inflated Negative Binomial Models (ZINB) using the pscl package (Jackman, 2020) and the lmtest package (Zeileis & Hothorn, 2002). The significance of the disturbance effect as a fixed effect in our models was evaluated using the likelihood ratio test (LRT; Bolker et al., 2008), enabling us to discern the effect of Ctenomys disturbances on bird species abundance.

To evaluate the differences in bird species richness and diversity in areas with varying levels of C. mendocinus activity, as well as to analyze the variations in the rank-abundance curves among these areas, we adopted a comprehensive analytical approach. We analyzed the entire dataset, incorporating all species recorded throughout the study period without omitting any due to low abundance. This comprehensive dataset underpinned the calculation of diversity indices, such as species richness, Shannon diversity, and Simpson diversity. Our approach ensured that the resultant diversity indices accurately reflect the full ecological complexity of both disturbed and undisturbed sites, without bias from species omission. As every patch type (disturbed by Ctenomys and undisturbed) had a different number of transects (56 and 53 respectively), we performed coverage-based rarefaction curves to estimate sampling completeness at each patch type (Chao et al., 2014). Sample completeness refers to the proportion of total abundance or number of individuals in an assemblage that belong to the species represented in the sample (Chao & Jost, 2012).

In order to quantify the species diversity in disturbed and undisturbed sites, ⁰D, ¹D and ²D metrics was used, which are parameterized by the Hill numbers 0, 1 and 2, respectively (Hill, 1973; Jost, 2006). The values of q are referred to as the “order” of the diversity measure. ⁰Dα represents species richness (q = 0); ¹Dα, Shannon diversity index, represents the number of common species (q = 1); and ²Dα, Simpson diversity index, represents the number of dominant species in a community (q = 2) (Jost, 2006). We used the function mcpHill (R package “simboot”) to compare diversities of patches highly disturbed by Ctenomys and undisturbed (Pallmann et al., 2012; Scherer & Pallmann, 2017; Canty & Ripley 2019). We assessed statistical differences in species diversity by applying multiple contrast tests based on the resampling technique (Westfall & Young, 1993).

We characterized species diversity at both types of patches, building two types of integrated rarefaction and extrapolation curves (sample-size and coverage-based), with 95 % confidence intervals, obtained by a bootstrap method based on 1,000 replications using iNEXT package (Hsieh et al., 2016). For sample-size curves, the maximum reference sample size (226 individuals) was used as a base sample size, and for coverage-based curves we standardized the samples to the maximum coverage for reference samples (100 %). Pielou’s evenness index (J) was calculated as a measure of equitability, i.e. the proportion of the diversity observed in each patch was assessed in relation to the maximum expected diversity (Magurran, 2004). These analyses were supplemented with bird species rank abundance curves for both types of patches using the Biodiversity R package (Kindt & Kindt, 2019).

Results are presented as mean ± standard error (SE) and, for null hypothesis testing, statistical tests were considered significant at α < 0.05. All analyses were performed with R Software 3.6.1 (R Core Team 2019). Results are expressed as the untransformed mean ± SE.

We recorded a total of 319 birds at patches disturbed and undisturbed by Ctenomys. They belong to two orders, seven families, eight genera and 11 species (Table 2).

Birds (all species) were significantly less abundant (93 individuals, mean number of birds per transect: 1.66 ± 0.93) at disturbed patches with high levels of activity by Ctenomys compared to undisturbed patches (226 individuals, 4.26 ± 0.88 birds/transect) (GLMM, LRTtest χ2 = 9.41; df = 1; p = 0.0021). We found that total abundance was also influenced by the random factor (month; p < 0.0001). The most abundant recorded bird species included the insectivores Muscisaxicola cinereus, Geositta cunicularia, Oreophollus ruficallis and the granivore Thinocorus rumicivorus. Since these species exhibited a migratory abundance pattern, we compared their abundances between disturbed and undisturbed sites during periods when their numbers were sufficient for statistical comparison. The most frequently recorded bird species was M. cinereus, with observations of 62 individuals in undisturbed patches versus 50 in disturbed patches during May. Following this species were G. cunicularia in May, and T. rumicivorus and O. ruficallis in February, with each species’ populations exceeding 19 individuals in a given month. The abundance of M. cinereus and T. rumicivorus was significantly different between undisturbed and disturbed patches by Ctenomys, but not that of O. ruficallis or G. cunicularia (Table 3). M. cinereus was more abundant at undisturbed than disturbed patches (mean 6.20 individuals/transect vs 3.85 respectively), and T. rumicivorus followed the same pattern (0.47 vs 0.06 respectively). Ctenomys did not affect the probability of absence of these four species (Zero-inflated models; Table 3). Mimus patagonicus, Geospizopsis plebejus, Geositta rufipennis, Asthenes dorbignyi and Sicalis auriventris were considered rare species because they were observed only in undisturbed patches and with low abundance (abundance < 4 individuals in all transects). For some cases, for instance for G. rufipennis, we observed only one individual over all study periods; and in the case of S. auriventris, we recorded three individuals just once on one transect.

At undisturbed patches, we recorded 226 individuals among 11 species, and at disturbed patches we recorded 93 individuals among 6 species. Sampling completeness reached 99 % for undisturbed patches and 100 % for disturbed ones.

Using Hill numbers (i.e., the effective number of species), our evidence confirmed that undisturbed patches were significantly more diverse than those disturbed by Ctemomys (⁰D, p = 0.01; ¹D, p = 0.04; ²D, p = 0.05) (Table 4 and Figure 3). At patches undisturbed we found one species represented by one individual (singleton species) and one species represented by two individuals (doubleton species). However, at disturbed patches we did not find any singletons or doubletons. The sample-size rarefaction and extrapolation curves (base sample size= 226) showed a consistent pattern in all indexes, with the diversity curve for undisturbed patches above the curve of disturbed patches for all three parameters (Figure 3).

By standardizing both patches to the same sampling effort (sample coverage 1.00), we found that, at both types of patches, the coverage-based curves give the same ordering of diversity (⁰D > ¹D > ²D) with higher diversity at patches without disturbance than at disturbed patches (Figure 3). Pielou’s evenness shows that patches free from disturbance by Ctenomys were more similarly distributed than disturbed patches (J’ = 0.77 and J’ = 0.68, respectively).

The rank-abundance curve showed that Muscisaxicola cinereus was the most abundant bird species at both types of patches. At undisturbed patches, Geossita rufipennis had the lowest abundance, while Thinocorus rumicivorus had the lowest abundance at disturbed ones (Figure 4). The declining species-rank curve indicates that there is great dominance of the most abundant species at both types of patches. A. dorbingnyi was recorded only at undisturbed patches, as well as the least abundant species S. auriventris, M. patagonicus and G. rufipennis (Figure 4).