Does grazing exclusion in Australia’s rangelands affect biomass and debris carbon stocks?

Study sites

The study assessed 46 existing grazing experiments that could be located across Australia’s rangelands (Fig. 1, Table 1). These experiments were identified using a combination of literature review and extensive consultation with land managers and scientists involved in their establishment and management. The experiments were mostly paired grazed control plots and fenced exclosures to prevent grazing by specific animals, but a subset of six experiments included paired plots with rotational grazing treatments and exclosures with control plots of continuous grazing (Table 1). The experiments had been maintained for 9–98 years (mean of approximately 40 years) at the time of measurement. For each of the exclosure plots, any breaches in the fence were recorded. Where grazed control plots had not previously been established for grazing exclosures, new plots with a known grazing history were located immediately adjacent to the exclosure area for measurement.

Fig. 1.

Locations of the study sites and associated vegetation types. The distribution of rangelands was obtained from DCCEEW (2023). Note that the ‘non-rangeland site’ was included because it had vegetation and climate considered representative of Eucalyptus rangelands.

Within the experimental sites, four main vegetation types were defined as dominated by Eucalyptus species, Acacia aneura, other Acacia species (possibly including, but not dominated by A. aneura), and other species.

Data collection

Carbon stocks in live aboveground biomass (AGB), belowground biomass (BGB) and standing dead biomass of trees and shrubs were determined at all sites, and C stocks in CWD and litter were measured at a subset of sites (n = 36 and 26 respectively).

Plot measurements were made between July 2022 and November 2023. Generally, the entire plot area was measured, except trial sites with plots of >1 ha. For these larger plots, live biomass and standing dead C stocks were measured using transects (100 m × 5 m), with the number of transects measured increasing with the variability of the soil and/or vegetation. Within each plot/transect, stem diameters of Eucalyptus, ‘Other trees’, and large Acacia trees were measured at 1.3 m height, whereas stem diameters of all other woody vegetation, including shrubs, mallees and multi-stemmed acacias, were measured at 0.1 m height. In plots/transects with high frequencies of small plants (<10 cm diameter at 1.3 m), counts of plants in each functional type (live and dead plants) in 2-cm diameter classes were taken instead of measuring the diameter of all small plants. AGB of individual woody plants (AGB_i) and live BGB of individual woody plants (BGB_i) were determined using allometric equations developed for these plant functional types (Paul et al. 2016, 2019). Standing dead AGB_i and BGB_i were determined using allometric equations developed for these plant functional types (Paul and Roxburgh 2024).

For senescing and senesced trees with stem diameters (at 1.3 m) of >55 cm, individual tree biomass estimates were reduced by multiplying the fraction remaining of aboveground structure by the biomass based on allometric relationships developed for Australian rangelands (M. Piper, pers. comm., January 2023). Specifically, dead plants were given a score that quantifies the remaining proportion of biomass relative to a live plant, as follows: 1, if it still has small branches in the crown; 0.8, if it has only large branches in the crown; 0.6, if only the main stem and few large branches remain; 0.4, if only the main bole is remaining or most of stem height remains; and 0.2, if more than half the height of the main stem is missing (M. Piper, pers. comm., January 2023).

Stand-level live AGB, BGB and standing dead mass (Mg dry matter (DM) ha⁻¹) were calculated by summing all individual live and dead biomass estimates respectively, per plot and dividing by the plot area. Stand-level basal area was calculated by summing the basal area of each tree and dividing it by the plot area. Since the height at which stem diameters were measured depended on plant function types, allometric equations (Paul et al. 2016) were used to calculate a diameter at 0.1 and 1.3 m for each tree, and hence stand basal area values for 0.1 and 1.3 m were calculated for each plot. Unless otherwise stated, stand basal area values are provided for a height of 0.1 m.

At the 26 sites where litter was sampled, at least 24 litter samples were collected from random locations within each plot/transect by using 0.0625 or 0.25 m² quadrats. All material with a diameter of <25 mm was collected, and all samples were combined into a single sample for each plot. Litter was sorted into foliage or bark/twigs/capsules from woody plants, discarding any grass and soil present in the samples. Litter was then oven dried at 70°C to a constant weight and weighed. Litter stocks were calculated in Mg DM ha⁻¹ for each sample from oven-dry litter mass and quadrat area.

CWD was measured at 36 sites by using an approach based on the line intersect method (van Wagner 1968) (Eqn 1). In plots where CWD was measured, the sum of transect lengths was at least 300 m, and CWD was defined as all woody material with a diameter of >25 mm where it crossed the transect line.

where d is the diameter (cm) of each piece of CWD where it crossed the transect, l is the length of the transect (m), g is the weighted average (by piece diameter) of the proportion of wood remaining, and ρ is the wood density (oven dry mass per pre-dried volume, Mg m⁻³). The ρ used was 0.55, which was the average of many studies from Australian ecosystems (Kingston and Risdon 1961; Attiwill 1979; Stewart et al. 1979; Gifford 2000; Ilic et al. 2000; Raymond and Muneri 2000; Miranda et al. 2001; Raymond and Muneri 2001; Slijepcevic 2001; Bridges 2005; Searle and Owen 2005; Ximenes et al. 2008; Callister and England 2010; Sillett et al. 2010; Ximenes et al. 2016; Falster et al. 2021). Dry-mass values were converted to carbon mass by multiplying live or standing dead mass by 0.5 (Gifford 2000; Mathers et al. 2003; Possell et al. 2015; Neumann et al. 2020), CWD by 0.49 (Stewart and Burrows 1994; Chee 1999; Woldendorp 2000; Harms and Dalal 2003; Roxburgh et al. 2006b; Hilbig 2009; Volkova and Weston 2013; Aponte et al. 2014; Possell et al. 2015), and litter by 0.51 (Hutson 1985; Chee 1999; Gifford 2000; Woldendorp 2000; Russell-Smith et al. 2009; Possell et al. 2015; Volkova and Weston 2015).

Climatic data for each site were obtained from Queensland Government’s Silo database (www.longpaddock.qld.gov.au/silo/) by using the R package cropgrowdays (ver. 0.2.0.; see https://CRAN.R-project.org/package=cropgrowdays) and R software (R Core Team 2023). Daily precipitation, minimum and maximum temperatures, and evapotranspiration were used to calculate monthly and annual values. For example, the mean annual precipitation (MAP) for each site was derived by taking the sum of all the annual values since the exclosure was established, divided by the number of years since the exclosure was established. Experimental sites were initially defined as arid (≤300 mm MAP) or non-arid >300 mm MAP. In Australia, arid zones are often defined as areas with ≤250 mm MAP (DCCEEW 2024); therefore, this was used as an additional alternative cut-off for arid and non-arid zones in subsequent analyses. Further, an aridity index defined by Eamus et al. (2000) was also derived for each experimental site. This index was the ratio of monthly mean daily pan evaporation rate to monthly rainfall. When monthly rainfall was zero, a value of 0.1% of annual rainfall was used for that month to avoid division by zero.

Statistical analysis

Differences among vegetation types or between grazed and control treatment means were examined using linear mixed models. A natural logarithmic transformation was applied to each C stock so that the residuals were centred at zero and approximately normally distributed. Treatment effects were examined when combining data, as well as splitting the data into the four vegetation types (Eucalyptus species, Acacia aneura, other Acacia species and Other), and/or into different regions (e.g. Queensland, New South Wales, Northern Territory, Western Australia, South Australia, and groups of these), and into the individual sites that contained replication. This was performed while including all data within these groups, as well as after excluding sites where exclosures had been breached at some point in time since establishment, and when dividing the dataset between arid and non-arid MAP. When testing for differences between treatments at individual sites that contained replication, the replicate was used as a random effect. These were fitted using the R package nlme (ver. 3; see https://CRAN.R-project.org/package=nlme; Pinheiro et al. 2018) and R software (R Core Team 2023).

For the total dataset, linear mixed models were used to test for correlations between C stocks and climatic variables, or the time since exclosure. ‘Site’ was included as a random effect. Linear mixed models were also used to test whether the relative difference between the C stock of a given exclosure treatment, and the grazed control plot at the same site, was correlated with climate variables, the time since exclosure or stand basal area of the grazed control plots.

To complement the mixed models, a non-parametric test was used to test whether a given exclosure treatment differed from the grazed plots. Using the available pairs of exclosure plots and grazed plots at each site or replicate, the observed test statistic was calculated as the proportion of pairs where the C stock of the exclosure was greater than that of the grazed plot. The null hypothesis was that there was no difference between the exclosures and grazed plots. The observed test statistic was compared with the proportion of plots where the exclosure C stocks were greater than the grazed plots after randomly assigning the treatment (exclosure or grazed) to each plot. This was repeated 100,000 times to calculate the probability of the observed test statistic being greater or smaller (i.e. two-tailed test) than the statistic calculated after randomly assigning the treatments, which was a distribution of probabilities centred on 0.5. As for the mixed models, this test was performed when combining data, as well after splitting the data into groups relating to vegetation type and region.

We also performed a meta-analysis on the effects of grazing exclusion using the R packages metafor (ver. 4.6-0; see https://CRAN.R-project.org/package=metafor; Viechtbauer 2010) and foresplot (ver. 3.1.3; see https://gforge.se/packages/url), following the methods described by Borenstein et al. (2009). R codes for the meta-analysis were based on those described in Zhang et al. (2020). The control mean (Cmean), treatment mean (Tmean), their corresponding standard deviations (Csd and Tsd), and sample size (n) were calculated and are provided in the Supplementary material (Tables S13-S26). Overall, 27 of the 46 grazing sites were unreplicated in the grazed and/or exclosure treatments. For treatments without replicates, the biomass-C or basal area was multiplied by the average coefficient of variation derived from the rest of the dataset for the parameter where standard deviations were available (van Groenigen et al. 2011). We used multilevel random effect meta-analysis to account for variable dependencies, such as when an experiment had rotation and exclusion treatments that shared a set of common controls. The statistical analysis was conducted on natural log-transformed response ratios (ln RR), where the ratio was the ungrazed exclosure mean divided by the grazed control mean. The effect sizes were back-transformed from ln RR to response ratio on the natural scale (RR) to aid interpretation. Treatment was considered to represent a significant (P < 0.05) increase or decrease compared with the control when the 95% CI of the response variable (i.e. AGB-C, BGB-C, basal area, CWD-C or litter C) did not overlap with zero.

C stocks in relation to vegetation type and climate

Most C stocks (except litter C) were significantly higher in the Eucalyptus plots than in the ‘other’ plots (Fig. 2). Biomass C stocks were positively correlated with mean annual precipitation (Supplementary Fig. S1) and potential evaporation (Fig. S3), but not correlated with mean maximum temperature in January (Fig. S2), usually the hottest month. CWD C and litter C followed the same trends. The low pseudo-R² (mR²) but relatively high conditional pseudo-R² (cR²) suggests that much of the variability in C stocks was explained by the variability within sites rather than the differences in climate among sites (Figs S1, S3).

Fig. 2.

Box plots comparing (a–e) the C stocks and (f) basal areas (at height of 0.1 m) of different vegetation types when combining data from exclosures and grazed control plots. The y-axes are ln-transformed. Different letters indicate significant differences between individual vegetation types (P < 0.05). AGB, aboveground biomass; BGB, belowground biomass.

Influence of domestic grazers on C stocks

There were no significant grazing effects on C stocks when combining all vegetation types (Fig. 3), or when considering individual vegetation types (Supplementary Table S12). Nor were there any significant differences between treatments at individual sites that contained replication (except for litter C stocks at Site 19, and several stocks at Site 10, Table S12). Similarly, there were no significant differences when focusing on individual vegetation types within a region, such as A. aneura in Queensland + New South Wales (except for litter C stocks), or Acacia and A. aneura in the Northern Territory (Table S12). Whether the grazers were cattle, sheep or combinations of these also had no significant effect (Table S12). The treatment effects were also not changed by excluding all sites where exclosures had been compromised at some point during the experimental period (Table S12). The non-parametric tests were usually consistent with the results from the mixed models (Table S12).

Fig. 3.

Box plots comparing (a–e) the C stocks and (f) basal areas (at height of 0.1 m) of exclosures and grazed control plots. The y-axes are ln-transformed. AGB, aboveground biomass; BGB, belowground biomass.

Results from the meta-analysis showed high variability and low frequency of significant grazing effects on biomass parameters across the vegetation types (Fig. 4, Table S1). Although mean response ratios were often >1, indicating that grazing exclusion increased the response variable relative to the grazed control treatment, owing to the high variability among grazing response observations, the 95% confidence intervals for response ratios generally included one for AGB-C, BGB-C and CWD-C, whether considering each vegetation type, or averaged across vegetation types. Basal area was higher, on average between 1.08 and 1.59 times, with grazing exclusion, but a consistent increase in basal area with grazing exclusion was observed only in A. aneura-dominated vegetation. However, there was some indication of a consistent trend whereby the mean response ratio in Eucalyptus-dominated vegetation was more positive with grazing exclusion than in the other vegetation types for AGB-C and BGB-C. Total litter C increased on average between 1.16 and 1.96 times with grazing exclusion. The response was largest within Acacia aneura-dominated vegetation, with litter C increasing between 1.27 and 3.08 times with grazing exclusion. When excluding plots from the analysis where exclosure fences had some damage, the trends remained similar.

Fig. 4.

The impact of grazing exclusion on the response variables aboveground biomass carbon (AGB-C), belowground biomass carbon (BGB-C), basal area (at height of 0.1 m), coarse woody debris carbon (CWD-C) and litter C by vegetation type. Blue lines represent the 95% confidence intervals. Red lines show the 95% confidence intervals for the overall effect of exclosure on the response variables. A response ratio of one indicates that grazing exclusion did not affect the response variable measured, whereas a ratio greater than one indicates that grazing exclusion increased the response variable relative to the grazed control treatment. The I² statistic estimates the amount of variability in effect size that can be attributed to heterogeneity among true effects (genuine variation among experiments). I² indices of zero indicate that the variability within experimental trials was high.

In non-arid regions (MAP > 300 mm), AGB-C, BGB-C, basal area and litter C stocks were consistently higher where grazing was excluded (Fig. 5, Table S4). The 95% CIs were generally large. When the impact of grazing exclusion was evaluated using 250 mm as a cut-off for defining aridity, these trends were generally weaker or no longer observed (Fig. S6, Table S5). However, the impact of grazing exclusion evaluated using an aridity index defined by Eamus et al. (2000) found similar trends to the 300 mm precipitation cut-off (Fig. S7, Table S6).

Fig. 5.

The impact of grazing exclusion on aboveground biomass carbon (AGB-C), belowground biomass C (BGB-C), basal area (at height of 0.1 m), coarse woody debris C (CWD-C) and litter C by arid and non-arid regions as defined by mean annual precipitation above and below 300 mm. Blue lines represent the 95% confidence intervals. Red lines show the 95% confidence interval for the overall effect of exclosure on the response variable. A response ratio of one indicates that grazing exclusion did not affect the response variable measured, whereas a ratio greater than one indicates that grazing exclusion increased the response variable relative to the grazed control treatment. The I² statistic estimates the amount of variability in effect size that can be attributed to heterogeneity among true effects (genuine variation among experiments). I² indices of zero indicate that the variability within experimental trials was high.

Across the total dataset, relative differences between treatments were highly variable (Fig. 6). As implied by the low frequency of significant effects mentioned above, most relative differences were close to zero or centred around zero; however, there were many pairs (counts provided in Fig. 6) with relative differences >0 (exclosure C stock > grazed C stocks). These relative differences in C stocks were not correlated with any climate variables, or the duration of the experiment on the basis of linear mixed modelling (P > 0.05, data not shown). However, the relative differences in live C stocks declined as the basal area of the grazed control plots increased (Fig. S4).

Fig. 6.

Frequency distributions showing the relative differences between ungrazed plots and grazed plots of (a–e) C stocks or (f) basal areas (at height of 0.1 m), i.e. (Cstock_ungrazed – Cstock_grazed) / Cstock_grazed. For sites with multiple exclosure treatments, each treatment is considered (not the average from each site). The legends show the number of ungrazed plots where relative differences were either >0 or <0. AGB, aboveground biomass; BGB, belowground biomass; CWD, coarse woody debris.

Additional random effect meta-regression was used to assess whether the effect sizes of grazing exclusion on C stocks were correlated with climate variables (mean annual precipitation, mean maximum temperature in January, mean annual evaporation), the duration of the experiment or latitude. These results were generally in agreement with the linear mixed modelling. The 95% CIs of the response gradients or slopes of the regression were generally large and were a mix of negatives on the lower limit and positives on the upper limit (Figs S7–S11, Tables S7–S11), indicating that there were no clear linear relationships between vegetation response to grazing exclusion and climatic variables, latitude and age of the experiment within the ranges we have investigated. The exceptions were for AGB-C and BGB-C in response to age of experiment (Fig. S10, Table S9). For every additional year of experiment, total AGB-C increased in response to grazing exclosure with the response gradient of 0.013 [95% CI: 0.003, 0.023], corresponding to a response ratio of 1.013 [1.003, 1.023]. This equates to a mean effect of between ~0.3% and 2% increase in AGB-C with every year of grazing exclusion. For every additional year of experiment, BGB-C increased in response to grazing exclosure with the response gradient of 0.014 [0.004, 0.023], corresponding to response ratio of 1.014 [1.004, 1.024] or a mean effect of between 0.4% and 2.4% increase in BGB-C with every additional year of experiment. The trends remained the same after excluding the experimental plots with damaged fences from the meta-regression analyses (data not shown).

Influence of domestic grazers and rabbits or macropods on C stocks

Across the total dataset, although treatments where rabbits and domestic livestock were excluded often had higher mean live AGB-C, BGB-C and CWD-C stocks than did domestic livestock exclusion alone, increases were significant only for basal area (Fig. 7). Exclosures where rabbits were excluded (in addition to domestic grazers) had significantly higher live above- and belowground C at the Old Mans Plains site, and for the Acacia vegetation type (Table S12). There were no clear trends when considering other vegetation types or states (e.g. Western Australia, Queensland) individually (Table S12). At the Boorara site, exclosures where domestic animals + macropods were removed had lower live above- and belowground C stocks than when only domestic animals were excluded. For the Eucalyptus vegetation type, basal area was highest in exclosures where domestic animals + macropods were excluded. There were no significant effects of excluding domestic animals and macropods (Fig. 8). These results are consistent with the meta-analysis results (Fig. 9, Table S2). When the types of grazers excluded were considered, grazing exclusion significantly affected only litter C and basal area. Specifically, litter accumulation was seen in experimental sites where only domestic animals were excluded (1.41–2.34). Basal area was highly variable but was on average between 1.05 and 3.66 times higher with grazing exclusion when domestic animals, macropods and rabbits were excluded, but this consistent positive trend was lost when we excluded exclosures that had been compromised at some point during the experimental period (data not presented). The trends remained similar for the rest of the variables when we excluded plots where the exclosure that had been compromised at some point during the experimental period.

Fig. 7.

Box plots comparing (a–e) the C stocks and (f) basal area (at height of 0.1 m) of grazed control plots and exclosures with or without rabbits. The y-axes are ln-transformed. Letters indicate significant differences between grazing and exclosure treatments. AGB, aboveground biomass; BGB, belowground biomass; Dom exc, domestic animals excluded; Dom + rabbits exc, domestic animals and rabbits excluded.

Fig. 8.

Box plots comparing (a–e) the C stocks and (f) basal area (at height of 0.1 m) of grazed control plots and exclosures with or without macropods. The y-axes are ln-transformed. Letters indicate significant differences between grazing and exclosure treatments. AGB, aboveground biomass; BGB, belowground biomass; Dom exc, domestic animals excluded; Dom + macropods exc, domestic animals and macropods excluded.

Fig. 9.

The impact of grazing exclusion on aboveground biomass C (AGB-C), belowground biomass C (BGB-C), basal area (at height of 0.1 m), coarse woody debris C (CWD-C) and litter C by the type of grazer excluded. Blue lines represent the 95% confidence intervals. Red lines show the 95% confidence interval for the overall effect of exclosure on the response variable. A response ratio of one indicates that grazing exclusion did not affect the response variable measured, whereas a ratio greater than one indicates that grazing exclusion increased the response variable relative to the grazed control treatment. The I² statistic estimates the amount of variability in effect size that can be attributed to heterogeneity among true effects (genuine variation between experiments). I² indices of zero indicate that the variability within experimental trials was high.

Influence of vegetation and climate on C stocks

Eucalyptus-dominated vegetation had the highest live AGB-C stocks (mean of 24, and maximum of 115 Mg C ha⁻¹), and these sites also had higher annual precipitation and lower temperatures than did the remaining vegetation types (Figs S1, S2). AGB-C stocks from other studies in Eucalyptus-dominated rangeland vegetation were 6–80 Mg C ha⁻¹ (Burrows et al. 2002; Cook et al. 2005; Harms et al. 2005). The Acacia-dominated vegetation type averaged the lowest live AGB-C stocks (mean 1.8, and maximum of 8.6 Mg C ha⁻¹) and these were generally on the drier sites.

The ‘other’ and A. aneura vegetation types averaged intermediate live aboveground C stocks with means of 2.5–6.4 Mg C ha⁻¹ and maximum of 8.5–30.6 Mg C ha⁻¹ respectively. Previous studies have reported up to 49 Mg C ha⁻¹ for A. aneura (Witt et al. 2009b; Witt et al. 2011; Daryanto et al. 2013; Williams et al. 2024). The correlations between live C stocks and annual precipitation and evaporation observed here (Figs S1, S3) are consistent with the strong influence of precipitation on biomass and soil C stocks in Australian rangelands (Fensham et al. 2012; Allen et al. 2013; Bray et al. 2016; Henry 2023).

Mean CWD-C was 1.4 Mg C ha⁻¹, with a maximum of 9.4 Mg C ha⁻¹ in A. aneura vegetation. This may reflect the practice of felling or pruning A. aneura to feed livestock during dry periods or tree clearing or thinning to enhance forage growth, which can increase CWD (Pringle et al. 2021). Previous studies in rangeland vegetation dominated by Eucalyptus, Acacia, A. aneura, Callitris or shrubs have recorded a mean CWD of about 3.4 C Mg ha⁻¹, and a maximum of 36 Mg C ha⁻¹ (Burrows 1976; Woldendorp et al. 2002; Murphy et al. 2003; Harms et al. 2005; Daryanto et al. 2013; Prior et al. 2017).

The litter C pool was generally the smallest of the C pools examined, ranging from 0.005 to 2.31, with a mean of 0.41 Mg C ha⁻¹, and did not vary among vegetation types. These values are comparable to those in other Australian rangeland studies, with litter C ranging from close to zero up to 8.7 Mg C ha⁻¹ (Burrows 1976; Roxburgh et al. 2006a; Fontaine et al. 2012; Bray et al. 2014; Hunt 2014; Dalgleish et al. 2015; Williams et al. 2024).

Influence of grazing on C stocks, and interactions with climate

A key finding was that relative differences between grazing treatments were highly variable, and as a result of that variability, were generally not significant. The variability in responses to grazing found here is consistent with that in previous studies. For example, a higher AGB-C was reported for a cohort of trees with a diameter of <15 cm in Queensland A. aneura experiments, whereas larger cohorts (>15 cm) contained more biomass in the grazed treatments, leading to similar total AGB-C (Witt et al. 2009b; Witt et al. 2011). Tree densities were higher within exclosures at the Kidman Springs experiment (Bastin et al. 2003). However, no grazing effects on AGB-C were found in the central Australian Acacia experiment at Central Mount Wedge (Foran et al. 1982) or the Wapweelah A. aneura experiment in western New South Wales (Daryanto et al. 2013). Similarly, there was little difference in woody biomass between moderate and heavily grazed plots in the Wambiana experiment in Queensland (Eucalyptus vegetation) (Bray et al. 2014).

As the basal area of the grazed control plots increased, the increase in live C stocks of exclosures relative to grazed controls decreased (Fig. S4). This might be expected because when live C stocks are already high, the C stocks may be limited more by soil and climatic conditions than grazing, and therefore, removing grazing animals may have a lesser influence.

Meta-regression analysis suggested that response ratios for AGB-C and BGB-C consistently increased linearly with time since the exclosures were established (Fig. S10), suggesting that significant effects of grazing exclusion are more likely after long durations of grazing exclosure. However, the linear mixed model analysis did not find this correlation. The incongruence suggests uncertainties driven by the lack of replication in 27 of 46 grazing experiments.

The influence of grazing on woody vegetation may be masked by the effects of other factors that have greater influences on woody vegetation in rangelands, such as precipitation, temperature or pre-existing biomass stocks (Fensham et al. 2005; Fensham et al. 2012). Temporal changes in precipitation can have greater effects than grazing, at least on non-woody vegetation (Foran et al. 1982) or on woody vegetation cover (Fensham et al. 2005). Here, we found no significant linear relationships between site-level climate averages and biomass response to grazing exclusion, but AGB-C, BGB-C, basal area and litter C were consistently, albeit variably, higher with grazing exclusion in non-arid regions with MAP of >300 mm. Similarly, the mean response ratio for the Eucalyptus-dominated vegetation type, which experienced the highest MAP, was more positive than that in other vegetation types for AGB-C and, BGB-C with grazing exclusion. Given that regeneration is facilitated by precipitation (e.g. Hodgkinson 1991), it is likely that there are more opportunities for woody vegetation to regenerate and survive in regions with MAP of >300 mm.

Exclusion of domestic grazers can increase soil water infiltration and, subsequently, soil water-holding capacity (Witt et al. 2011), which could potentially enhance woody vegetation growth and recruitment while reducing mortality rates. The presence of grasses has been associated with lower water stress and higher recruitment of A. aneura than in locations where grasses were heavily grazed (Anderson and Hodgkinson 1997). Although this example may not hold for all vegetation types and locations, given that precipitation and soil water availability can have a strong temporal and spatial influence on rangeland vegetation (Foran et al. 1982; Fensham et al. 2005; Fensham et al. 2012), combining C stock measurements with additional data on soil moisture availability may help explain some of the observed variability in grazing effects. Furthermore, grazing intensity may modify the effect of a broad ‘grazing’ category on woody plants in arid and semi-arid rangelands. Without more detailed data, distance to water point could be used as a proxy for grazing intensity in future studies, although changes in water infrastructure during the experimental period would need to be considered. As an example, sheep grazing in a chenopod-dominated rangeland increased adult shrub mortality close to water and reduced recruitment, resulting in a decline in the shrub population, with the largest decline closer to water (Hunt 2010).

We were unable to obtain data on the stocking density of grazing animals during the experimental period, either domestic livestock or other grazers, including macropods and feral animals. However, six studies had rotational grazing treatments in addition to total exclusion plots. These studies showed more negative or less positive trends in biomass C in rotational grazing treatments than in grazing exclusion (Fig. S6, Table S3). The stocking density of grazing animals can be difficult to obtain (Allen et al. 2013), but would allow for more accurate measures of grazing intensity than the discrete treatments used here. This should include the stocking density of non-domestic animals, such as rabbits, feral goats and macropods, which may have a strong confounding influence or even a major influence on rangeland biodiversity in some regions (Mutze 2016).

Other factors that could be considered include fire and the palatability and accessibility of tree and shrub species for grazing at various stages of their lifecycle. It is worth noting that grazing can lead to stunted growth forms of some species, such as A. aneura (e.g. ‘cattle mulga’ Everist 1949), which could also potentially influence biomass C stocks. However, no such effects on A. aneura were observed in the 46 grazed A. aneura plots at 23 sites within the paired-plot network.

Implications

This study shows that responses to grazing exclusion are variable and difficult to predict, not that there are no grazing effects on C stocks. The effects of grazing exclusion on most C stocks ranged from moderate reductions to large increases in C stocks or basal area following grazing exclusion (Figs 4, 6). AGB-C, BGB-C, basal area and litter C were significantly, albeit variably, higher with grazing exclusion in regions with annual precipitation of >300 mm. The large variability within and across sites, as well as the small sample sizes, and the fact that woody vegetation C stock sizes are small in rangelands, all contributed to the difficulty in quantifying grazing effects in rangelands. It is important to note that although a sample size of 46 experiments includes all existing long-term exclosure experiments that were still intact and could be located and measured, it is a small sample relative to the large area of Australian rangelands.

Furthermore, the underlying drivers of the variability in C stocks could not be identified. As discussed above, future studies may benefit from combining C stock data with information on soil moisture availability, total grazing pressure, including grazing intensities of major native and feral animals such as macropods and rabbits, and repeated measurements. This information, although difficult to source broadly across the rangelands, could provide the foundations for a model capable of predicting where and when grazing effects might be significant and their impact on regional C stocks. The impact of fire, if grazing is removed from large areas, should also be considered. Given that rangelands cover 80% of the Australian continent, even if grazing has a small influence per area, or has an influence only in specific locations, altering grazing could have large implications for regional and national C accounting, warranting further investigation.

It is also important to note that maximising C stocks may not necessarily be the most ecologically sustainable objective (Witt et al. 2009b; Witt et al. 2011). For example, higher stand densities in terms of leaf area can increase water use (Chen et al. 2014) and thereby reduce soil water availability, potentially exacerbating the effects of low precipitation, and create ‘boom–bust’ cycles in response to climate variability. Increases in woody C stocks nearly always lead to reductions in pasture biomass, thereby negatively influencing livestock production (McIvor and Gardener 1995; Jackson and Ash 1998; Scanlan 2002; Dyer and Stafford Smith 2003; Witt 2013; Cowley et al. 2014).