1. Droughts are projected to increase in magnitude, frequency and duration in the near future. In rangelands, the provision of valuable ecosystem services such as forage supply for livestock productivity is intimately linked to rainfall patterns, which makes it particularly vulnerable to droughts. Nonetheless, rangelands can differ in their sensitivity to droughts as shown by strong differences in the impacts of inter-annual precipitation changes on vegetation productivity in different sites. The aim of this study was to assess the sensitivity to droughts of nine rangelands located across a broad aridity gradient in Argentina, South America.

2. We experimentally imposed comparable droughts under field conditions by reducing a fixed proportion of each incoming precipitation event within-year during three consecutive years and tracked changes in total aboveground and forage productivity.

3. We found that arid and semi-arid rangelands were more severely impaired in their forage provision by drought than mesic rangelands, i. e. that sensitivity to drought declined as aridity decreased. Forage productivity decreased on average by ca. 50%, in arid and semi-arid rangelands, whereas mesic sites did not exhibit significant changes between drought and control treatments. The negative impact in forage productivity of arid and semi-arid rangelands was mainly driven by the productivity reduction of few key plant species at each site. In seven of the nine rangelands, we found detrimental effects of drought on forage productivity during the first experimental-drought year, and in five of them the impact was further accentuated until the end of the experiment, which indicates how serious can these events be.

4. Synthesis and applications: Our main findings indicate that the drought-induced impacts on forage provision are higher as aridity increases. This pattern highlights the urgent need to implement strategies to mitigate the detrimental consequences of drought, particularly in arid and semiarid rangelands, where forage provision is strongly associated with human well-being. Management approaches focused on key forage species, such as reducing the grazing pressure during drought periods according to these species’ productivity dynamics, can attenuate impacts on vulnerable ecosystems, preserving the rangelands’ integrity while maintaining high long-term productivity levels.

Aridity gradient

In order to evaluate the effect of drought on forage productivity, we conducted a coordinated drought experiment (with the appropriate permits and licenses for fieldwork) in nine rangelands scattered across a natural aridity gradient (Figure 1), with mean annual precipitation (MAP) ranging from 170 to 950 mm yr^-1 and from 5 to 21°C mean annual temperature (Table 1).  Our study did not require for ethical approval. We studied natural rangelands, typically subjected to extensive livestock grazing, including grass steppes, grass-shrub steppes, and highly productive grasslands that encompass a large range in  plant species diversity, ANPP, forage productivity (FP), soil type and texture, land-use history, and livestock type (Table 2; Table S4; Figure S1) (Oyarzabal et al., 2018; Peri et al., 2021). Considering the aridity index (MAP/Potential Evapotranspiration; the lower, the drier) the nine rangelands were equally distributed into arid, semi-arid, and mesic ecosystems with two of them which possessed an aridity index value edge of the category classified as “arid, semi-arid” (Le Houérou, 1996; Table 1).

Experiment set-up and plot design

All the experimental sites are part of an international collaborative research network consisting of coordinated distributed drought experiments, The Drought Network (IDE; http://wp.natsci.colostate.edu/droughtnet/). The sites followed a common experimental protocol which allows reliable comparisons among contrasting ecosystems (Knapp et al., 2017) and included experimental rainfall manipulations to understand the ecological drought impacts (Hoover & Rogers, 2016). The simulation of  droughts was induced through a passive well-tested design based on rainout shelters that intercept a fraction of the incoming precipitation (Yahdjian & Sala, 2002), a fraction which differ among sites to emulate a drought that occurred only once in the last 100 years (Knapp et al., 2017). To determine the percentage of precipitation to intercept, we applied the standardized protocol disposed by The Drought Network and ran the Precipitation Manipulation Tool. This is a software that allows to upload the longest precipitation time series of each site (https://drought-net.colostate.edu/terrestrial-precipitation-analysis-package; Lemoine et al., 2016). Antecedent conditions did vary as sites differed in their previous pre-treatment precipitation (Table S1). Thus, we reduced annual precipitation from 60% in arid sites to 50% in mesic ones (Table S1). Intercepted precipitation was collected in gutters and directed away by pipes in order to avoid water infiltrating into the experimental drought plots. The rainout shelters were constructed with transparent plastic tiles, placed above the plant canopy (1.20 to 1.60 m height) in order to minimize impacts on micrometeorological conditions (Figure 1 see photos). Except for a decrease of 11 to 20% of incident radiation (compared to the measurement outside the shelter), which are relatively low compared to other materials and did not imply significant reductions for plant canopy light interception (Yahdjian & Sala, 2002), no other side effects were detected.  The rainout-shelter design employed in this study has proven to accomplish the desired reduction in water input and the expected effects on soil moisture (Gherardi & Sala, 2013; Yahdjian & Sala, 2002), and a considerable number of studies have set up drought experiments that imposed passive reductions in precipitation using this rainfall shelters design with similar results (see for example Byrne, Adler, & Lauenroth, 2017; Knapp et al., 2002; Siebert et al., 2019). The experimental plots (of at least 3 m × 3 m) were installed in autumn on different years (2015 to 2017) following a completely randomized block design that included control plots of the same size. At each experimental site, the experiments were generally replicated three times (n = 3), although some sites had four replicates (n = 4). All experimental sites were fenced to prevent grazing on experimental plots.

Data collection

We took pre-treatment measurements (Time 0) describing the plant community, soil and long-term plant species composition and abundance. Also, ANPP and soil properties of each site were collected to be considered as predictors of drought sensitivity (Table 2; Table S2, S3 and S4). During the experiment, we annually estimated plant cover of each species in permanent plots of 1 m² in the center area of each treatment plot, and aboveground biomass at the peak of the growing season, inside a core-sampling area in each drought or control plot, excluding a buffer area to minimize edge effects. For the Río Mayo and Chacra Patagones sites, the estimation of ANPP was made non-destructively (Table 2), through allometric equations that related the relative cover of each plant species present with their aboveground biomass (e.g. Flombaum & Sala, 2007). For the rest of the sites, ANPP was estimated by clipping to the soil surface all aboveground biomass within two frames of 0.2 m x 0.5 m during the peak of biomass productivity (Sala & Austin, 2000), preventing to repeat clipping the same area each year. The harvested material was sorted into green, senescent (yellowish) and dead components for the main plant functional groups: grasses, forbs, legumes and shrubs; then dried to constant weight in a stove for 48 h to obtain dry matter. We repeated plant measurements during three consecutive years in drought and control plots, except in San Claudio, where we failed to harvest during the third experimental year because of a national COVID-19 quarantine and in Chacra Patagones during the first year because of logistic complications. For all sites, we calculated forage productivity (FP) following Easdale & Aguiar (2012): we multiplied the plant cover of each species by the ANPP of the corresponding plant functional group (grasses, forbs, legumes, shrubs) to obtain species ANPP (ANPP_SPP), and then we multiplied ANPP_SPP by a forage aptitude factor (unitless, detailed for each species in each site in Table S2) based on Easdale & Aguiar (2012) guidelines and considering the literature (Ambrosino et al., 2021; Guevara et al., 2002; Ojeda et al., 2018; Oñatibia et al., 2015). We then asked for the approval of local experts of each site for the classification of the forage aptitude factors. The forage aptitude factor classified plant species according to three categories of consumption of grazers: highly consumed (forage aptitude factor=1), moderately consumed (forage aptitude factor=0.5) and avoided by grazers (forage aptitude factor=0). Then, FP values for each site and year were used to assess drought sensitivity using an absolute metric, which has been used in previous studies (Koerner et al., 2015; Raynor et al., 2020; Smith et al., 2017; Moran et al., 2014; Wilcox et al., 2017; Wilcox et al., 2015) and a relative index, following Equations 1 and 2, respectively.

Data Analysis

We analyzed the Data from all sites in R software (RTeam, 2019). We used linear regressions models (lm function) to describe changes of FP and ANPP along the aridity gradient and to fit the data of the proportion of FP to ANPP for each site of the gradient. We used linear mixed-effects models with separate repeated measures analysis of variance (ANOVA) to evaluate differences between treatments (Drought and Control) for FP for each site, with year, treatment, and their interaction as fixed effects, while sampling block was considered as a random
effect (random intercept). Pre-treatment values were incorporated and considered in the analysis using the offset function which corrects for initial differences. To analyze FP sensitivity to drought along the aridity gradient, we tested linear and non-linear regressions for relative and absolute sensitivity and then selected the best-fitted models for each response variable based on the Akaike Criterion. In particular, for the absolute sensitivity index, we concluded that the non-linear model was the best model based on the AIC (AIC linear model: 21.2 vs AIC non-linear model -50.86), which takes into account the trade-off between fit (residual error) and parsimony (i.e., model complexity in terms of number of parameters), and the determination coefficient (R² linear model: 0.32 vs R² non- linear model: 0.51). We used the lme function and nmle packages (Bates et al., 2015). Models followed the assumptions of homocedasticity of variances and normal distribution of residuals, assessed through Levene and Shapiro–Wilks tests, respectively.

Refer to the README files for Data supporting each figure