Livestock grazing and aridity reduce the functional diversity of biocrusts

Livestock grazing and climate change are two of the most important global change drivers affecting ecosystem functioning in drylands. Grazing and climate are known to influence the cover and composition of biocrusts, which are substantial components of dryland soils globally. Much less is known, however, about how these global change drivers affect the functional diversity of biocrust communities in these ecosystems. Here, we evaluate the role of increasing aridity and grazing intensity in driving the functional diversity of biocrusts. We collected data on multiple biocrust functional traits and community composition, recent and historic grazing intensity, and vascular plants at 151 sites from drylands in eastern Australia. We then used structural equation modelling and a fourth corner analysis to examine the combined effects of aridity and grazing on biocrust functional diversity and individual functional traits. Aridity had a significant direct suppressive effect on biocrust functional diversity. Effects of grazing by livestock, kangaroos and rabbits on functional diversity were predominantly indirect and suppressive, mediated by a reduction in biocrust cover. Grazing did, however, promote functional diversity via an increase in vascular plant richness, with a concomitant increase in biocrust richness. The overall effect of grazing on biocrust functional diversity however was negative. Fourth corner analyses revealed that livestock grazing had a significant negative effect on the ability of biocrusts to stabilise the soil. Aridity had strong negative effects on biocrust height and their ability to absorb water and capture sediment. Few significant relationships were detected between enzyme-related traits and environmental variables. Our findings provide novel evidence that the combination of increasing aridity and intensified livestock grazing will reduce the functional diversity and capabilities of biocrust communities, with resultant declines in ecosystem functioning.


Introduction
Biocrusts inhabit most terrestrial ecosystems but are particularly abundant in drylands as a consequence of their high tolerance to temperature and moisture extremes (Weber et al. 2016). In these ecosystems, biocrusts play important roles in moderating key ecosystem functions such as water infiltration and nutrient cycling (Bowker et al. 2013;Delgado-Baquerizo et al. 2013). Biocrust attributes such as cover and taxonomic richness have been reported to be highly vulnerable to climate change and land use intensification (Ferrenberg et al. 2015;Concostrina-Zubiri et al. 2016). Much less is known, however, about how the functional diversity of biocrusts will respond to the combined effects of aridity and land use intensification.
The way in which biodiversity, commonly quantified as taxonomic diversity (e.g. number of species), interacts with ecosystem functioning has been long studied. In these studies, however, taxonomic diversity is a surrogate for functional diversity, relying on a correlation that is frequently, but not consistently, positive (Díaz and Cabido 2001). For example, studies that measure both functional and taxonomic diversity often show that both respond differently to the same environmental drivers (de Bello et al. 2006;Ernst et al. 2006;Villéger et al. 2010). Since functional effect traits directly measure the mean effects of taxa or functional types on ecosystem functioning, conclusions drawn from the link between functional diversity indices and ecosystem functioning are theoretically more valid than those drawn from simple taxonomic diversity.
Livestock grazing and climate change are two of the most important drivers of ecosystem function in drylands . Overgrazing by livestock alters ecosystem functions, particularly in drier environments (Eldridge et al. 2016b), and pressures from grazing are predicted to intensify by up to 70% by 2050 (FAO 2017). Livestock grazing reduces crust cover directly by trampling (Eldridge 1998) and indirectly, via its effects on vascular plant cover (Tabeni et al. 2014). Climate is also likely to have substantial effects on biocrusts by altering composition Ochoa-Hueso et al. 2016) and the ability of crusts to produce radiation-resistant pigments (Belnap et al. 2004;Reed et al. 2012;Ferrenberg et al. 2015). Improving our understanding of how grazing and climate affect biocrust function is critically important if we are to manage drylands sustainably.
Despite the mounting literature on the effects of livestock disturbance and climate change on biocrust attributes, there are three major knowledge gaps. First, most studies have focused on biocrust cover and taxonomic richness, but little is known of the combined effects of climate change and disturbance in driving functional diversity of biocrusts. Second, most studies provide evidence of the vulnerability of biocrusts to disturbances at the local scale, but much less is known about how climate and grazing affect biocrust functional diversity at regional scales, which are aligned with the scales at which aridity and grazing operate ). Finally, we lack studies exploring the role of long-term climate changes (e.g. aridity) and land use intensification across environmental gradients (Breal world^situations), as most information comes from laboratory-based studies simulating short-term climate changes (e.g. drought or sudden increases in temperature).
A recent study showed a positive relationship between rainfall and the functional diversity of biocrusts on the Iberian Peninsula (Concostrina-Zubiri et al. 2014). Others have provided evidence of multiple interactions among biocrust traits and environmental variables (Bowker et al. 2011;Michel et al. 2013), or similarities in traits among functional groups of biocrusts (Gavazov et al. 2010;Kidron and Tal 2012;Mallen-Cooper and Eldridge 2016). Here, we aim to build a system-level understanding of the direct roles of aridity and grazing in driving biocrust function, and indirectly, via multiple environmental drivers and biocrust attributes. To do this we used two different modelling approaches; a structural equation model, to examine the combined effects of aridity and grazing on biocrust functional diversity, and a fourth corner model, to investigate the effects of these drivers on specific traits of species (Brown et al. 2014).
With limited previous research on biocrust function, it is difficult to make a comprehensive set of a priori predictions about environment-trait interactions derived from the fourth corner model. Nevertheless, Concostrina-Zubiri et al. (2014) showed that drier regions support biocrust taxa with thicker attachment structures. Consequently, we expected that our drylands would support taxa with longer roots (or rhizines in lichens). We also expected grazing to hinder the establishment of taller taxa, as they are more susceptible to trampling, but that aridity would promote communities with taller taxa, as they are better able to withstand the high levels of erosion and regular depositional events associated with more arid regions (Read et al. 2014).

Study area
Our study area covered an area of about 0.4 M ha in eastern Australia (31.5°S -36°S, 145°E -148.3°E; Fig.  1). The climate in this region is mostly semi-arid and drysubhumid (Aridity Index = 0.19 to 0.63; FAO 2013), and rain falls uniformly throughout the year, with slightly more rainfall in winter in the southernmost sites. Average annual rainfall varied along a northwesterly gradient from 850 mm in the southeast to 320 mm in the northwest. Average annual temperatures (~18°C) varied little across the area (BOM 2017). The 151 sites were non-contiguous fragments of woodland dominated by the tree species White cypress pine (Callitris glaucophylla). These communities have been used extensively for livestock grazing and forestry over the past~200 years (Thompson and Eldridge 2005). The soils are relatively deep sands to clay loams that support a moderate cover of biocrusts (Fig. 1).

Data collection
In order to capture a range of grazing intensities, we selected sites from a range of tenures including roadside reserves used for moving livestock (Travelling Stock Reserves), conservation reserves, exclosures and private land. We also used Arc GIS to stratify sites to capture a range of distances from permanent water, which is regarded as a useful surrogate of grazing intensity (Andrew 1988).
At each site we established a 200 m transect that ran perpendicular to the nearest livestock watering point. Along this transect we placed five 25 m 2 (5 m × 5 m) quadrats (henceforth 'large quadrat') 50 m apart (i.e. at the 0 m, 50 m, 100 m, 150 m and 200 m points). Within each large quadrat, we selected a uniform area of biocrust for placement of a smaller 0.25 m 2 quadrat ('small quadrat'). Within each small quadrat we estimated total crust cover. At the 50 m position at all sites we also counted the number of individuals of each biocrust morphospecies within the small quadrat using a 20× hand lens, and collected representative samples of all taxa present. These were placed into a paper bag and transported back to the laboratory, where they were identified using dissecting and compound microscopes and standard taxonomic keys described in Mallen-Cooper and . To assess abundance, we counted the number of individuals in a sample. However, we recognise that the assessment of biocrust abundance can be problematic because counting the discrete crust entities (the method we used) equates a nearmicroscopic moss with a crustose lichen the size of a small dinner plate. This could potentially skew the A B D C E Fig. 1 Images of short (a) and tall (b) mosses, and crustose (c) and squamulose (d) lichens from semi-arid Australia, and site map (e) showing the location of the study sites in eastern Australia. Photographs: Heino Lepp functional dispersion centroid towards small mosses and away from large lichens with continuous thalli. Biomass would be more representative than counts in this case, particularly as most traits were measured per gram.
We identified and estimated the cover of groundstorey plants and shrubs in the five large quadrats at each site (Eldridge et al. 2016a). To estimate recent grazing intensity, we first identified and counted the faecal pellets of all mammalian herbivores (cattle, sheep/goat, kangaroo, rabbit) in the large quadrat at each site. We were unable to differentiate between sheep and goat dung, or between European rabbit (Oryctolagus cuniculus) and European hare (Lepus europaeus) dung. Samples of dung were also collected at 20 sites, dried and weighed in order to calculate the relationship between dung counts and dry mass for each herbivore. We used this relationship to calculate the dry mass of dung at each site as a measure of recent grazing intensity (see Eldridge et al. 2016a). For each site we recorded the width and depth of all livestock tracks crossing the 200 m transect, calculated the cross-sectional area of tracks, and used this as a measure of historic grazing history. We extracted values of the Aridity Index (AI), which measures the relationship between precipitation and potential evapotranspiration, from the CGIAR-CSI Global-Aridity and Global-PET Database (http://www. cgiar-csi.org, Zomer et al. 2008). We then used this index to calculate aridity, as 1 -AI, so that increasing values of aridity corresponded with increasing dryness.

Biocrust trait data
Biocrust organisms have highly diverse morphologies, yet biocrust traits must be applicable to all constituent organisms, or at the very least, the macroscopic organisms. In a previous study we devised and measured eight functional effect traits in individual macroscopic biocrust species (Mallen-Cooper and Eldridge 2016): (1) absorptivity, (2) sediment capture, (3) height and (4) root (or rhizine) length relate to erosional and hydrological functions (Table 1; Bowker et al. 2010); (5) p h o s p h a t a s e ( P H O S ) a n d ( 6 ) N -a c e t y l -βglucosaminidase (NAG) activity provide a measure of the degree to which each biocrust species contributes to phosphorus and nitrogen cycles respectively; (7) β-Dcellobiosidase (CB) and (8)  The biocrusts of south eastern Australia are known to contain many rare species (Eldridge et al. 2006). Unfortunately, measuring traits across large environmental gradients is only feasible in the more abundant species. Rare species, however, are unlikely to contribute greatly to ecosystem functioning, according to the mass ratio hypothesis (Grime 1998).

Data analysis
We used two complementary approaches to analyse the data: a structural equation model (SEM; Shipley 2000) and a fourth corner model (Brown et al. 2014). To date, there are no standardised methods of measuring functional diversity (Pavoine and Bonsall 2011). The Laliberté and Legendre (2010) measure of functional diversity (Functional dispersion -FDis), however, has garnered considerable support in recent field and simulation studies due to several desirable properties (Paquette and Messier 2011;Mason et al. 2013;Gagic et al. 2015). First, the index can be weighted by species abundances, which has been shown to improve explanatory power (Gagic et al. 2015). Second, it varies independently of taxonomic richness and evenness, a necessary condition identified by Mason et al. (2005). Third, it tolerates any number and type of traits. Functional dispersion, however, does not allow for intraspecific variation, which accounts for a sizeable fraction of the total trait variation in real world situations (Siefert et al. 2015).
For each site, the functional diversity index (FDis) used here calculates the mean distance of species in multivariate trait space to a centroid that is weighted by species' abundances. To compute functional diversity, we firstly generated a Gower dissimilarity matrix between species for which we had sufficient trait data (Supplementary Material S1). Gower dissimilarity was used because it tolerates missing values. Principle coordinates analysis (PCoA) was performed on this matrix in the FD package version 1.0-12 in R (Laliberté et al. 2014), which derived a value of functional diversity for each of the 129 sites that contained biocrust samples. A site where only one species is present is allocated a value of zero for functional diversity. The 22 sites where no species were recorded were excluded from these analyses because the aim of this study was to determine how the traits of biocrusts vary, rather than the environmental conditions controlling the presence or absence of biocrusts.
We then used structural equation modelling to build a systems understanding of the direct and indirect effects of grazing and aridity on functional diversity. In our a priori model both aridity and grazing had direct effects on functional diversity and indirect effects, mediated by changes in perennial plant cover, plant richness, biocrust cover and biocrust richness. We used perennial plant cover rather than total cover because perennial plants are persistent and provide substantial ecosystem stability, and would be expected to support communities with an extensive crust cover. The effects of recent and historic grazing were combined in our models into a single composite variable ('grazing'), which increased with increasing grazing intensity. The SEM allowed us to partition direct and indirect effects of one variable upon another and to estimate the strengths of these multiple effects. To improve normality, the six predictor values were standardized (z-transformed) prior to analyses. Overall goodness of fit probability tests were performed to determine the absolute fit of the best models. The goodness of fit test estimates the longterm probability of the observed data given 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 was conducted using AMOS Software Version 20. The stability of these models was evaluated as described in Reisner et al. (2013).
Fourth corner analysis is being used increasingly in ecological studies. Unlike other methods, which are more indirect (such as community-weighted means), fourth corner analysis combines traits, environmental variables and species abundances into the same predictive model (Brown et al. 2014). It also provides a detailed insight into the effects of drivers like aridity and grazing on specific ecosystem functions. Fourth corner analysis describes the process of using three matricessites by environmental variables, sites by species' abundances, and species by functional traitsto derive a 'fourth' matrix of functional traits by environmental variables (Legendre et al. 1997). We followed the method of Brown et al. (2014), whereby a predictive species distribution model is fitted with environmental variables, species traits and their interaction as explanatory variables. The environment-trait interaction produces the fourth corner coefficients, which represent the pairwise relationships between each functional trait and environmental variable. To calculate the fourth corner coefficients, we used the mvabund package version 3.11.9 in R (Wang et al. 2012). The same predictors used in the structural equation model were used in this model, excepting biocrust cover and richness. Only species for which we had complete trait data were included in the fourth corner analyses because the calculations do not allow for missing values. Using the enzyme traits we had data for 104 sites and 16 species, and for physical traits 115 sites and 24 species. We therefore analysed enzyme and physical trait data sets separately in order not to compromise the power of the four physical traits. Both models were fitted with a Table 1 Summary of biocrust traits, including the number of species measured, the number of replicates per species and a brief description of the functional role(s) of each trait (PHOS = phosphatase activity, NAG = N-acetyl-β-glucosaminidase activity, CB = β-D-cellobiosidase activity, BG = β-glucosidase activity) The functional roles of each of these traits are also briefly described PHOS phosphatase activity, NAG N-acetyl-β-glucosaminidase activity, CB β-D-cellobiosidase activity and BG β-glucosidase activity. a 'Root' includes moss roots and lichen rhizines binomial distribution. A least absolute shrinkage and selection operator (LASSO) was used to determine which environment-trait interactions lowered the Bayesian Information Criterion (BIC) of the model. Interactions that did not lower the BIC were set to zero and were deemed non-significant. Livestock tracks and all dung variables were log-transformed in order to satisfy statistical assumptions. Graphical outputs were produced in the R packages ggplot2 version 2.2.0 (Wickham 2009) and lattice version 0.20-34 (Sarkar 2008). For the functional diversity index used in the SEMs, abundance was measured on a log (10) scale. For the fourth corner analysis, abundance was binary transformed (i.e. presence or absence) in order to improve model fit.

Results
We identified a total of 59 biocrust species across 151 sites. However, due to a shortage of trait data, only 40 of these species from 129 sites could be included in the structural equation model (Supplementary Material S1). Functional diversity was unrelated to biocrust species richness (Spearman's ρ = 0.006; P < 0.01) after removal of a cluster of zeroes; Fig. 2).
Grazing and aridity effects on functional diversity Our structural equation model explained 32% of the variance in biocrust functional diversity (χ 2 = 0.074, df = 1, P = 0.79; Fig. 3). We detected no direct effects of grazing on functional diversity, but indirect effects on biocrust richness mediated by changes in biocrust cover and plant richness. Perennial plant cover had a weak positive effect on biocrust functional diversity and aridity a direct suppressive effect (Fig. 3). Overall, the standardized total effects were strongly negative for aridity (−0.22) and strongly positive for biocrust richness (0.58) and cover (0.22), and to a lesser extent, plant richness and cover (0.12; Fig. 3).

Grazing and aridity interactions with species traits
There were a number of significant interactions among biocrust traits and environmental variables (Fig. 4, Supplementary Material S2). Aridity was negatively correlated with sediment capture, height and absorptivity of water. Thus, as aridity increased, biocrust communities not only became functionally less complex, but also shorter, less able to capture sediment and less able to absorb water. Grazing herbivores had predominantly negative associations with traits, with the exception of Fig. 2 Relationship between the functional diversity (FDis) and species richness of biocrusts (n sites = 129, n traits = 8, n species = 40). The point at (1,0) represents 21 sites which contained 1 species and were therefore assigned a value of 0 functional diversity. When these zeroes are removed, Spearman's ρ = 0.00635 rabbits, which enhanced the height of biocrust communities and their ability to capture sediment.
Few significant interactions were detected between enzyme traits and environmental variables, which may be a consequence of relatively low statistical power.

Discussion
Our study provides evidence that the functional diversity of biocrusts is highly vulnerable to increases in aridity and grazing. All effects from grazing on functional diversity were indirectly driven via changes in other biocrust attributes such as cover and taxonomic richness. Our results support the general notion that livestock grazing has overall negative impacts on ecosystem functions, though grazing by livestock, rabbits and native herbivores (kangaroos) showed a positive indirect effect via vascular plant richness, which did not counteract other indirect negative effects from grazing on biocrust functional diversity. The increase in plant richness is likely to be skewed towards exotic species, which tend to be enhanced by introduced grazers (Oduor et al. 2010). Aridity, however, showed strong and significant negative direct effects on the functional diversity of biocrusts, likely because of the reduced number of functional groups able to survive very arid conditions. Aridity also had negative effects on traits associated with erosion resistance, invertebrate habitat and hydrology. Our results imply that the combination of increasing aridity with intensified livestock grazing could have detrimental consequences for the functional composition of biocrust communities; with implications for the adequate functioning of drylands, the largest biome on Earth.
Indirect effect of grazing on functional diversity Interestingly, we did not observe a direct effect of grazing animals on biocrust functional diversity. Rather, any effect from grazing on functional diversity was indirectly driven via biocrust cover and richness. Livestock are well-known to directly reduce cover by trampling (e.g. Eldridge 1998), though other indirect mechanisms, which were not measured in the present study, are possible and would contribute to the direct effect in our simplified model. Our study suggests that grazing does not directly affect the suite of functional traits in a given biocrust community, but rather reduces the overall cover of these organisms, indirectly reducing the ranges of functional traits. In drylands, cover, or any measure of abundance, and species richness are almost always positively correlated, although a theoretical consensus explaining this relationship has yet to be reached (Bock et al. 2007). The accompanying increase in functional diversity that follows an increase in species richness is likely due to the 'niche complementarity effect', whereby a species-rich community is more likely to contain a greater range of functional traits, and thus use resources more efficiently, than a relatively species-poor community (Loreau 1998). Another more controversial explanation is the 'selection effect' (Díaz and Cabido 2001). We found a weak relationship between livestock tracks and 'root' length. These results are consistent with results from Liu et al. (2009) who found that heavy grazing reduced the ability of biocrusts to stabilise the soil. Increasing grazing by sheep reduced the capacity of biocrusts to capture sediments but had no effect on biocrust height, which is one property that influences sediment capture (Eldridge and Rosentreter 1999). This finding suggests that sheep are affecting a different property of biocrusts that influences their ability to capture sediment, perhaps by favouring species that grow in a dispersed, scattered arrangement rather than those that grow in dense clumps (e.g. Eldridge 1998). Unlike Read et al. (2014), who found that shorter biocrust taxa were indicators of heavily grazed sites while taller taxa were indicators of protected sites, we found no relationships between the intensity of livestock grazing and biocrust height. Interestingly, we found no relationship between absorptivity and grazing. This result contrasts with Concostrina-Zubiri et al. (2016), who found grazed plots had a lower abundance of fruticose lichens, which had particularly high water absorptivity and retention values.
Strong suppressive effects of aridity on functional diversity Aridity had a strong negative effect on the functional diversity of biocrusts and reduced the prevalence of biocrust taxa with extreme values of some traits. The extremes of temperature and moisture which characterise arid regions are important environmental filters, as not all biocrust functional groups are able to establish and survive in such harsh conditions (Coe et al. 2014). This process of habitat filtering constrains the range of various traits, resulting in lower values of functional diversity (Cornwell and Ackerly 2009). Increases in aridity have been reported to increase the cover of biocrusts ), but our results suggest that the likely increase in functional diversity linked to increases in cover with aridity will not offset the strong reduction in functional diversity related to aridity per se. This interesting result suggests that increases in aridity from climate change may promote the cover of biocrust communities with strongly reduced functional diversity. Ecosystem functioning in drylands is dependent to a large extent on biocrusts because as aridity increases and biocrust communities are less able to stabilise soils and cycle nutrients, dryland landscapes may become severely degraded.
We also found that increasing aridity strongly reduced sediment capture ability and generally reduced biocrust height. The lack of relationship between root length and aridity contravenes our hypothesis, revealing differential effects of aridity on the length and thickness of biocrust attachment structures. Although there are occasional resource pulses, arid regions generally have sparser vegetation, drier soils and stronger desiccating winds than more humid regions, leading to high rates of erosion . Biocrust communities in arid regions are thus likely to be assembled from taxa with a poor ability to capture sediment because those taxa which more readily capture sediments in their leaves and thalli are more likely to be buried. The negative relationship between aridity and height goes against our hypothesis. Read et al. (2014) found that short biocrust taxa were more likely to be found in areas protected from high winds, and attributed this to the high likelihood of burial by wind-borne sediment, thus preventing the establishment and persistence of these taxa in more exposed areas. Despite taller taxa being better equipped to deal with high rates of aeolian erosion, it is possible that there is a limitation on growth in arid conditions (e.g. productivity, soil nutrients, or desiccation). Since there is a reasonably strong positive correlation between height and sediment capture (Mallen-Cooper and Eldridge 2016), another possibility is that being tall increases a biocrust organism's likelihood of surviving burial but also increases the likelihood of being buried.

Concluding remarks
It is likely that the projected increase in aridity with climate change will be accompanied by a reduction in the functional diversity of biocrusts and in specific biocrust functional traits, irrespective of the reported increase in biocrust cover with aridity. Given that trait-based functional diversity is fundamentally linked to the provision of ecosystem goods and services (Díaz et al. 2007), a number of ecosystem services provided by biocrusts, such as soil stabilisation and nitrogen fixation, may be threatened under current climate change scenarios. It may be possible to artificially improve the resilience of biocrusts to climate change through 'assisted migration' (Young et al. 2016). This might necessitate introducing taxa adapted to very arid areas into areas that are expected to develop similar climatic profiles in time. However, such promising developments are still in their infancy.