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).