Soil carbon storage capacity of drylands under altered fire regimes

The determinants of fire-driven changes in soil organic carbon (SOC) across broad environmental gradients remains unclear, especially in global drylands. Here we combined datasets and field sampling of fire-manipulation experiments to evaluate where and why fire changes SOC and compared our statistical model to simulations from ecosystem models. Drier ecosystems experienced larger relative changes in SOC than humid ecosystems—in some cases exceeding losses from plant biomass pools—primarily explained by high fire-driven declines in tree biomass inputs in dry ecosystems. Many ecosystem models underestimated the SOC changes in drier ecosystems. Upscaling our statistical model predicted that soils in savannah–grassland regions may have gained 0.64 PgC due to net-declines in burned area over the past approximately two decades. Consequently, ongoing declines in fire frequencies have probably created an extensive carbon sink in the soils of global drylands that may have been underestimated by ecosystem models. Fire impacts soil organic carbon stocks, in addition to aboveground biomass, yet changes are not well constrained. This study shows that more soil carbon is lost from drier ecosystems than humid ones and that the carbon sink is increasing in savannah–grassland regions with declining burned area.

The determinants of fire-driven changes in soil organic carbon (SOC) across broad environmental gradients remains unclear, especially in global drylands. Here we combined datasets and field sampling of fire-manipulation experiments to evaluate where and why fire changes SOC and compared our statistical model to simulations from ecosystem models. Drier ecosystems experienced larger relative changes in SOC than humid ecosystems-in some cases exceeding losses from plant biomass pools-primarily explained by high fire-driven declines in tree biomass inputs in dry ecosystems. Many ecosystem models underestimated the SOC changes in drier ecosystems. Upscaling our statistical model predicted that soils in savannah-grassland regions may have gained 0.64 PgC due to net-declines in burned area over the past approximately two decades. Consequently, ongoing declines in fire frequencies have probably created an extensive carbon sink in the soils of global drylands that may have been underestimated by ecosystem models.
Fire-driven changes in soil organic carbon (SOC) arising from altered fire frequencies are hypothesized to be predicted by how much fire directly combusts SOC 1,2 and indirectly alters plant biomass inputs to soils and decomposition of residual SOC post-fire [3][4][5][6][7][8][9][10] . In drier and warmer ecosystems, which dominate global burned area 11,12 , most SOC is in the mineral horizon where heat rapidly dissipates 13 and little direct combustion of SOC occurs 8,14,15 . In these drier sites, fire-driven shifts in plant biomass inputs, especially from trees [16][17][18] , are thought to determine changes in SOC stored in the mineral horizon [19][20][21] . Consequently, increases in fire frequency may drive large SOC losses in climates with low precipitation and/or seasonal rainfall, where water constrains tree growth and post-fire biomass recovery [22][23][24][25] relative to ecosystems in climates where biomass recovery is faster. In addition to water availability, temperature and soil texture and mineralogy can modify post-fire decomposition rates 26-28 such that warmer climates and coarse-textured soils may allow for higher C losses because the residual plant material is more quickly decomposed 4 . Thus, we hypothesize that water availability, temperature and soil texture all act to modify the effect of repeated burning on SOC storage in the mineral horizon.
Global data to evaluate these hypotheses are lacking because there have yet to be studies examining repeated burning effects on SOC and plant biomass in parallel across broad climatic and ecological gradients, despite comparisons within individual ecosystems 2,19,29 . Thus, models used to simulate the effects of fire regime changes on ecosystems, such as fire-enabled Dynamic Global Vegetation Models 30 (DGVMs), lack a clear benchmark for evaluating how well they simulate SOC responses across environmental gradients 31 . Here we examine the factors that determine the magnitude of SOC losses or gains in the mineral horizon when fire frequencies change, evaluate whether DGVMs capture spatial patterns in fire effects on SOC storage and estimate the potential Article https://doi.org/10.1038/s41558-023-01800-7 found that while fire reduced SOC in semi-arid and dry sub-humid zones by 35 ± 9% and 23 ± 15%, respectively, it did not significantly decrease it in humid and hyper-humid zones (humid = 19 ± 35% lower and hyper-humid = 8 ± 14% lower, respectively; Fig. 1b).
Annual precipitation seasonality was the second most important environmental variable in the statistical model, with twice as large fire-driven declines in high vs mean seasonality sites (27 ± 8% lower vs 13 ± 9% lower, respectively, p < 0.001, at seasonality values of 69 vs 47, respectively; Fig. 1c and Supplementary Table 2; Pearson correlation coefficient between annual precipitation and aridity was relatively low (σ = 0.34)). Relative SOC losses were also greater in sites with cooler temperatures and coarser textured soils (p < 0.01 and p = 0.068, respectively), although these environmental variables were less important than aridity and precipitation seasonality according to the model selection analyses (Supplementary Table 2 and Extended Data Fig. 2). Taken together, water availability was the most important environmental factor explaining the relative change in SOC with altered fire frequencies.
To test the hypothesis that fire-driven changes in SOC could be attributed to changes in tree biomass inputs across sites, we focused on savannah-grasslands and analysed 74 plots across seven sites in our meta-analysis with data on soil δ 13 C (Methods). We used δ 13 C as a proxy for tree biomass inputs in these sites because C3 tree biomass has a lower 13 C than C4 grass biomass. Thus, SOC 13 C is commonly used to quantify C3 tree biomass inputs relative to C4 grass inputs in savannahs 19,20,34 .
Fire-driven changes in 13 C illustrate that larger decreases in SOC in drier climates were linked with lower tree biomass inputs to soils. Comparing δ 13 C values across sites and in unburned vs burned plots illustrated that 13 C was higher in burned plots (F 1,40 = 50.9, p < 0.001, mixed-effects model with site as random intercept; F 1,40 , F statistic with degrees of freedom), illustrating frequent burning caused a shift from C3 tree-to C4 grass-derived biomass inputs: the proportion of SOC from C3 trees was, on average, 39 ± 27% lower in frequently burned plots relative to unburned plots (Supplementary Table 3). The losses of C3-derived inputs positively correlated with the losses in SOC stocks across sites, pointing to changes in woody biomass inputs to soils driving changes in SOC storage, but the relative magnitude of change varied across sites (r 2 = 0.71, p = 0.01; Fig. 2a and Supplementary Table 3). In contrast, total SOC stocks from C4-derived inputs were unchanged by fire across sites (p > 0.50). Thus, while grass biomass inputs to SOC were impact of recent regional changes in fire frequencies on SOC storage.
To evaluate the determinants of fire effects on SOC, we conducted a meta-analysis to identify the environmental variables relating to how multi-decadal alterations of fire frequency impact SOC storage in the mineral horizon using data from experiments in 53 sites containing 434 replicate plots. Within these sites, we compared the effect of repeated burning at different frequencies relative to unburned plots or plots burned at lower frequencies over the same period (Supplementary  Table 1 and Methods). We focused our analyses on ecosystems that account for the majority of both total burned area and recent changes in fire frequency (savannahs, grasslands and seasonal woodlands and forests) 12 (Extended Data Fig. 1).
Globally, our meta-analysis demonstrated that the most important climatic and edaphic variables explaining fire effects on SOC in the mineral horizon were aridity, precipitation seasonality, mean annual temperature and silt content, with larger relative changes in SOC in drier and cooler environments on coarsely textured soils (r 2 = 0.82, p < 0.001; Supplementary Table 2 and Extended Data Fig. 2). Using Akaike information criterion-based model selection on mixed-effects meta-regression models, we identified the important environmental variables based on their inclusion in the top model and their relative importance calculated from summing the weights of the models in which the variable occurred. We then fit the top model to the data to illustrate the influence of these environmental variables on fire effects (Supplementary Table 2 and Methods) and used the first and third quartiles of the covariates to compare 'low' and 'high' values, respectively, of environmental covariates. Variables related to experimental design (for example, duration that fire frequencies were altered and the absolute fire frequency) and overall ecosystem type were also incorporated to better isolate environmental variables (Extended Data Fig. 2 and Methods).
Relative to unburned plots, SOC in burned plots was 17 ± 10% lower in sites with mean aridity and 37 ± 23% lower in sites with high aridity (p < 0.001, aridity index = 0.63 vs 0.31, respectively; Fig. 1a  Article https://doi.org/10.1038/s41558-023-01800-7 robust to fire treatment, declines in woody plant inputs determined the magnitude of SOC losses. Water availability was important in explaining the variability in 13 C changes across sites. Within these seven savannah-grassland sites, aridity and mean annual precipitation strongly covaried (σ = 0.95), and the limited sample size restricted our ability to conduct model selection including multiple variables. Consequently, we present analyses for aridity but point out that relationships are equivalent for mean annual precipitation and thus refer to the gradient as one of 'dryness'. The driest sites experienced the strongest fire-driven declines in C3-derived inputs (aridity: r 2 = 0.58, p = 0.029; mean annual precipitation: r 2 = 0.64, p = 0.019; Fig. 2b), consistent with fire causing the largest relative changes in SOC in drier climates. Thus, fire-driven changes in tree biomass inputs into soils helps explain SOC responses to fire across sites.
To further assess the relative changes in SOC across fire frequency treatments and attribute the changes to shifts in biomass inputs, we conducted a field-sampling campaign across the six longest-running fire-manipulation experiments (that experienced 53-64 years of altered fire frequencies) in savannahs that span semi-arid to humid zones (Methods). SOC was on average 26% lower (range of 13-44% lower) in the highest-frequency plots relative to the fire exclusion plots (F 1,34.1 = 9.1, p = 0.005; Fig. 3; 0-20 cm depth). Across all sites, fire reduced the proportion of SOC derived from C3 plants (F 1,34 = 35.4, p < 0.001). Finally, within a site, SOC in a plot positively correlated with 13 C, illustrating larger SOC stocks had a greater proportion of SOC derived from C3 plants (mixed-effects model, F 1,40.4 = 9.4, p = 0.004).
To evaluate estimates of fire effects on SOC at the global scale, we analysed whether an ensemble of seven fire-enabled DGVMs were able to recreate the biogeographical trends in fire effects found in our empirical findings (model details in refs. 30,31 and Methods; Extended Data Fig. 4 for global maps). Fire-vegetation models are able to represent spatial and seasonal patterns of burnt area, and the inclusion of fire improves the simulated vegetation distribution 30,31 . We used model experiments that were similar in concept to our field experiments: comparing simulations of SOC in a 'world with fire' (allowing fire frequency and impact to emerge as a function of climate, vegetation type and fuel load) vs a 'world without fire' (fire modules are turned off). Rather than compare model-based estimates of SOC fluxes with data at individual sites, we compare the within-model relationships between fire effects and water availability gradients with the predictions from our empirical model across the same gradients. Generally, the DGVMs predicted that areas experiencing the largest differences in burned area between the c, Scatter plots between the mean SOC and proportion from C3 plants within a plot within a site. Lines are linear regressions within a site. For a and b, the box boundary represents the interquartile range, the centre line represents the median and the whiskers represent the minimum and maximum values. Sites span a broad climate gradient and are ordered by mean annual precipitation: IBGE = hyper-humid (n = 7), CdrC = humid (n = 9), KrgP = semi-arid (n = 12), KrgrSk = semi-arid (n = 12), KrgrSt = semi-arid (n = 9), KrgM = semi-arid (n = 12).
Article https://doi.org/10.1038/s41558-023-01800-7 fire vs no fire simulations also experienced the largest changes in SOC (Extended Data Fig. 5). However, the models were inconsistent in their simulated patterns of SOC sensitivity to fire across aridity and precipitation seasonality gradients ( Fig. 4 and Extended Data Figs. 6 and 7). In savannah-grasslands, which account for 70% of global burned area, only two of the seven models (LPJ-GUESS-BLAZE and CTEM models) correctly recreated the empirically determined relationships between the sensitivity of SOC to fire and both precipitation seasonality and aridity ( Fig. 4 and Extended Data Figs. 6 and 7). Although we cannot isolate the role of model-based differences in simulating burned area across climates because we did not run constrained simulations with a forced fire frequency, these results suggest a DGVM ensemble is probably biased towards underestimating fire-driven changes in SOC in drier regions.
To estimate the potential area over which frequent burning could limit SOC storage in mineral soils, we scaled up our statistical model of fire effects on SOC to savannah-grasslands. To estimate what areas may be either losing or gaining SOC, we extrapolated observed trends in burned area for approximately two decades 12 to identify areas of increasing or decreasing fire frequency and used the environmental covariates and SOC content derived from other global maps (Methods) to estimate potential SOC changes. Across 2.3 million km 2 where burned area is tending to decline, SOC has potentially risen by 23% (Fig. 5a). In 1.38 million km 2 where burned area is tending to rise, SOC has potentially declined by 25% (Fig. 5b). The causes of the changes in burned area are described in detail in other studies 12 , but namely arise from fire suppression due to population expansion and landscape fragmentation in savanna-grassland regions. By multiplying these relative values with total SOC stocks, we estimate reductions in burning from 1998-2015 resulted in a gain of 1.78 PgC, while more frequent fires resulted in a loss of 1.14 PgC, for a net change of 0.64 PgC, or a flux of 0.038 PgC yr −1 .
While previous research has highlighted the theoretical capacity of savannah-grassland soils to serve as a C sink 21 , subsequent studies have argued that variance limits broad extrapolations; our study, for the first time, identifies the environmental variables that explain the wide variability in SOC responses to fire across drylands as a whole. This information can inform management choices that implement nature-based climate solutions 4 , which estimate SOC contributes >50% to the total accrual potential in savannah-grassland ecosystems. Current estimates of potential SOC change from nature-based solutions in grasslands focus on adjustments to grazing regimes (optimizing intensities and plant composition, totalling 0.30 PgCO 2 equivalent (PgCO 2 e) yr −1 ) and avoided conversion (0.23 Pg CO 2 e yr −1 ) (refs. 35,36). Our estimated sink in areas with declining burned area is equivalent to ~0.38 Pg CO 2 e yr −1 , with the caveat that we do not account for changes in other greenhouse gas emissions and we focus on both savannahs and grasslands (for example, CH 4 and N 2 O, although it is unlikely these change much in dry savannah-grasslands). The impact of fire on SOC was considered to be variable enough to be omitted from past estimates 37 . Given that we have identified the environmental conditions that explain a large portion of such variance, we propose fire management should now be integrated into estimates of nature-based climate solutions in savannah-grasslands. Article https://doi.org/10.1038/s41558-023-01800-7 We demonstrate the relative SOC sequestration when fire is excluded is roughly double in sites with high precipitation seasonality and aridity (we used first and third quartiles of the covariates to compare 'low' and 'high' values). SOC sequestration comes with a potential cost of tree encroachment, which can reduce biodiversity-a key consideration that may be offset through management of browsing herbivores, which we did not consider here. For example, browsers-or lack thereof-can determine decadal trends in tree cover in savannahs 38,39 , such that the restoration of their populations may help abate deleterious effects of woody encroachment caused by fire suppression. The carbon-biodiversity tradeoffs are difficult to ascertain because little work has been done assessing cross-site patterns of the sensitivity of dryland biodiversity to changes in fire, and even less on how that relates to SOC, across environmental conditions. Such assessments will be key given that fire has been shown to both decrease 40 , increase 41 or not change 42 biodiversity. Thus, further work is needed to understand the relationship between SOC and biodiversity to manage fire for nature-based solutions. Along those lines, our dataset has large gaps in Europe and Asia, which have unique biogeographic assemblages of plants and differences in fire behaviour, which should be addressed in future studies to test the generality of our findings. However, our dataset does cover the bioclimatic conditions and soil properties representative of the area over which we are making predictions.
Although statistical upscaling provides a benchmark for evaluating the spatial distribution of SOC changes and the general order of magnitude, improving the process-based models should be a priority. Paths forward include simulating potential direct effects of fire on SOC, which may be minimal in some systems with limited organic horizons but become more important in forests. Second is accurately representing tree growth, overestimates of which may inflate resilience of trees in dry areas relative to the greater reductions we observe empirically 16 . Third is capturing the potential for tree-grass coexistence and changes in tree cover within savannah, which is key for how fire frequency changes SOC in a landscape that remains savannah 17,19 . Given ecosystem models are used to estimate both historical and future SOC changes under altered land use 43 , ensuring they capture observed changes arising from fire (a key process used to manage agricultural lands 12,44 ) is critical.
We did not focus on systems with either intense crown wildfires such as boreal spruce forest or smouldering ground fires in peatlands. Although these ecosystems can store large amounts of SOC, they burn less frequently than savannah-grasslands (albeit with increasing frequency) 4 , and the factors that determine direct combustion of SOC are well characterized 4 . Meanwhile, savannah-grasslands have lower per area SOC stocks, but given their large spatial extent and frequent burning, SOC stocks underlying savannah-grasslands that burn equate to ~14.5 PgC (ref. 4). Our sink estimate is small relative to the residual land C sink (3.1 ± 0.6 PgC yr −1 ) (ref. 45), but given we observed continued declines in SOC with the long duration of experiments in our dataset (~60 years), SOC in drylands probably have a relatively long-lasting capacity to sequester carbon.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-023-01800-7.

Study compilation and overview
Our methods were similar to previous meta-analyses of how fire affects soils by calculating a response ratio of soil carbon (C) in plots with different fire treatments within sites and then comparing responses across sites 10 . To obtain data from the literature, we searched for studies that measured the response of mineral soils to repeated burning (completed in May 2020) on Google Scholar. We included 'decadal', 'fire frequency', 'soils', 'repeated burning', 'carbon', 'long-term' to isolate papers with long-term repeated burning manipulations; searches were conducted in English and yielded 156 articles. We used a threshold of plots experiencing at least two fires that had been running for a decade or longer. We made one exception for a site that had been running for nine years because it was in an ecosystem that was not well represented in the dataset. We focused on the mineral soil because of our focus on understanding fire effects in drylands, where it dominates soil C storage pools. We analysed data from the uppermost soil layers ( < 20 cm depth) because these are the most biologically active and likely to be the most responsive to burning. There is evidence that fire can alter soil C > 20 cm deep in some ecosystems 8,10,46 , suggesting our estimates could be conservative.
The majority of our sites were from fire-manipulation experiments where fire treatments were prescribed (43 out of 53). The fire frequencies ranged from one fire every 17 years to one fire every year. The low fire treatments were usually complete fire exclusion. In some cases the authors note an incidental fire that burned through one of the fire exclusion plots, but these fires were mostly rare. One site was running for only 9 years 47 , but we included this study because it came from a grassland that had received >5 fires and was one of the only sites in an arid region. Fire treatments were replicated at the landscape scale in all but 16 sites. Independent replicates of the fire treatments were generally defined based on the application of different fires (for example, fire breaks separated the plots and the managers burned each plot separately).
In several sites, the plots had experienced other types of land use before the establishment of the experiment. For example, many of the forested sites in the southeastern United States were established on abandoned agricultural land or tree plantations. Other types of disturbance also occurred during the experiment such as intermittent drought and herbivory. We assume that there were no overarching biases in the land use history of the experiment that would drive our trends. Herbivory was clearly biased, however, with the African savannah sites exposed to browsing and grazing. We assume that these characteristics are important natural processes in the ecosystem ( just as intermittent drought in many ecosystems) and thus included the sites in our analysis.
In addition to the literature search, we incorporated data from our own surveys of seven sites. Several of these sites replace previous measurements because we wanted to extend the timescale over which fire had been manipulated (for example, in some cases plots had been surveyed in the 1990s and we re-surveyed them).

Different sites and calculations from a previous meta-analysis and field surveys
Matopos studies: weighted the soil C under grasses vs trees by the tree cover. Texture data are in refs. 47,48. Kruger: we used data from our own re-survey across all the plots in 2015. Morton, Hitchiti, Cedar Creek: re-analysed the soil samples in our lab. Cedar Creek duration extended. Limestone and Chimney Springs: surveyed soils with our own sampling to extend the dataset by two decades. University Missouri: added data from Pellegrini et al. 2021 16 . Wet Peachester: now from refs. 49 because this extended the dataset.

Soil texture
We compiled data on soil texture using either (1) measurements within the study site, (2) extrapolations based on discrete soil classification (for example, clay loam, sand and so on; http://www.fao.org/tempref/ FI/CDrom/FAO_Training/FAO_Training/General/x6706e/x6706e06. htm) 50 . When percent texture did not add to 100, we used average clay, then silt and then took the difference with sand. In the cases where only clay was reported, we used the reported value of clay, mean silt in the texture class and then subtracted out sand from 100. In the cases where texture was provided but class was not given, we used texture to assign it to a class.

Climate variables
We used WorldClim 32 to obtain climate variables for each site. In our model selection, we focused on variables related to growing season length and variability and water availability. Specifically, we analysed mean annual temperature, annual temperature range, mean precipitation (mean annual and broken down into the wettest and driest quarters) and precipitation seasonality (coefficient of variation of precipitation calculated with the standard deviation across months within a year divided by the annual precipitation and multiplied by 100).
To integrate water and evaporative demand, we used an Aridity Index 51 . Aridity was calculated from the WorldClim data based on data from 1950-2000. The aridity index is given as: Aridity Index = (mean annual precipitation) / (mean annual potential evapotranspiration) Potential evapotranspiration was calculated via the Hargreaves equation 52 (mm per month) using the WorldClim data as input: We divided the index by 10,000. Aridity classes are calculated from the UNEP World atlas of desertification, 2nd addition. Arid and semi-arid: 0 < AI ≤ 0.5; dry sub-humid: 0.5 < AI ≤ 0.65; humid: 0.65 < AI ≤ 0.75; hyper-humid: AI > 0.75.

Imputed standard deviations
For nine of the sites, we imputed standard deviations. We did this by first calculating the coefficient of variation across all the sites with reported standard deviations and then multiplying the coefficient of variation by the mean for the sites being imputed. Imputing standard deviations using values from studies included within the meta-analysis has been shown to be the most accurate 53 .

Meta-analysis statistics
All statistics were run in R 54 . We calculated log response ratios on soil C concentrations in the high-frequency plots divided by the concentrations in the unburned plots for each site (RR site ).

Nature Climate Change
Article https://doi.org/10.1038/s41558-023-01800-7 To determine variable weights (var), we used the inverse of the variance within a site. We calculated the variance within each site by combining the standard deviation (σ), sample size (n) and sample means (μ) for each treatment within each site. We then calculated the overall effects, 95% confidence intervals and their significance using multivariate meta-analysis models (rma. mv) in the metaphor package 55 .
To evaluate the important predictors of fire effects on soil C, we performed model selection on linear meta-analysis models using the glmulti package 55 with the log response ratio as the predicted variable and environmental variables as the predictor variables. Models were fit using maximum likelihood estimation via the function rma.glmulti. For this analysis we only considered first-order effects and not interactions between environmental variables. Variables included are listed in the 'Weighted variable importance from model selection' in Supplementary Table 2 and bolded are those included in the top model (criterion for top model in the following paragraph).
To evaluate the top model, we first extracted the models with the lowest corrected Akaike Information Criterion within a value of 2. We calculated variable importance by summing the weights of the models that the variable occurred in. Sometimes a cut-off of 0.8 is used to delineate important vs unimportant variables. Here, however, the top model contained several variables that were not above this cut-off, making it less straightforward to evaluate whether a variable was important. Consequently, we tested for variable significance in the top model (which included the variables that had an importance <0.80), which is reported in the main text, using a meta-analysis via linear mixed-effects model (rma) in package metafor with a maximum likelihood estimator. Significance of moderators was evaluated using an omnibus test (Q M test) using a chi-square distribution. Significance of individual coefficients are tested using a standard normal distribution. We allowed for aridity and precipitation to both be in the top model because they had a relatively low Pearson correlation coefficient (σ = 0.34).
We compare 'low' and 'high' values of environmental covariates using the first and third quartiles of the covariates, respectively.

Survey of savannah sites for soil carbon and δ 13 C
In seven of the sites included in the meta-analysis, δ 13 C was measured in combination with soil C (Breaks, Cedar Creek, Hitchiti, IBGE, Satara, MatoposClay and MatoposSand). These sites spanned tropical and temperate regions in North and South America and Africa. The savannahs all contained C4 grasses, which allowed us to use δ 13 C values to partition tree vs grass biomass contributions to soil C (although some contained a mix of C3 and C4 grass). The absolute number of soil samples collected varied across sites because of differences in tree cover, but we used plot-level averages in the analysis. Duration of fire frequency experiments ranged from 25-64 years. The ratio of 13 C to 12 C is assumed to be relatively unchanged by fire compared to the difference between C3 vs C4 photosynthetic pathways 20,34,56,57 .
We used isotopic mixing models to calculate the proportion of soil C derived from grasses vs trees. We performed these calculations in savannah sites where C4 grasses comprised a large proportion (if not all) of herbaceous biomass. We assumed a C4 signature of −15‰ and C3 of −28‰ for the sites where plant isotope values are not available. In other cases, site-specific values were used 19 .
We generated a two end-member mixing model to calculate the relative contribution of trees vs C4 grasses, where x denotes the proportion coming from trees.
We did not necessarily expect the functional form of the curve between soil δ 13 C and total soil C to be linear so we fit generalized additive models with a penalized spline. The maximum degrees of freedom was set to k = 3, where k is a parameter determining the flexibility of the spline. All fitting was done using the mgcv package in R 58 .
We tested the overall effect of fire on δ 13 C values across sites using a mixed-effects model with site as a random effect. Within each site, we tested for significant effects using linear models (Supplementary Table 3).
For this analysis we calculate the relative contributions of C3 trees vs C4 grasses in terms of total stocks of SOC. To calculate this, we use bulk density values reported for the site, standardized it to a depth of 0-20 cm and multiplied total soil mass per area by the soil C concentration.

Field sampling in six fire-manipulation experiments
We sampled fire-manipulation experiments in the Cedar Creek Savanna Fire Experiment in Minnesota, USA, the Reserva Ecologica do IBGE in Brasília, Brazil, and four sites across the Experimental Burn Plots in Kruger National Park, South Africa. These are some of the longest-running fire-manipulation experiments in the world and described in detail in ref. 57. These sites span a gradient from semi-arid savannahs to mesic savannah-forest ecotones with a range in total precipitation and its seasonality.

Experimental design
In all the sites, replicate plots >1 hectare in size have received different prescribed burning frequencies. Burns are conducted at the end of the dry season in the tropical sites and early spring in the temperate sites. Fires are typical of savannahs with low to moderate intensities that result in little adult tree mortality but tend to topkill saplings [59][60][61] .
In all cases the 'low' fire frequency treatment is fire exclusion that was initiated on a savannah landscape. The intermediate frequency approximates the historical average for the region (fire return intervals of 3-5 years in Brazil, 3 years in South Africa, 3 years in the United States) and a high-frequency treatment is greater than the historical mean (fire return intervals of 2 years in Brazil and 1 year in South Africa; the plots in the United States were burned 3 out of every 4 years).

Sampling
We sampled the top 0-20 cm of the mineral soil horizons. Bulk density was measured concurrently and calculated within each fire treatment, thereby controlling for potential compaction under more frequent burning. Soils were dried, sieved to 2 mm and analysed for total C via combustion on Elemental Analysers at Stanford University, Princeton University and the University of Cape Town. Acid digests on soils with high carbonate concentrations were performed to isolate SOC from inorganic C.
Soil cores were taken at 2-5 locations within each replicate plot within each treatment in a site. When necessary, we stratified sampling underneath either a tree canopy or in the grassy matrix. The five cores were then homogenized within each vegetation type. Bulk densities were measured on these samples and used to calculate total stocks. We weighted our calculated soil variables based on the tree cover in a site (for example, 20% tree cover meant that the soil C value under trees contributed only 20% to the plot-level average and the soil C values under grasses contributed the other 80%).
We tested for significant fire effects on SOC, the contribution of C3 plants to SOC and the relationship between SOC and the contribution Article https://doi.org/10.1038/s41558-023-01800-7 of C3 plants using linear regressions followed by an Analysis of Variance (ANOVA).

Fire Model Intercomparison Project simulations
The empirical spatial pattern and climate relationships of fire's impact on soil carbon was compared to that simulated by a set of fire-enabled global vegetation models. Simulated global soil carbon output from seven global fire-vegetation models was obtained from a set of standardized simulations provided by the Fire Model Intercomparison Project (FireMIP 30,62 ): Two sets of simulations were used: first a fully transient simulation with changing climate, population density, land use and [CO 2 ] and another identical sensitivity experiment in which fire was switched off. The difference between both runs indicated the long-term impact fire has on soil carbon as simulated by each model. Climate, land use, population density and [CO 2 ] forcing data was standardized so that inter-model differences can be traced back to differences in model structure and parameterization and not to differences in forcing data. Of the models we used, CLM-Li and LPJ-GUESS-SIMFIRE-BLAZE had nitrogen cycles.
A spin up until the slowest C pool was in equilibrium was performed for each run, during which climate was recycled over the 1901-1920 period and all other forcing data were kept constant at their initial value. After spin up, the models were forced by time varying land use, population density and [CO 2 ] (1700-2013) and climate . Two models slightly deviated from this protocol by starting their transient simulations later (CLM: 1850 and CTEM: 1861), but this should have minimal impact on the results used here. Soil carbon output averaged over the last two decades of each model's simulations was used for analysis. More detailed information about the fire on/off simulations can be found in refs. 31,62. For our statistical analyses, we merged the WorldClim dataset with the model output for direct comparability between the spatial trends in our empirical data. Because we are interested in the factors that determine the impact of fire on soil C, we filtered the data to only include grid cells that had non-zero burned area.
All the DGVMs we used in this study explicitly simulate litter burning. The models do not simulate SOC burning except for CLM-Li. CLM-Li simulates SOC burning only over peatlands. However, for consistency with the other models, we did not consider its impact on SOC in the FireMIP simulations.

Upscaling calculations
We determined the distribution of savannah-grasslands using the World Wildlife Fund (WWF) ecoregions but excluded the Montane Grasslands and Shrublands category because this included steppe (https:// www.worldwildlife.org/biomes/montane-grasslands-and-shrublands). Maps were downloaded on 1 January 2023.
We included the main environmental variables determined in our main meta-regression: percent silt, precipitation seasonality, mean annual temperature, mean annual precipitation and aridity. Climate data were acquired from the same sources as we described above. Percent silt data were taken from the Harmonized World Soils Database v.1.2 (accessed 20 May 2020). We used a 17-year duration of fire treatments based on the 1998-2015 GFED4 records.
We then applied our statistical model to each gridcell across savannah-grasslands. We determined the potential significance of cells using the 95% confidence intervals of the model fit. We present a comparison between global data extrapolated over and data used to fit the statistical model in Supplementary Fig. 1. For our calculations, we restricted our analysis to include only environmental conditions covered by our dataset used to build the statistical model.
Areas that were experiencing either gains or declines in burned area were determined using remote sensing data on fire occurrences from 1998 to 2015 from a previously published dataset 12 . For the declines in burned area, we re-calculated the response ratios in terms of fire exclusion/high frequency to illustrate the potential gain in soil C and transformed them to percent differences; this was repeated for the cases where fire frequency increased. These were then multiplied to the total SOC stocks in each grid cell calculated from the Harmonized World Soils Database v.1.2 63 to transport a percentage to a stock change.
As a comparison between the extrapolation-based SOC flux estimates and empirical measurements, we used the SOC fluxes from the six sites from our field campaign where we measured bulk densities to a standardized depth of 20 cm. Field measurement fluxes were 0.21 ± 0.10 MgC ha −1 yr −1 (mean ± standard deviation) while model-based estimates were 0.37 ± 0.26; the median was 0.31 MgC ha −1 yr −1 (model-based estimates follow a log normal distribution).

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
Data are freely available on Figshare at https://doi.org/10.6084/ m9.figshare.23899530. Source data are provided with this paper.
Article https://doi.org/10.1038/s41558-023-01800-7 Extended Data Fig. 1 | Distribution of semi-arid regions and savannahgrasslands. Distribution of semi-arid regions and savannah-grasslands in yellow with points identifying the locations of our sites in the meta-analysis. We determined the distribution of savannah-grasslands using the WWF ecoregions but excluded the Montane Grasslands & Shrublands category because this included steppe (https://www.worldwildlife.org/biomes/montane-grasslandsand-shrublands). Maps were downloaded on 1/1/2023. We also added the MODIS map of semidry shrub and savannah 64 . References for the specific experimental sites are as follows 65-99 . Article https://doi.org/10.1038/s41558-023-01800-7 Extended Data Fig. 3 | Fire frequency is important in sites but does not explain the relative changes in soil carbon across sites. Fire frequency is not important at explaining the relative changes in soil carbon (C) across sites but it is important within sites. a) Each point is a site, displaying the percent difference in mineral soil C concentrations in the high frequency vs. unburned plots. b-c) comparisons using the intermediate fire frequency treatments on 34 sites. These 36 sites have intermediate frequencies of burning (generally relative to what is believed to be the historical or 'natural' frequency), allowing us to evaluate the impact of increased burning or decreased burning relative to an intermediate level. b) soil C in the intermediate fire frequency plots vs. fire exclusion plots. c) soil C in the high-frequency plots vs. intermediate fire frequency plots. For both b-c Negative values illustrate lower C concentrations in the higher frequency treatment. Circles in b-c represent the means and error bars are 95% confidence intervals. Sample sizes for b-c: grassland: n=3, woody savannah: n=13, broadleaf forest: n=10, warm needleleaf forest: n=6. between precipitation and potential evapotranspiration. Precip. seasonality = Precipitation seasonality, which is the coefficient of variation of monthly precipitation within a year multiplied by 100. Sources for these values are explained in detail in the SI, but briefly the climate data come from WorldClim and soils data from the Harmonized World Soils Database v.1.2 (accessed 5/20/2020).