Contrasting Effects of Aridity and Grazing Intensity on Multiple Ecosystem Functions and Services in Australian Woodlands

Global change is expected to reduce the provision of multiple ecosystem services in drylands, the largest biome on Earth. Understanding the relative importance of climate change and overgrazing on ecosystems services is critical for predicting the effects of global change on ecosystem well‐being. We generated a system‐level understanding of the effects of climate (aridity intensity) and land use intensification (herbivore grazing intensity) on four regulating ecosystem services (C‐storage, N‐availability and P‐availability, and organic matter decomposition) and one provisioning service (plant production) in wooded drylands from eastern Australia. Climate change and grazing intensity had different effects on multiple ecosystem services. Increasing aridity from 0·19 (dry subhumid) to 0·63 (arid) had consistent suppressive effects on C‐storage, N‐availability, decomposition and plant biomass services, but not on P‐availability. The magnitude of these suppressive effects was greater than any effects due to grazing. All sites showed evidence of kangaroo grazing, but the heaviest grazing was due to cattle (dung: range 0–4545 kg ha−1; mean 142 kg ha−1). Any effects of grazing on ecosystem services were herbivore specific and ranged from positive to neutral or negative. Sheep, and to a lesser extent cattle, were associated with greater N‐availability. Rabbits, however, had a greater effect on P‐availability than aridity. Our study suggests that increases in livestock grazing may fail to sustain ecosystem services because of the generally stronger negative effect of increasing aridity on most ecosystem services in our model dryland. These services are likely therefore to decline with global increases in aridity. Copyright © 2017 John Wiley & Sons, Ltd.


INTRODUCTION
One of the terrestrial systems most at risk of environmental degradation is drylands, which occupy about 45% of Earth's land mass (Prăvălie, 2016) and support about 38% of its people (Reynolds et al., 2007). Drylands occur disproportionately in developing countries (Prăvălie, 2016), are agriculturally marginal and support many socially disadvantaged groups that rely heavily on grazing and farming for their livelihoods (Powell et al., 2004). Two of the most important drivers of ecosystem function in drylands are overgrazing and increased aridity (Maestre et al., 2016). Understanding the roles of land use intensification (overgrazing) and climate change (increasing aridity) is of paramount importance if we are up to predict how ecosystem services might change in the face of a well-being.
Overgrazing by livestock is emerging as one of the greatest threats to global biodiversity conservation and ecosystem sustainability (Fleischner, 1994;Vázquez & Simberloff, 2003). Over the past three centuries, the total global area of land dedicated to grazing has increased more than sixfold, from 524 to 3451 million hectares, largely as a result of population growth and colonisation (Goldewijk, 2001(Goldewijk, , 2005. Consequently, almost a quarter of the global terrestrial land mass is now under grazing by domestic livestock (Asner et al., 2004). About 60% of land dedicated to grazing is considered drylands, where evapotranspiration exceeds rainfall (Maestre et al., 2012;Prăvălie, 2016) and where the effects of human-induced changes such as overgrazing are likely to have the largest effect on the provision of ecosystem services . Carbon storage and N-availability in plant biomass and forage supply may be greatest in moderately grazed drylands (Oñatibia et al., 2015). However, overgrazing reduces vegetation cover and biomass (Jones, 2000), leading to a decline in the spatial heterogeneity of litter cover (Daryanto et al., 2013b). Reduced litter cover can lead to reduced decomposition and therefore smaller pools of soil C and N.
The effects of overgrazing on ecosystem functions and services may differ, however, across different climatic conditions. For example, Rabbi et al. (2015) found that land management did not influence the capacity of Australian drylands to soil carbon storage, a key ecosystem service, owing to the strong limitations in water availability, the main limiting factor in these ecosystems. Because of this, the effect of overgrazing cannot be evaluated in isolation and needs to consider the relative effects of multiple environmental factors, which are likely to influence how ecosystem services respond to changing land use intensity. This has not been widely reported in current literature.
The second driver, aridity, has negative effects on a large number of ecosystem functions and services via multiple pathways. Increased aridity reduces the diversity and abundance of soil microbial communities that carry out multiple soil functions (Maestre et al., 2012Delgado-Baquerizo et al., 2016b), reduces C-storage and Navailability by suppressing plant production and therefore nutrient inputs into the soil and alters the concentration of some soil enzymes, which control nutrient production processes (Li & Sarah, 2003;Delgado-Baquerizo et al., 2013a. However, the suppressive effect of aridity on ecosystem functions and services is highly dependent on the specific service. Increasing aridity, for example, has been shown to enhance soil P, largely by increasing the exposure of P-rich parent material and enhancing the amount of P bound to soil carbonates, which are abundant under the most arid conditions (Delgado-Baquerizo et al., 2013b). This may also have feedback effects on plant community composition by promoting plants that can grow in soil with lower N:P ratios (Güsewell & Bollens, 2003), thereby reducing ecosystem resilience. Single functions, therefore, should not be seen in isolation because overall ecosystem functioning will depend on a range of functions or services operating together (multifunctionality sensu Gamfeldt et al., 2008). In drylands, resilience of ecosystems to ongoing global environmental change such as increased aridity and/or drought is enhanced by increasing plant taxonomic and functional diversity and woody plant cover (Gaitan et al., 2014;Valencia et al., 2015), or by maintaining a diverse and abundant soil microbial community . The maintenance of high levels of diversity in drylands is critical in order to maximise the provision of a range of ecosystem services on which a range of biota, including humans, depend (Naidoo et al., 2008;Bellard et al., 2012).
We know relatively little about the relative strengths of the two main drivers, grazing and aridity, on ecosystem services in drylands (although see Ibáñez et al., 2014 for a modelling approach). This lack of knowledge hampers our capacity to manage change in drylands (e.g. Prăvălie, 2016) or to determine the interactive effects of continued grazing at current levels under a regime of reduced rainfall. We developed an a priori model, based upon the known relationships among the various predictor and response variables ( Figure 1a; Table I), and included in our model woody plant cover and soil health. Soil health was included because it is an essential component of human and environmental well-being (Pankhurst et al., 1997). It is also particularly important in drylands because the soils tend to be shallow and low in available nutrients and are often sparsely vegetated and susceptible to wind and water erosion. Both (e) (f) Figure 1. Structural equation models for (a) the a priori model and (b-f) five ecosystem services. Grazing is a composite variable comprising recent grazing by all herbivores, and historic grazing by livestock. Standardised path coefficients, embedded within the arrows, are analogous to partial correlation coefficients, and indicate the effect size of the relationship. Continuous and dashed arrows indicate positive and negative relationships, respectively. The width of arrows is proportional to the strength of path coefficients. The proportion of variance explained (R 2 ) is shown in each figure. Only significant pathways are shown in the models. Model fit: χ 2 = 8·48, d.f. = 9, p = 0·49. 1 p = 0·09; woody cover and soil health have been shown to mediate the effects of both aridity and grazing on the provision of ecosystem services and functions (e.g. Maestre et al., 2016). We predicted that increasing aridity (path 3) would have negative effects on ecosystem functions associated with C-storage, N-availability and plant production (Maestre et al., 2012(Maestre et al., , 2016Zhang et al., 2016;Delgado-Baquerizo et al., 2016a), but that P-availability would be advantaged by increased aridity due to the re-distribution of P-rich sediments resulting from increased degradation (Delgado-Baquerizo et al., 2013b; hypothesis 1). Second, we expected to detect substantial herbivore-specific effects on services and functions (path 4). Specifically, kangaroos (i.e. native species) were predicted to have few effects on functions and services because they have different grazing patterns from livestock and have coevolved with plants and soils (hypothesis 2). Conversely, we expected to detect substantial increases in soil P-availability with increased grazing by rabbits or livestock because rabbits create substantial soil disturbance, which exposes P-rich subsoil (hypothesis 3), and high levels of livestock grazing results in the concentration of P-rich dung (Duncan et al., 2008;hypothesis 4). Finally, we predicted that increased livestock grazing would be associated with substantial reductions in plant production (hypothesis 5).

Study Area
The study was undertaken in south-eastern Australia ( Figure 2). Landforms in the study area are predominantly plains of coarse colluvium, and low ridges and valleys and slopes to 3%. The soils are dominated by well-drained gradational red loams and red-brown duplex soils with neutral to slightly acidic surfaces (pH 5-7; Thompson & Eldridge, 2005). The climate across the area is Mediterranean and typically semi-arid (Table II). Mean annual rainfall varies from 385 to 460 mm and ranges from being evenly distributed throughout the year in the east and central parts of the study area to 20% more rainfall during the six cooler months in the south and south-west. Average temperature in the study area ranges from 9°C to 12°C for May to September and from 18°C to 24°C for October to April, and the mean annual temperature between 15°C and 18°C.

Site Selection, Woody Cover, Grazing Intensity and Aridity
Between September and November 2013, we surveyed 151 sites that were characterised by the community dominant Callitris glaucophylla Joy Thomps. & L. A. S. Johnson ( Figure 2). The study sites were first identified using Arc GIS and pre-inspected to ensure that they were more than 250 m from any road. In order to sample across a full spectrum of grazing intensities, we selected some low intensity and long ungrazed sites from conservation reserves, road verges with intermittent grazing, commercial forests, conservation reserves and long-term grazing exclosures as well as high-intensity grazed sites (Table II). At each site, we positioned a 200-m-long transect within which were placed five 25-m 2 (5 m × 5 m) plots (hereafter 'large quadrat') every 50 m (i.e. 0, 50, 100, 150 and 200 m). A smaller (0·5 m × 0·5 m) quadrat (hereafter, 'small quadrat') was located in the centre of the large quadrat. Along the line transect, the projected cover of trees (>3 m tall) and shrubs (0·5-3 m) was recorded at 100 points every 2 m to provide a measure of woody (tree + shrub) cover. We also measured the width and depth of all livestock tracks crossing the 200-m transect.
We used the total cross-sectional area, averaged to a 100-m transect (cm 2 /100 m) of stock tracks as a proxy of historic (long-term) grazing. In order to assess recent grazing intensity at the sites, we identified and counted dung/pellets in both of the small and Table I. Hypothesised mechanisms underlying the grazing a priori meta-model shown in Figure 1a Path Hypothesised mechanism 1 (À) Increasing aridity reduces woody cover in semi-arid landscapes (Westerband et al., 2015). 2 (+/À) Ecosystem multifunctionality peaks at 41-60% of relative woody covers (RWC), but aridity alters the effect of RWC on multifunctionality. The RWC-multifunctionality relationship was linear positive in dry-subhumid sites, but it shifted into hump shaped and ended a negative relationship in most arid conditions (Soliveres et al., 2014). 3 (À) Aridity has negative effects on ecosystem multifunctions or services (Maestre et al., 2012;Delgado-Baquerizo et al., 2016a), because it reduces soil microbial diversity and abundance that promote ecosystem multifunctions Delgado-Baquerizo et al., 2016b) and it reduces organic C, total N and enzyme activities such as phosphatase activity (Li & Sarah, 2003;Delgado-Baquerizo et al., 2013a; (+) but aridity increases total P (Delgado-Baquerizo et al., 2013b). 4 (À) Grazing has significant negative impact on overall ecosystem multifunction and multiservices (Zhang et al., 2016) by dampening positive effect of shrubs (Eldridge et al., 2013. It reduces plant biomass and litter cover (Eldridge et al., , 2016b, and therefore soil organic C (Daryanto et al., 2013a), and increases bare soil (Daryanto et al., 2012). (+/À) Under high intensity of grazing microbial biomass-C, phosphatase and dehydrogenase activities increase owing to herbivore urine and dung at intercanopy but reduces β-glucosidase activity, organic C and total N at under plant canopy soil (Prieto et al., 2011;Olivera et al., 2014). 5 (À) Soil organisms, one of the soil health indicators, are very sensitive to climate (Doran & Zeiss, 2000); negative correlation between aridity and plant cover may enable soil erosion (Delgado-Baquerizo et al., 2013b). 6 (À) High stocking rate increases soil compaction and bulk density (Pulido et al., 2016), and therefore reduces infiltration and macroporosity (Castellano & Valone, 2007;du Toit et al., 2009), as well as nutrients and stability indices (Eldridge et al., 2013). 7 (+) Healthy soils with higher plant, litter and crust cover will increase organic inputs into soil, which support diverse microbial community and enrich the soil nutrient pool, will promote overall ecosystem functions.
large plots separately by herbivore type. Kangaroo (Macropus spp.), rabbit (i.e. rabbits and hares; Oryctolagus cuniculus L. and Lepus europaeus Pallas) and sheep (which included sheep Ovis aries L. and goats Capra hircus L.) dung/pellets were counted in the small quadrats, and cattle dung (Bos taurus L.), sheep and kangaroo dung/pellets counted in the large quadrats. For cattle, we counted dung events rather than individual fragments, that is, we considered a number of small fragments to have originated from one dung event, if the fragments were within an area of a few metres. Dung and pellet samples of each type were collected, oven dried at 40°C for 24 h and weighed to estimate the mass of individual pellets, or in the case of cattle, dung events. Average mass of dung was then used to calculate the total mass of each type of dung per hectare. We used algorithms, developed previously for the study area (Eldridge et al., , 2016b, to calculate the total ovendried mass of dung per hectare per herbivore on the basis of the number of pellets recorded in the field. This total oven-dried mass of dung was used as our measure of recent grazing intensity for each herbivore (Eldridge et al., , 2016b. Where dung from the same herbivore was assessed in both the large and small quadrats, we derived an average mass per hectare on the basis of the large quadrat for that herbivore type. We used dung as a measure of recent grazing because dung persists for about 3 years in the field before it is decomposed. Thus, it is a useful proxy of short-term grazing by herbivores. Stock tracks provide a long-term integrated history of continued use by livestock. Aridity was calculated as 1 À AI, where AI = precipitation/potential evapotranspiration using Food and Agriculture Organization's global aridity map (http://ref. data.fao.org/map). Subtracting AI from 1 changes the direction of the index such that larger values are more arid.

Assessment of Soil Health and Plant Production (Biomass)
We assessed the status and morphology of the soil surface within the small quadrats using rigorous, field-based protocols (Tongway, 1995). Within the small quadrats, we measured 11 soil surface attributes: surface roughness, crust resistance, crust brokenness, crust stability, surface integrity (cover of uneroded surface), cover of deposited material (e.g. sand), biocrust cover, plant basal cover, litter cover, litter origin and the degree of litter incorporation (Table S1). From these data, we calculated a soil health index as the mathematical mean of the 11 attributes following standardisation (z-transformation).
To estimate plant production, we took oblique photographs of all small quadrats, then clipped, oven dried (45°C for 24 h) and weighed all above-ground material rooted within the small quadrats at the 50-m mark along the transect. Biomass values for all small quadrats were then estimated and calibrated against the 151 actual biomass values using the photostandards. Using this double-sampling approach, the predictive power, on the basis of a second-order polynomial, was R 2 = 0·69.

Soil Sampling and Analyses
During the field survey, we collected about 500 g of soil, from the surface 5 cm, from the centre of each small quadrat. A total of 755 soil samples were analysed in this study, from 151 sites. The soil samples were air dried and passed through a 2-mm sieve to remove roots, organic debris and stones prior to chemical analyses. Total C and N were assessed using high-intensity combustion (LECO CNS-2000; LECO Corporation, St. Joseph, MI, USA), available (Olsen) P according to Colwell (1963). Labile carbon was assessed by measuring the change in absorbance when slightly alkaline KMnO 4 reacts with the most readily oxidisable (active) forms of soil C to convert Mn (VII) to Mn (II; Weil et al., 2003). Ammonium and nitrate concentrations were measured using flow injection analysis (QuikChem® 8500; Lachat Instruments, Milwaukee, WI, USA) following the extraction with 0·5 M K 2 SO 4 . For the activity of the enzyme β-glucosidase, a mixture of 1 g of air-dried soil and 33 ml of sodium acetate buffer (pH < 7·5) was shaken at 200 rpm on an orbital shaker for 30 min, and 800 μl soil slurry was sampled and 200 μl substrate of 4methylumbelliferyl β-D-glucopyranoside solution was added to the slurry. The 1000-μl (1 ml) solution was incubated at 25°C for 3 h, and the activity (nmol activity g À1 dry soil À1 h À1 ) was measured at the 365-nm excitation wavelength and 450 nm of emission wavelength in a microplate reader. The same procedure was used, but with different substrate solutions, for an additional three enzymes. Thus 4-methylumbelliferyl β-D-cellobioside was used for cellobiosidase, 4-methylumbelliferyl N-β-D-glucosaminide for N-acetyl-β-glucosaminidase and 4-methylumbelliferyl phosphatase for phosphatase activity (Bell et al., 2013).

Ecosystem Services Scoring
In this study, we used the term 'ecosystem function' to define the 11 individual functions that we measured. The term 'ecosystem service' or service refers to an assemblage of functions grouped according to similar functions that they perform, for example, N-availability. The 11 separate functions were grouped into five services, as follows: C-storage (total C and labile C), N-availability (total N, ammonia and nitrate), P-availability (Colwell P only), decomposition (extracellular enzyme activities of β-glucosidase, cellobiosidase and β-N-acetylhexosaminidase and phosphatase) and plant biomass (plant biomass only). Decomposition, C-storage, and N-availability and P-availability are regulating services, and biomass is a provisioning service (MEA, 2005). Carbon storage and N-availability and Pavailability are proxies of ecosystem processes such as C sequestration and N and P mineralisation. In the case of available P, note that while total P is considered a soil property (mainly provided by bedrock), available P is the result of P mineralisation and solubilisation of soil organic matter and minerals in the bedrock, a process conducted by plant roots and microbial communities. The index of each service was calculated by averaging the standardised (z-transformed) value of each function for each site into a single metric (multifunctionality sensu Maestre et al., 2012;Soliveres et al., 2014;Delgado-Baquerizo et al., 2016aZhang et al., 2016). This allowed us to transform all services to a common scale of standard deviation units (Maestre et al., 2012). Integrating or averaging the standardised multiple functions into a single metric for overall ecosystem enables compensation of the decrease in one or several function by the increase of one or several functions (Gamfeldt et al., 2008;Quero et al., 2013); therefore, it is a highly informative quantitative measure of the overall ecosystem performance.

Statistical Analyses and Modelling
Structural equation modelling (SEM) was used to build a system-level understanding on the effects of grazing and aridity on the ecosystem services. An a priori model (Figure 1a; Table I) was applied to the five ecosystem service indices. In the model, grazing and aridity were included as the main factors on the ecosystem services, directly and indirectly via soil health and tree cover. In our models, the a Woody (shrub + tree) cover can exceed 100% where shrubs occur beneath trees. Aridity = 1 À Food and Agriculture Organization aridity index.
effects of recent and historic grazing were combined into a single composite variable ('grazing'). Increases in this composite variable corresponded to increasing grazing intensity. The use of composite variables does not alter the underlying SEMs but collapses the effects of multiple, conceptually related variables into a single combined effect, aiding the interpretation of model results (Grace, 2006). The SEM allowed us to partition direct and indirect effects of one variable upon another and to estimate the strengths of these multiple effects. We included woody cover, aridity, soil health and grazing in our models because we considered them to be the most informative predictors of ecosystem services out of a suite of many potential predictors. We are aware, however, that other variables not considered here could have been included, such as trampling by cattle or the number of pits dug by rabbits. These may have provided additional insights into the mechanisms underlying grazing effects on ecosystem services. To improve normality, we standardised (z-transformed) values for soil health, woody cover, aridity and grazing 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 long-term 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 (IBM Corporation, Armonk, New York, USA). The stability of these models was evaluated as described in Reisner et al. (2013).

RESULTS
Aridity (1 À aridity index) at our study sites ranged from 0·19 (dry subhumid) to 0·63 (arid), with a mean value of 0·36 (semi-arid; Table II). Average soil health index and woody cover values were 44·1% and 41·9%, respectively. Recent grazing intensity data, as measured by the mass of herbivore dung, varied markedly, from no grazing at some sites to more than 4·5 t ha À1 of dung at one site grazed by cattle. Kangaroos were the only herbivore to be recorded at all 151 sites. Plant biomass varied markedly (0·19 to 3·16 t ha À1 ) with a mean of 1·14 t ha À1 , typical of these semi-arid woodlands. Soil enzyme concentrations ranged from 63 (cellobiosidase) to 112·9 nmol g À1 h À1 (phosphatase). Total N (0·15%) and C (2·08%) reflected typically low values found for eastern Australian rangelands.

Effects of Grazing on Multiple Ecosystem Services
Grazing had moderate to strong direct and positive effects on all services (Figure 1b-f), but also strong, indirect negative effects on N-availability, C-storage and plant biomass, via a suppression of the positive effects of soil health on these three services (Figure 1b-c, e). Conversely, the indirect effect of grazing was positive on P-availability via reductions in soil health (Figure 1d). The standardised total effects (STEs) from the SEM, that is, the sum of all direct and indirect effects, indicated a general positive effect of grazing on N-availability, which was due mainly to sheep and cattle (Figure 3a), but an overall neutral effect on Cstorage ( Figure 3b). Grazing had an overall positive effect on P-availability, owing mainly to rabbit grazing (Figure 3  c). Conversely, for plant biomass, the overall effect of grazing was negative, although there was evidence of an increase in plant biomass resulting from increased rabbit grazing (Figure 3d). The overall positive effects of grazing on decomposition resulted from increases in historic grazing (track; Figure 3e). Overall, grazing effects were positive for C-storage, and N-availability and P-availability but negative for plant biomass and neutral for decomposition.

Effects of Aridity on Multiple Ecosystem Services
Increasing aridity was associated with reduced soil health and lower cover of woody plants, but its effects on various services were mixed (Figure 1). There were significant negative and direct effects of aridity on N-availability and C-storage (Figure 1b, c) and decomposition (Figure 1f). Increasing aridity suppressed the positive effect of soil health for N-availability and C-storage (Figure 1b-c) and plant biomass ( Figure 1e). However, increasing aridity had indirect positive effects on P-availability via reductions in woody plants and soil health (Figure 1d). In general, the sum of all direct and indirect effects or the STE of aridity was negative for four of the five services, and the only positive effect was observed for P-availability (Figure 3).

DISCUSSION
Our results indicate that aridity had a strong suppressive effect on all functions and services, except P-availability, whereas the effects of grazing ranged from positive to negative or neutral. The effect of aridity was strong, despite the relatively short gradient in our study (0·19 to 0·63). Further, different herbivores had different effects on specific services, with livestock associated with greater N-availability and decomposition, and rabbits with greater P-availability and plant biomass. Management of critical services and functions associated with plant production and C-storage, and N-availability and P-availability need to consider not only the effects of grazing by different herbivores but also the effects of changing climates on these attributes.

Grazing Effects on Plant Biomass
Plant production is the attribute most strongly influenced by livestock grazing (herbivory: Milchunas & Lauenroth, 1989;Charles et al., 2016;Eldridge et al., 2016aEldridge et al., , 2016b, so reductions in biomass with increasing grazing are expected. The STEs indicated that the overall effects of kangaroo grazing on plant biomass were largely neutral (Figure 3d), consistent with studies of kangaroo effects on soil health (Eldridge et al., , 2016b. Livestock grazing resulted in reductions in plant biomass, but increasing rabbit grazing was associated with increased plant biomass. This might at first seem counterintuitive, but rabbit grazing has been shown to suppress small forbs and grasses at the expense of large, high biomass Mediterranean weeds (Myers & Poole, 1963). Examination of sites in our study with a high intensity of rabbit grazing (i.e. sites with a mass of rabbit dung > 100 kg ha À1 ) showed that they were dominated by annual or biennial weedy plants of high biomass and cover such as Echium plantagineum L., Stipa scabra Lindl., Hordeum leporinum Link and Arctotheca calendula Levyns, typically exotic plants (Table III). Most importantly, as discussed more fully later, rabbits had a positive effect on the amount of available P in the soil (Figure 3). Soil P is the main limiting soil nutrient in the largely weathered old soils from Australia (Lambers et al., 2008); thus, by bringing back soil available P to the surface, rabbits may promote the productivity of these sites.

Grazing and Aridity Effects on C-storage
Grazing and aridity had direct and contrasting effects on Cstorage, with positive effects of grazing and negative effects of increasing aridity (Figure 1c). We also found strong, indirect negative effect of grazing and aridity via the suppression of the positive effects of soil health on this service.
The STEs for C-storage indicated that there is an overall neutral effect of grazing, but a substantial negative effect of aridity on C-storage (Figure 3b), consistent with other studies (Delgado-Baquerizo et al., 2013b;Rabbi et al., 2015). Studies of grazing exclosures have shown few differences in total C-storage after 15 years between ungrazed and grazed sites (Nosetto et al., 2006). Similarly, Rabbi et al. (2015) showed that land use in drylands in eastern Australia had a relatively neutral effect on C-storage, explained only 1·4% of the total variation in C-storage. This contrasted with aridity, which, along with soil clay content, explained 64% of the variation in C-content. Ecosystem services other than C-storage, however, are likely to be more strongly influenced by differences in land management.

Grazing and Aridity Effects on N-availability
Overall, grazing had a positive effect on N-availability, that is, increased nitrate, ammonium and total N; and this was due mainly to sheep/goat and cattle grazing (Figure 3a). Increased grazing intensity was associated with direct positive effects on N-availability, but indirect suppressive effects, via reductions in soil health (Figure 1b). The addition of nitrate and ammonium from urine and dung, particularly at livestock resting camps, could partly account for the direct effects of grazing on N-availability. Although much of the nitrate is lost via volatilisation (Powell et al., 1998), addition of dung and urine (i.e. ammonia) provides a source of readily available nutrients for plants and microbes (Augustine et al., 2003;Schrama et al., 2013). In support of this, we found that grazing intensity increased the NH 4 :NO 3 ratio ( Figure S1), although the strength of the relationship was weak (R 2 = 0·08). The mechanism underpinning the indirect suppressive effect of grazing on N-availability may relate to reductions in soil surface roughness and integrity, biocrust cover, the depth and incorporation of the litter layer (Eldridge et al., , 2016b. Biocrusts are essential components of soil health, and cyanolichens and cyanobacteria in biocrusts are known to fix atmospheric C and N, accounting for the strong effect of increased soil health on both Cstorage and N-availability, particularly in arid and semi-arid systems . In our study, increases in this ratio were linked to increases in the intensity of kangaroo grazing but reductions in cattle grazing.

Grazing Effects on P-availability
Phosphorus cycling was one of the ecosystem services less affected by increasing aridity and grazing intensity. Unlike C and N, P is an abiotically derived element, and its availability has been shown to increase when P-rich parent material is exposed, often via soil erosion (Delgado-Baquerizo et al., 2013b). Interestingly, rabbit grazing had the strongest stimulatory effect on P-availability ( Figure 3c). Production of phosphatase is extremely costly in terms of N and C. Thus, microbes and plants only produce phosphatase when it is really needed. If rabbits are providing a directly available form of P, such as PO 3 À to plant and microbes by, for example, exposing bedrock or dung, then it is likely that phosphatase production will be inhibited. Indeed, P is a relatively large component of rabbit dung, with five times more P in rabbit than in sheep or cattle dung (http://www. crossroadsrabbitry.com/rabbit-manure-info/). This explains the decoupling of phosphatase from inorganic P-availability.
Kangaroo dung is known to contribute relatively high levels of total P to floodplain systems (Kobayashi et al., 2011), but in our study, kangaroo effects on P-availability were neutral. This may be due to slow breakdown of pellets in our system compared with floodplains, owing to low levels of soil moisture (Davis & Coulson, 2016). Thus, increased P could have resulted from the concentration of rabbit dung, such as what occurs in rabbit latrines (Dixon & Hambler, 1993), or the localisation of dense patches of dung in litter dams on sloping surfaces following overland flow (Mitchell & Humphreys, 1987). The most parsimonious explanation, however, is that P-rich subsoil is exposed during the construction of the extensive communal burrow systems of rabbits. The relatively low to neutral pH values of these soils would have made P more available because there is little soil calcium to bind onto the P (Lajtha & Bloomer, 1988). Overall, therefore, increases in rabbit grazing are likely to lead to increases in P, with resulting changes in the stoichiometry of P and N in some local areas from Australia (Delgado-Baquerizo et al., 2013b; Figure S2).

Organic Matter Decomposition
We detected small declines in decomposition, our measure of enzyme activity, with increasing aridity and grazing. Aridity suppressed, and grazing slightly increased, all enzyme functions except phosphatase (Figure 4). Previous studies have shown that overgrazing reduces phosphatase and beta-glucosidase, but that effects likely depend on the  patch type in which measurements are made (Zhang et al., 2016). In our study, the decomposition service, and the four individual enzymatic functions, was associated with livestock tracks, our measure of historic grazing. Grazing is typically associated with declines in enzyme activity (e.g. Prieto et al., 2011;Olivera et al., 2014), so it is somewhat counterintuitive that increased grazing in our study was associated with increased enzyme activity. The STEs indicated that historic grazing (i.e. livestock tracks), rather than recent grazing (i.e. dung from livestock), was linked to increases in the decomposition function (Figure 3e). Our measure of historic grazing could be a proxy for increasing soil texture given that livestock tracks would be more pronounced in finer soils. Equally plausible is that increased historic grazing is linked to larger pools of herbivore urine and dung and therefore greater levels of decomposition. Also, herbivores break down organic matters such as plant litter through hoof action, and this may enhance the decomposition process. Overall, however, grazing-linked increases in decomposition were matched by declines due to increasing aridity, which would reduce decomposition rates and therefore nutrient cycling functions Delgado-Baquerizo et al., 2016b).

Stronger Negative Effects of Aridity on Services than Grazing
Aridity levels are predicted to increase into the next century and lead to substantial shifts in plant and microbial processes in drylands Maestre et al., 2016). This will likely reduce Earth's capacity to support essential ecosystem functions and services associated with the storage and availability of C and N, and the production of forage for livestock (Maestre et al., 2016). In our study, increasing aridity was associated with an increase in P-availability, but reductions in the other four functions, with the greatest reduction in Cstorage (STE = À0·54; Figure 3b) and N-availability (STE = À0·31; Figure 3a). Increasing aridity was also associated with indirect suppression of N-availability, Cstorage and plant biomass via reductions in soil health (Figure 1b-c, e). Declines in N-availability with aridity were matched by strong increases due to grazing. Predicted reductions in grazing capacity with increased aridity are therefore likely to lead to global reductions in Navailability. The shift from free-range grazing to feedlots will likely lead to reductions in C emissions and may also reduce farm-level N deposition, but the positive effects of reduced N-availability will likely be more apparent under less arid conditions (Giese et al., 2011). The effect of aridity was also to suppress the negative effects of woody plants on plant biomass. This could occur by removing competition for light or soil moisture, allelopathic effects that exist in some Eucalyptus species (Zhang & Fu, 2009), or suppression resulting from below-ground resource competition from C. glaucophylla (Harris et al., 2003). Changes in land management may not lead to increased levels of ecosystem services owing to the strong negative effect of aridity on most services. Thus, services are likely to decline over the next century as aridity increases.

Conclusions
The effects of grazing on ecosystem services are herbivore specific and vary from positive to neutral or negative. Critical functions associated with decomposition and nutrients cycling declined with increasing aridity, and these effects were of a greater magnitude than any effects due to grazing.
Our study suggests that changes in land management may fail to compensate for the negative effects of aridity on all functions other than P-availability. Thus, strategies to manage ongoing climate change are likely to be a priority of governments as we move towards a drier world.

SUPPORTING INFORMATION
Additional Supporting Information may be found online in the supporting information tab for this article. Table S1. Description of attributes and their relevance for assessing soil health. Attributes used to assess the overall measure of soil health. Figure S1. Structural equation model for NH 4 :NO 3 including the standardised total effects of aridity, measures of grazing, soil health and woody cover. Figure S2. Structural equation model for N:P including the standardised total effects of aridity, measures of grazing, soil health and woody cover.