MSAVI2 and NDVI
The MSAVI2 and NDVI were used as a proxy for ANPP at a larger spatial scale than the ANPP estimated at the plant or quadrat scale. Five field NDVI spectral reflectance sensors (SRS sensors, Decagon Inc, Pullman, WA, USA) were used to estimate NDVI. An upward-looking hemispherical sensor provided reference values of sky irradiance against which the canopy radiance values of the four downward-looking sensors (one sensor per treatment) were normalized. The SRS sensors were positioned above the plots at a height of approximately 9 m above the ground, covering an area of about 80.6 m2. The downward-looking hemispherical sensors were installed facing the canopy at a 45° angle. A data logger (Em50, Decagon Inc, Pullman, WA, USA) recorded data from each sensor every 1 min and expressed as 30-min averages. We calculated NDVI as:
NDVI = (ρNIR - ρred) / (ρNIR + ρred) = (Rn/In - Rr/Ir) / (Rn/In + Rr/Ir)
where ρNIR and ρred are the percent reflectance in the near infrared (NIR) and red respectively, Rn and Rr are the radiation reflected by the canopy in the NIR and red, and In and Ir are the incident radiation in the NIR and red, respectively. We assumed that the percent reflectance is the ratio of reflected radiation to incident radiation in the specified waveband.
We calculated daily averages, using values ​​acquired only during the noon hour, to reduce data variability due to changes in sun-sensor-surface illumination geometry throughout the day. The data were recorded from July 2019 to July 2020. A technical problem with the downward-looking hemispherical sensor placed in the +N plot prevented obtaining the data for this treatment.
The Copernicus Sentinel-2 mission consists of two polar-orbiting satellites providing multi-spectral imaginery with high resolution (10 m to 60 m) and revisiting time (5 days). Twenty-three Level 1C from July 2019 to June 2020 (image dates are given in Fig. S2) corresponding to tile T19GCK were obtained from the Copernicus Open Access Hub. Images were visually inspected to avoid cloud coverage over the study area. Level 1C products consists of atmospherically corrected TOA (top-of-atmosphere) reflectance images, projected in the UTM WGS84 system and Images converted to BOA (bottom-of-atmosphere).
For each image, the MSAVI version 2 index (MSAVI2) was estimated (Qi et al., 1994). This index is a better alternative to widely used NDVI for arid and semiarid environments as it reduces the effect of soil signal when vegetation cover is sparse.
MSAVI2 =\(\frac{\left[2\ \times\ \rho_{b8}+1-\ \sqrt{\left(2\ \times\ \rho_{b8}+1\right)^{2}-8\ \times\ \left(\rho_{b8}-\ \rho_{b4}\right)}\right]}{2}\)
where ρb8 and ρb4 are the at-the-ground reflectance of bands 8 and 4, respectively. Then, for each date, data for all pixels corresponding to field plots were extracted, including between 5 and 7 pixels for each plot (the central coordinates of each plot are shown in Table S1). All the image processing and data extraction for Sentinel-2 imaginery was done in the SNAP software (ESA, European Union). The Sen2Cor plugin was used for the TOA to BOA correction, and the MSAVI2 processor for index estimation.

Statistical analysis

To evaluate the effect of the treatments on each soil variable (total nitrogen, inorganic nitrogen, available phosphorus, carbon content, C/N ratio and pH), ANOVAs using linear models (LM) were used. ANOVAs were also tested to assess the effect of the treatments on the percent cover of each species and leaf nitrogen content. To evaluate the effect of the interaction between treatments and years on the number of dead plants, a generalized linear model (GLM) with binomial distribution and test-χ2 was tested. Differences in tussock size between treatments were evaluated with an ANOVA for each species.P. humilis was excluded from this analysis due to the lack of data in the C and +W. Generalized least squares models (GLSs) were used to evaluate the effect of the interaction between treatments and years on the ANPP of each shrub species. In this model, a temporal correlation structure between years was used for each sample (corAR1). Then, the shrub ANPP, integrating the ANPP of the three species of shrubs, the grass ANPP and the total ANPP (shrubs + grasses) were evaluated with ANOVAs, using LMs, to test the effect of the treatments in the different years. Multiple regression analysis, using LMs, were performed to evaluate the effect of the treatments on the relationship between annual precipitation, or precipitation plus irrigation in the case of +W and +NW treatments, and ANPP (shrub ANPP, grass ANPP and total ANPP). Generalized additive models (GAMs) were used to adjust the vegetation indices (NDVI and MSAVI2) through year for each treatment. In addition, an ANOVA with LM test was used to test the effect of the treatments on daily NDVI.
All statistical analyses were performed with R software version 4.1.2 (R Development Core Team, 2023). The GLSs were carried out using the “gls” function of the R package “nlme” version 3.1-162 (Pinheiro et al., 2023; Pinheiro & Bates, 2000). GLMs were performed using the “glm” function of the R package ”lme4” version 1.1-15 (Bates et al., 2015). GAMs were carried out using the “gam” function of the R package “mgcv” version 1.8-31 (Wood, 2004). When necessary, all models were adjusted using variance models. Model selection was based on the Akaike’s information criterion (corrected for small sample size, AICc) (Burnham & Anderson, 2002). Tukey’s post-hoc analysis were used for multiple comparisons of LMs and GLSs models when the F-test or χ2 tests were significant, using the “glht” function of the R package “multcomp” version 1.4-8 (Hothorn et al., 2008). The multiple comparisons of GAMs models were performed with the “emmeans” function of the R package “emmeans” version 1.8.6 (Lenth, 2023). The “visreg” function of the R package “visreg” version 2.7.0 (Breheny & Burchett, 2017) and the “ggplot” function of the R package “ggplot2” version 3.3.5 (Wickham, 2016) were used to graphically represent the results.
Results