Data analysis
We performed all data analyses using the statistical program R (R Core
Team 2016). We checked for normality using the Shapiro-Wilk test and
made square root transformations of percentage of cover for bare soil,
darkly-pigmented cyanobacterial cover, lichen cover, and moss cover to
meet normality assumptions. Data on lightly-pigmented cyanobacterial and
late successional cover (the sum of darkly-pigmented cyanobacterial
cover, lichen cover, and mosses cover) did not need transformations. To
evaluate changes in community successional maturity, we compared the
percentage of cover of the biocrust functional groups between control
and treatment using a blocked ANOVA with site as a random effect.
Using the community ecology vegan package in R (Oksanen et
al . 2018), we summarized biocrust community composition with a
two-dimensional nonmetric multidimensional scaling (NMDS) ordination
using Bray-Curtis distance. This method allowed us to visualize and test
the direction and magnitude of community composition changes from
control to treatment; we extracted the NMDS coordinates and calculated
the difference (treatment minus control) along each axis. We performed a
one sample t-test on the values for each axis, the p-values were then
combined in a single Fisher’s C statistic (Shipley 2000), which follows
a chi-squared distribution with 2k degrees of freedom.
To further test for dissimilarities in the communities between the
control and drought treatments, we used a dissimilarity index test,
called a temporal beta-diversity index (TBI), treating the observations
between control and treatment as if they were surveys through time. TBI
is calculated as a common pairwise beta diversity index but allows a
test of significance at each site, which differs from most beta
diversity analyses (Legendre 2019). The dissimilarity index used for
community composition for each of the 25 sites was the percentage of
difference (Odum 1950), also known as the Bray-Curtis index and
equivalent to the Sorensen index, but for species abundances, which tend
to be more informative than incidence data (Legendre 2012). Further, the
dissimilarity index was partitioned into finer indices of losses (B) and
gains (C) of abundances-per-species. The mean of the differences between
the B and C statistics is also computed across all sites in a study. The
gains versus losses (C – B) difference across all sites were tested for
significance using a paired t-test to observe specific patterns obscured
in the overall test. All calculations were implemented in the TBI()
function, available in the R package adespatial (Dray et
al . 2020). We used the dissimilarity index and its components of losses
(B) and gains (C) as indices of community resistance to climate
disturbance, where a greater value indicates a greater degree of change
in the community. Multiple linear regressions were used to evaluate how
much variation of the community resistance was explained by the control
Shannon alpha diversity, aridity, and its interaction. To avoid
violations of assumptions, the normality of model residuals was verified
by the QQ-plot, model linearity was visually evaluated by inspecting the
residuals vs fitted values plot, and we identified outliers using the
outlierTest function from the car package. Finally, plots were
made using the ggplot2 package (Wickham 2016).