1.    Biocrusts are critical components of drylands where they regulate a wide range of ecosystem functions, however, their response to the worldwide phenomenon of shrub encroachment and to livestock grazing, the most extensive land use in drylands, are not well studied. Grazing by livestock and increases in shrub cover could influence biocrust communities directly via trampling or shading, or indirectly, by altering biotic interactions amongst biocrust taxa. The extent of these changes in biocrust cover, diversity and composition are poorly known.
2.      We used linear models and structural equation modelling to examine the direct effects of grazing and shrubs on biocrust community composition and the indirect effects mediated by changes in species interactions.
3.      Biocrust richness and cover increased with increasing shrub cover at the site level. This pattern occurred despite the negative response we found (lower cover and richness) under shrub patches vs open areas, which was consistent irrespective of grazing level. Functional diversity and evenness were similar between shrubs and open at low grazing intensity, but at high grazing functional diversity was greater in the open. Competition between biocrust species was an important driver of their community assembly irrespective of shrub cover, grazing intensity or patch type. Structural equation models showed that the effects of grazing and shrub cover on functional evenness, functional diversity and richness were controlled by biotic interactions within the shrub microsites. In the open, however, these effects were either direct or mediated by changes in cover.
4.      Biocrust cover, species richness and functional diversity increase with shrub cover at the site scale, despite the negative effects at the microsite level. We demonstrate here that drivers of community assembly differ markedly at small spatial scales. Though biocrust communities were directly driven by environmental filtering in the open, biotic interactions played a fundamental role in their assembly when growing beneath shrubs.

At each site, five shrub and five open (interspace) patches (hereafter ‘microsites’) were randomly selected, and within each of these microsites, a 5 cm-diameter circular subplot was sampled under the canopy or in the open (N = 320 in total). Each subplot was carefully removed from the ground, placed in a paper bag, and stored in a cool place until biocrust cover by species was measured. Once in the laboratory, each taxon was identified to species level using keys in Filson and Rogers (1979), McCarthy (1991), Catcheside (1980), Scott (1985) and Scott and Stone (1976) and more recent generic revisions. Nomenclature followed Buck and Vitt (2006) for mosses, McCarthy (2006) for liverworts, McCarthy (2015) for lichens, and where appropriate, more recent taxonomic revisions. From these data, species richness and cover were directly obtained. Changes in species composition were assessed using non-metric multidimensional scaling ordination (nMDS) and the Bray-Curtis distance measure on log-transformed abundance data with the PRIMER/PERMANOVA statistical package (PRIMER-E Ltd., Plymouth Marine Laboratory, UK).

Functional attributes were assigned to each of our biocrust taxa based on eight functional traits: root (rhizine) length, taxon height, sediment capture, absorptivity, and the activity of four enzymes associated with carbon (b-glucosidase, b-D-cellobiosidase), nitrogen (N-acetyl-b-glucosaminidase) and phosphorus (phosphatase) cycling. These functional traits are linked to competitive ability and response to the environment in biocrust taxa, but also to their effects on important ecosystem functions such as nutrient cycling, water infiltration and resistance to erosion (see Bowker et al., 2010; Mallen-Cooper & Eldridge, 2016). From this information, we calculated two complementary metrics of functional diversity: RaoQ (sum of pairwise distances in multidimensional trait space), and ii) functional evenness (a measure of the extent to which different trait strategies are evenly distributed in the community; see Pakeman, 2014 for a review). Since RaoQ reflects pairwise differences in multidimensional trait space, it should be a good measure of functional similarity between two given species. Thus, assuming small fitness differences and hierarchical competition, species should differ as much as possible (i.e., high RaoQ) to be able to coexist (limiting similarity).

Functional metrics were calculated using the FD package in R (Laliberté et al., 2014). Functional dispersion and richness were also considered, but removed in further analyses due to their strong correlations with RaoQ (ρ = 0.97 [functional dispersion]) and species richness (ρ = 0.75 [functional richness]). Functional trait information was available for 23 of the 53 species we found, and these represented, which represented 56.1 ± 3.08% (mean ± SE) of the total cover within the 320 subplots. All functional diversity metrics were abundance-weighted.

Changes in biocrusts interactions

To calculate the degree of competition intransitivity, the inverse Markov-chain approach of Ulrich et al. (2014) was used as implemented in the free software Transitivity.  This approach provides a metric for intransitivity (I), which ranges between 0 (total competition hierarchy) and 1 (total competition intransitivity). The metric is based upon the number of competition reversals, i.e., how many times the competitor dominant(s) lose (see Laird & Schamp, 2006; Ulrich et al., 2014 for details). With this approach the number of species to analyse cannot be exceed one less than the number of sites, therefore in our case, four species. Thus, we analysed the degree of intransitivity on the four dominant species (as we had five replicates of each microsite × site combination). These were three lichen (Collema coccophorum, Psora crystallifera and Psora decipiens) and one moss (Didymodon torquatus) species. Together with the levels of intransitivity, the software calculates a measure (Match) of how well the simulated Markov chains fit the observed abundances. This Match metric is an indicator of the importance of competition for the assembly of the studied species (values close to one indicate that competition is strongly driving the observed abundances). The Match metric also indicates how reliably the level of intransitivity has been estimated. Only in three out of the 64 microsite × site combinations was this metric below 0.6 (60% of the observed abundances explained by the simulated Markov chains). Therefore, it can be assumed that competition was a dominant interaction among the studied species and that the level of intransitivity was reliably estimated (see Ulrich et al., 2014; Soliveres et al., 2015 for further details on the methods).

To account for changes in composition, we applied this approach to two different datasets: one including the four dominant species across all sites (therefore removing changes in composition), and one including the four dominant species at each site (therefore accounting for effects of changes in composition + environment). The intransitivity metrics obtained using the overall dominant species or those most dominant at each site rendered similar results (ρ = 0.41; P = 0.017; N = 33), although in the former case, this metric could only be calculated in half of the 64 microsite × site combinations because these species were not present in all subplots. Thus, we only report results from the metrics calculated for the four dominant species at each site.

Overall, we analysed changes in biocrust communities by examining total cover, species richness and composition, functional composition, the importance of competition for community assembly, and the degree of competition intransitivity. For the taxonomic and cover measures, we used our full database. For the functional diversity metrics, we used a reduced database of the 23 taxa for which we had functional trait data. Finally, to assess the interactions among different biocrust species, we used the four dominant species.

Statistical analyses

Changes in biocrust cover, species richness, functional diversity (functional evenness and diversity [RaoQ]) and the level of intransitivity among the four dominant species in each subplot) were analysed by using linear models with site level shrub cover (0-50%), microsite (shrub vs open) and grazing level (low vs high), and their interactions as fixed factors using the ‘lme4’ package (Bates et al., 2015) within R statistical software (Version 3.4.4, R core Team 2018).

To further interpret compositional changes, the degree of association of biocrust species in relation to patch type and grazing was measured with Indicator Species Analysis in R (De Caceres, 2013) using a data matrix of 56 species and 64 combinations of microsite by site. Indicator values combine information on relative abundance and frequency of species, and the indicator value is maximal (IV=100%) when all individuals of a given species are restricted to a particular microsite (e.g. shrub), and all samples from that particular microsite contain an occurrence of that species. Data were randomized among the factors and a Monte Carlo randomization procedure performed with 1000 iterations in order to determine the statistical significance of the indicator values.

Finally, we used Structural Equation Modelling (SEM) to examine the direct and indirect effects (i.e. mediated by biocrust cover and the way species interact) of increasing grazing intensity and shrub cover, on biocrust richness, functional diversity, and composition. In order to evaluate possible interactions among the importance of each driver and their direct and indirect effects for a given microsite, we fitted separate models for open and shrub patches. Structural equation modelling tests the plausibility of a causal model, based on a priori information (Supporting Information Appendix S3: Fig. S1), in explaining the relationships among different variables. In our model we predicted that grazing would have direct effects on individual crust attributes (e.g. functional composition, richness; Eldridge et al. 2013; Bowker et al. 2013; Concostrina-Zubiri et al., 2017) but also indirect effects via changes in intransitivity (Bowker et al., 2010; Soliveres et al., 2015) and crust cover (as larger patches are more likely to support more species which are likely more competitive, therefore modulating biotic interactions). Similarly, the changes in heterogeneity and productivity induced by shrub patches and increasing density at the plot scale should directly affect the composition of biocrust communities (Maestre et al., 2002), but also their cover and how they interact. Overall goodness of fit probability tests (χ² and Bollen-Stine) were performed to determine how well our a priori model structure fit the data. The goodness of fit test estimates the match between the variance/covariance matrix of the observed data and that expected under the a priori model structure. Thus, high probability values indicate that these models are highly plausible causal structures underlying the observed correlations. All SEM analysis were conducted using AMOS Software Version 20.