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