Human land occupation regulates the effect of the climate on the burned area of the Brazilian Cerrado

Human activities and climate change are transforming ﬁ re regimes globally. The interaction between these two drivers is poorly understood, yet critical if we aim at predicting how biomes will respond to novel ﬁ re regimes. In the Brazilian Cerrado, altered ﬁ re regimes are threatening its unique biodiversity and ecosystem functioning. Here, using geospatial data for the period 1985-2020 and a causal inference framework to design Bayesian statistical models, we demonstrate that a larger human presence in the landscape ( ≥ 40% land-use area) reduces the Cerrado ’ s burned area and hinders its responsiveness to climate; while climatic effects only become apparent in landscapes with little human presence, where hotter and drier conditions increase burned area. Finally, we ﬁ nd spatially heterogeneous burned area trends over time, with increases associated to climate change in landscapes that have remained mostly intact, and decreases caused by anthropic expansion. Both diverging trends have important implications for the conservation of the Cerrado as land-use expansion and climate change continue to unfold.

https://doi.org/10.1038/s43247-024-01521-5 Human land occupation regulates the effect of the climate on the burned area of the Brazilian Cerrado

Check for updates
Carlota Segura-Garcia1 , David Bauman 1,2 , Vera L. S. Arruda3 , Ane A. C. Alencar 3 & Imma Oliveras Menor 1,2 Human activities and climate change are transforming fire regimes globally.The interaction between these two drivers is poorly understood, yet critical if we aim at predicting how biomes will respond to novel fire regimes.In the Brazilian Cerrado, altered fire regimes are threatening its unique biodiversity and ecosystem functioning.Here, using geospatial data for the period 1985-2020 and a causal inference framework to design Bayesian statistical models, we demonstrate that a larger human presence in the landscape (≥40% land-use area) reduces the Cerrado's burned area and hinders its responsiveness to climate; while climatic effects only become apparent in landscapes with little human presence, where hotter and drier conditions increase burned area.Finally, we find spatially heterogeneous burned area trends over time, with increases associated to climate change in landscapes that have remained mostly intact, and decreases caused by anthropic expansion.Both diverging trends have important implications for the conservation of the Cerrado as land-use expansion and climate change continue to unfold.
Fire is a complex process that involves the interplay between anthropic activities, the climate, and vegetation structure that both shape and are in turn affected by this phenomenon 1 .Understanding the drivers of fire and how they interact is essential to tackling the many challenges that fire poses globally, from human health and the economy 2,3 to its influences on the Earth's biogeochemical cycles 1,4 , and the conservation of both fireadapted and fire-sensitive ecosystems 5 .Indeed, fire is a natural disturbance process playing a crucial role in shaping ecosystems 6 , but human activities are altering the Earth's fire regimes 7,8 .Over the last two decades, global burned area has declined in association to agricultural expansion and intensification, with greater declines concentrated in savannas 9 .In other regions, modern human activities have escalated the occurrence of fire 10 indicating that anthropic effects vary depending on the context 9 .Simultaneously, climate change is intensifying the risk of fire weather conditions 11,12 and the frequency of extreme weather events 13 , exacerbating fire activity and burned area.Therefore, globally, the two main drivers of fire, humans and the climate, could be pushing burned area in opposite directions, but there is still considerable uncertainty in how they interact 14,15 .Quantifying such interactions is fundamental to improve our understanding of feedbacks between the climate, fires and vegetation 4,16 .
Savanna ecosystems stand out as a paradigmatic case of these apparently contrasting processes 17 .Savannas are a predominant biome covering around 20% of the terrestrial surface with extensive presence of humans.These ecosystems play a relevant role on the Earth's carbon cycle contributing to around 25% of the total gross primary productivity 18 and have an -often neglected -high biodiversity and ecological value [19][20][21] .Savannas are composed of fire-adapted vegetation that depends on this natural disturbance for their continuity and the preservation of their biodiversity 22 .Therefore, departures from the historical fire regimes in either direction can have important ecological consequences potentially resulting in ecosystem shifts to other forms 23,24 , as well as effects on the Earth's carbon cycle that can feed back to the climate 4 .In these regions, the observed declines in burned area 17,25 are often linked to landscape fragmentation 17,26 , population density, certain agricultural practices or increasing grazing pressure that disrupt fuel continuity impairing fire spread.Active fire suppression is also a relevant cause of burned area reduction, known to lead to woody encroachment causing biodiversity losses 27,28 and increasing the risk of more severe wildfire events 29 .Regarding the climate, some studies forecast increases in savannas' burned area associated with warming and the lengthening of the dry season 16 .However, it has also been suggested that human controls on savannas' burned area may limit its responsiveness to climatic conditions 30,31 , but it is still unclear how much human activities regulate the climate effects on fire.
Here, we focus on the Cerrado, a Brazilian ecoregion comprising several savanna ecosystems and considered to be a biodiversity hotspot 32 .The Cerrado occupies over 2M km 2 , with around 50% of its area already converted to anthropic land uses 33 , and undergoing a warming and drying process 34 associated to climate change 35 .These savannas present spatially heterogeneous patterns of burned area 36 in accordance with its mosaic vegetation and widespread agricultural areas, with both increasing and decreasing burned area trends 25,26,37 .Much research has been devoted to studying climate effects on the intra-annual [38][39][40] and inter-annual 35,38 variability of burned area, identifying increasing burned area trends in relation to climate change 25,41 .Regarding anthropic area and landscape fragmentation, these seem to have a more complex association with fire activity 25,26,30 , while the burned area is concentrated in regions with larger natural remnants 36,39 .Disentangling the effects of climate change and human presence is essential to devise effective fire management policies and conservation plans that are targeted to the characteristics of the landscape and its vegetation composition.
Using a causal inference framework and a Bayesian approach, we aim to study whether and how the level of human presence in the landscape regulates the effects of the climate on the burned area of the Cerrado.We use the percentage of anthropic area -comprising agricultural, urban and mining land uses -as a proxy for the level of human presence in the landscape to, first, determine its effects on burned natural area -measured as the percentage area of Cerrado native vegetation that burned.Then, we quantify the interaction between the climate and human presence by modelling the effects of the climate on burned area in landscapes with different levels of human presence using Bayesian multilevel regression models.Finally, we explore the spatial distribution of the temporal trends in burned area, climate and anthropic expansion and build a model to determine the causal effects of the last two as predictors of the former.

Results
Limited burned natural area in human-dominated landscapes First, we modelled the causal effect of anthropic area percentage on the burned natural area percentage of the Cerrado.With this aim, we combined annual land use and burned area data on a 0.2°grid over the Cerrado for the period 1985 to 2020 (see Methods).As we were interested in studying the causal effect of the predictor -anthropic area -, we used a structural causal modelling framework 42 which, based on assumptions and a priori knowledge about our study system, provides a set of rules -known as the backdoor criterion 42,43 -to design models with the sufficient statistical adjustments to limit bias and confounding, thus enabling causal inference.
Therefore, we first designed a Directed Acyclic Graph (DAG), DAG1 (Supplementary Fig. 1), a form of heuristic causal model used to visually represent the causal relations between the different variables of our system and inform covariate selection (see Methods and Supplementary Methods).This DAG was designed based on author's expert knowledge about savanna fire and its drivers, the literature and exploratory data analysis, and reflects our assumptions about how the different variables affect each other (see Methods and Supplementary Methods).Hence, based on DAG1, we built model M1 (Equation ( 3)).This model included anthropic area percentage as the explanatory variable of interest for which we wanted to determine its total causal effects on burned natural area, along with protected area and latitude as control variables.These last two covariates were needed according to the backdoor criterion to enable the causal interpretation of anthropic area effects, but it is important to stress that their own modelled effects cannot be interpreted causally (see Methods).
With model M1, we found that the higher the presence of humans in the landscape in terms of area occupied, the smaller the fraction of burned natural area (Supplementary Fig. 2, slope coefficient: −0.130, [−0.135, −0.126], 95% Highest Posterior Density Interval, HPDI, % burned natural area % anthropic area −1 ), implying that human dominance of a landscape acts to limit the fraction of burned area in Cerrado's vegetation as hypothesised.
Finally, we calculated the R-squared coefficient adjusted for Bayesian inference 44 of model M1 (0.544 [0.542, 0.547] 95% confidence interval, Supplementary Table 1).However, it is important to bear in mind that the aim of this study was to determine total causal effects, not to maximise model predictability.Therefore, in this and subsequent models, the R-squared coefficient was calculated just as an indication of the proportion of variance explained by the model, but it should not be interpreted as a measure of the causal relation between the predictors and the response variable.
Climate effects on burned natural areas are conditional on human landscape dominance Then, we assessed the causal effects of climatic controls on the burned natural area of the Cerrado, and the extent to which these controls depend on anthropic factors such as anthropic area.Thus, we quantified the causal effects of temperature, vapour pressure deficit (VPD), same and previous year total precipitation and dry season precipitation (Supplementary Fig. 3) on the observed percentage of burned natural area (Supplementary Fig. 4) using climate data from different sources (see Methods).
Using the same DAG1 as before, we built two Bayesian multilevel models (M2 and M3, Eqs. ( 4) and (5), see "Methods" section) where we allowed the climate effect coefficients to vary by decile of anthropic area percentage (Supplementary Fig. 5) to determine the interaction of the climate covariates with the level of human occupation (Fig. 1, see R 2 coefficients in Supplementary Table 2).Using the backdoor criterion, we included latitude as a control variable to avoid confounding in model M2 for temperature and VPD (see "Methods" section), while no covariates were needed in model M3 for the different precipitation variables.
We found that the slope coefficients for all climate predictors showed a non-linear change in behaviour with the level of human presence in the landscape (Fig. 1).In cells with a large area of natural vegetationaround 40% or less of anthropic area -, we observed the expected effect of the climate on burned natural area: increases in temperature and VPD cause increases in burned natural area; the lower the dry season precipitation, the larger the burned area; and higher precipitation in the previous year brought about larger burned area in the following year.We found no strong causal effect of same-year total precipitation (Supplementary Fig. 6a).In contrast, where anthropic land uses dominate the landscape -~40% or more of the cell's area -, we observed no effect of these climatic variables on the percentage of burned natural area, suggesting strong controls of land use activities.
Those observations in the [90-100]% decile of anthropic area showed no effect of any climatic variables on burned natural area (Fig. 1).We used these results as a reference to test the presence of thresholds from which the effect of the climate in following deciles was different compared to this extreme case.With this aim, for each climate covariate, we generated pairwise contrasts by subtracting the posterior distribution of the slope coefficient obtained for the [90, 100] % decile from the corresponding posterior of every other decile (Supplementary Fig. 7).We considered that if the 90% HPDI of the pairwise contrast (the difference between the posteriors) did not encompass 0, then there was an effect of the climate on the burned natural area in that decile compared to the no-effect [90, 100] % reference.Indeed, we found that below some anthropic area percentage, temperature (50% anthropic area), VPD (40%), previous year precipitation (40%) and dry season precipitation (20%) had an effect on burned natural area.Hence, for each climatic variable there was a certain threshold value of anthropic area percentage at which burned natural area transitioned from being not responsive to becoming increasingly affected by climate.This transition was sudden, that is, non-linear.
Among these climatic variables, VPD and temperature seemed to have a more marked effect on burned area than dry season and previous year precipitation (Fig. 1), suggesting the former two may emerge as stronger climatic controls of burned natural area.We explored this further by comparing the predicted effect on burned natural area of each climate variable at their average value to their predicted effect at their average value plus (minus in the case of dry season precipitation) 2 standard deviations in each anthropic area decile (Supplementary Fig. 8).We observed that, in the less anthropic landscapes (particularly the [0, 10)% decile, see Supplementary Fig. 9), when increasing (decreasing) each climatic variable's value by 2 standard deviations, VPD was the variable with the largest increase in burned natural area, followed by temperature; while the change in effect of dry season and previous year precipitation was more limited.

Trends in burned natural area, climate change and anthropic expansion
We explored the spatial distribution of the temporal trends in burned natural area, the climate and anthropic area (Fig. 2) using a Bayesian linear regression model with year as the independent variable (model M4, Equation ( 6)).We observed a heterogeneous spatial structure of increases and decreases in burned natural area percentage over the period 1985-2020 (Fig. 2a).We also detected a widespread expansion of human land uses (Fig. 2b) particularly acute in some regions that were mostly natural at the beginning of the study period (Supplementary Fig. 5a).Regarding the climate, there was a generalised increase in temperature in the Cerrado (Fig. 2c), but we found no evidence of clear trends in precipitation, in agreement with previous studies 34 .
Then, we studied whether these trends in temperature (tT) and in anthropic area (tA) had a causal effect on the temporal changes in burned natural area (tBA), and whether these variables mediated each other's effects.With this aim, we built a second DAG (DAG2, Supplementary Fig. 10) representing the relations between these different trends in time.Since the extent of anthropic area expansion was contingent on the level of human occupation at the beginning of the study period, we grouped the cells in five intervals of 20% anthropic area in 1985.Hence, we modelled the causal effects of tT and tA allowing the slope coefficients to vary by interval using a Bayesian multilevel model (model M5, Eq. ( 7)) based on DAG2.In addition, we included an interaction term between tT and tA to allow the effect of climate change severity to depend on the level of change in anthropic area.Finally, we explored the implications of model M5 looking at the implied posterior predictions of tT and tA on tBA in the study period (Fig. 3 and Supplementary Figs.11 and 12).That is, we explored the values of tBA predicted by model M5 for different combinations of tA and tT.We found that, indeed, the effects on tBA exerted by tT and tA varied depending on the pre-existing levels of anthropic area in 1985 (Supplementary Figs.11 and 12, see R 2 coefficients in Supplementary Table 3).For those cells that were mostly anthropic in 1985 (>40% of anthropic area), we found that both predictors had limited to no effects on tBA.However, for those cells that were mostly natural at the beginning of the period (<40% of anthropic area), model M5 predicted a different effect for each predictor.These differing effects were even more evident in cells with <20% of anthropic area in 1985 (Fig. 3).
In these cells, the effects of the rate of temperature increase seemed to be mediated by the extent of anthropic expansion (Fig. 3a).Where human presence in the form of land uses markedly expanded (tA value of 1.36% per year, or 48.96% in the period, red line in Fig. 3a), burned natural area decreased over time (negative tBA) regardless of tT.In contrast, where anthropic area remained low and roughly at the same level (tA of 0%, blue line in Fig. 3a), we found that warming caused increases in burned natural area over time (positive tBA), and the more acute the temperature increase, the greater the increase in burned natural area over time.
This behaviour became even clearer when we fixed tT at a large (0.039 °C year −1 , or 1.40 °C in 36 years, red line in Fig. 3b) and small (0.027 °C year −1 , or 0.97 °C in 36 years, blue line in Fig. 3b) value and looked at the effect on tBA of varying tA (Fig. 3b).In both cases, when tA is small (limited anthropic expansion), burned natural area increases over time (positive tBA), and the larger the increase in temperature (tT), the larger the tBA (larger tBA value in the red line than in the blue line in Fig. 3b).However, when tA is large (large anthropic expansion) tBA becomes negative (burned area decreases over time) independently of tT value (similar tBA values in the red and blue predicted values).

Discussion
This study demonstrates that anthropic land use activities limit the responsiveness of burned natural area to climatic conditions in areas of the Brazilian savannas where humans have a substantial presence (Fig. 1).Moreover, we find that the interaction between these two drivers is nonlinear: humans completely limit burned natural area responses to climate up to a certain level of anthropic area percentage beneath which the expected relation between the climate and burned natural area is recovered.These mitigating effects extend to areas where anthropic land uses are not the major land cover (up to 30-40% of the area, Supplementary Fig. 7), indicating that human controls on fire stretch to landscapes where considerable natural vegetation remains.In addition, we find that a larger anthropic area is associated with an overall reduced percentage of burned area in natural land cover areas (Supplementary Fig. 2).
We use anthropic area percentage as an indicator of the presence of human activities and their influence in the landscape, as it encompasses a series of direct anthropic effects on burned natural area.Most anthropic area in the Cerrado consists of agricultural land 33 conforming fragmented landscapes 26 , which limits fire spread through fuel discontinuities 14 .Another cause of burned area limitation is active fire suppression, a common practice in the Cerrado enforced by policies 45 , although there have been substantial efforts to change these approaches [45][46][47] .Among other consequences, this reduction in burned natural area contributes to woody encroachment 27,28 , in turn leading to the loss of biodiversity 48 and potential ecosystem shifts 24 .
On the other hand, where anthropic area is much smaller, landscape fragmentation is reduced and, consequently, fuel bed continuity may no longer be a limiting factor to fire spread 14 .Furthermore, access to these areas may be more challenging, reducing fire suppression opportunities.Accordingly, in mostly natural land cover areas fuel load and its flammability may become more relevant determinants of burned area 49 , both of which are influenced by climate.Indeed, our results demonstrate that in these landscapes climate is a driver of burned natural area (Fig. 1).We find that a higher temperature and VPD causes an increase in burned natural area, as expected, since temperature increases the moisture holding capacity of the air, drying the fuel and making it more flammable 50 .In agreement with earlier studies, we observe that lower precipitation in the dry season contributes to larger burned areas as it also influences fuel flammability 30,51 , while higher precipitation in the previous year causes larger burned areas as it affects vegetation productivity and, hence, fuel load 52 .However, although these fires are occurring in mostly non-anthropic land cover areas, they are not necessarily naturally ignited (e.g. by lightning), as humans are still the main ignition source 5 .This fact does not conflict with the low presence of humans in these landscapes because most burned area is produced by few very large fires 53 , so a small number of ignitions can result in large burned areas 14 , especially under hot and dry conditions.Regarding the trends in time, we find a heterogeneous pattern of increases and decreases in burned area (Fig. 2), in line with previous studies 25,36,37 .When exploring the relation between trends in burned natural area, anthropic expansion and climate change, model predictions reinforce our results presented above.First, changes in burned natural area over time are conditional on the extent of anthropic area at the beginning of the period.Landscapes that were mostly anthropic in 1985 do not present strong trends because burned natural area was already limited and not responsive to climate (Supplementary Figs.11 and 12).As such, significant trends in burned natural area are concentrated in cells that were mostly natural in 1985.We find that, in these areas, trends are firstly subject to the extent of anthropic area increases (Fig. 3), which brings about a reduction in burned natural area -adding to 31.1% of the cells -, probably in relation to increasing fragmentation and human controls on fire.In contrast, those cells that have remained more intact are the ones where burned natural area increased over time, totalling 9.7% of the cells.Our model indicates that these increases are associated to the warming and drying of the Cerrado due to climate change.
Indeed, both VPD and temperature, which we identify as the strongest climatic drivers of burned natural area change out of those considered (Supplementary Figs. 8 and 9), have increased over the Cerrado in the past few decades and are expected to continue rising 34 .Such increases have been particularly acute towards the end of the dry season 34 when climatic conditions are most favourable for fire spread and when most human-started fires occur 36 .Our findings indicate that as the warming and drying conditions continue lengthening and exacerbating under climate change, we can expect these increasing trends in the burned areas of the more intact regions of the Cerrado to continue.Therefore, we identify those areas with larger natural vegetation remnants as being more vulnerable to climate change.Among these, the region of the Parque Estadual do Araguaia in the centrewest of the Cerrado stands out for its considerable increase in burned natural area over the period.This area is formed by seasonally flooded savannas particularly vulnerable to droughts 54 , suggesting that the type of vegetation may be important in the evolution of fire activity.Although savannas are fire-adapted ecosystems, changes in fire regimes towards larger, more frequent, and intense events under climate change may have degradation effects on these ecosystems, which could hamper their recovery following disturbances in the future 55 .In this sense, practices such as integrated fire management are proving to be effective strategies to prevent large and severe fires while contributing to preserve the Cerrado's vegetation mosaic 29,56,57 .
In conclusion, climate change and human land use activities in the Cerrado are altering the fire regimes of the region, both with important ecological consequences.Importantly, we show that human controls over the Cerrado's fire dynamics are substantial and act to limit burned area and to regulate the effect of climate and climate change, contributing to the growing recognition of the role of humans on shaping fire patterns 9,14,30,58 and emphasising the need to improve their representation in projections of fire activity and in climate-vegetation models 4,16 .We find that the greatest changes in burned natural areas over time occurred in the most intact areas, which experienced increases in burned area due to climate change, and in areas with more severe agricultural expansion, which brought about decreases in the proportion of burned natural area.Both directions of burned area change can have important deleterious effects on Cerrado vegetation.Future research could build upon these results by exploring whether these responses to human activities and the climate differ between vegetation formations or land use type -mainly croplands or pastures 26 .The widespread changes of fire activity in either direction, in combination with other consequences of large-scale agricultural expansion and climate change, highlight the multiple pressures the Brazilian savannas are subject to.Considering the importance of the Cerrado in providing local livelihoods, its major role in providing and distributing regional water resources 59 , and its rich biodiversity and high levels of endemism 32 , our findings stress the urgent need to extend conservation practices beyond protected areas and to deploy appropriate fire management policies 45,46 .

Study area
In this study, we focus on the Cerrado, a region of over 2M km 2 comprising the Brazilian savannas.These savannas are very biodiverse 32 and spatially heterogeneous, with vegetation ranging from grasslands to tracts of gallery forests and other woody forms 60 .The climate of the Cerrado is characterised by a strong seasonality, with a wet and a dry season, when most burned area occurs 36 .The average temperatures range between 18°C and 28°C increasing about a latitudinal gradient towards the Equator (Supplementary Fig. 3a).Precipitation varies from about 800 mm to 2000 mm of total annual precipitation, roughly increasing from East to West towards the Amazon (Supplementary Fig. 3b).
Over the last four decades, the Cerrado has experienced a rapid agricultural expansion affecting around 45% of its area 61 .In this study, we focus on the 36-year period from 1985 to 2020.It is important to note that already in 1985 large areas in the south of the Cerrado had been converted to pastures and croplands (Supplementary Fig. 5a).
As we use data from different sources with various spatial resolutions, we divide the Cerrado into a 0.2°grid (~22.2 km 2 ) and summarise the different variables at the cell level (see details below).In the analysis, we only include cells with more than 90% of their area within the Cerrado.

Land use data
We use annual maps of land use and land cover from MapBiomas (Collection 6.0, 33 ) for the period 1985 to 2020.This open data and open source project used a random forest algorithm to classify pixels from annual mosaics of Landsat images at 30 m spatial resolution into several categories.This dataset was developed using Google Earth Engine and is available through the MapBiomas toolkit on the same platform.
In this study, we first simplify the MapBiomas categorisation to either natural land cover pixels (herein natural pixels), anthropic land use pixels (anthropic pixels), or other pixels based on the legend codes of Collection 6.0 33 .Under natural land cover, we group all subcategories under "Forest", those of "Non Forest Natural Formation", and "Beach, Dune and Sand Spot".As anthropic land uses, we group all subcategories under "Farming", along with "Urban Area" and "Mining".The remaining categories -"Other non Vegetated Area", "Water" and "Non Observed" -are grouped as "other".
For each cell and year, we calculate the percentage of anthropic area as where N stands for number.Therefore, %anthropic area = 100−% natural area, as we are leaving out of the calculation the "other" categories because they are not susceptible to burning.

Burned area data
We obtain the burned area data from MapBiomas Fogo (Collection 1.0, ref. 62) for the study period, 1985-2020.This dataset consists of maps where pixels are identified as either burned or not burned, along with the month in which the pixel was mapped as such.These maps were created from mosaics of Landsat images at 30 m spatial resolution with a 16-day interval.
For each cell and year, we use the number of natural pixels classified as burned to calculate our main variable of interest, the percentage of natural area that burned (herein burned natural area), % burned natural area ¼ 100 Á N burned natural pixels N natural pixels : ð2Þ

Climate reanalysis data
We use various climatic variables from different data products because of their relation to burned area.We use the 2 metre temperature parameter from ERA5-Land Monthly Averaged climate reanalysis data at 9 km spatial resolution 63 .We use this variable along with the 2 metre dewpoint temperature from the same dataset to calculate vapour pressure deficit (VPD), a direct measure of the atmospheric demand for water.We calculate VPD as the difference between saturation pressure (e s ) and actual vapour pressure (e a ) following the procedure detailed in 64 .Finally, we calculate the annual average temperature and VPD for each cell and year.We use precipitation data from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS, version 2.0), a quasi-global rainfall dataset with 0.05°spatial resolution covering more than 35 years 65 .This data was acquired in the Pentad version (5-day precipitation accumulation in mm).To obtain the total precipitation per year and cell, we first calculate the monthly total precipitation per pixel, then we calculate the monthly precipitation per cell taking the spatial average, finally adding up the monthly averages to the annual precipitation.
For each calendar year, we calculate the total dry season precipitation as the sum of the monthly precipitation of the dry season months (usually around May to September).We identify the dry season months for each year as those in which precipitation is smaller than the potential evapotranspiration (PET).We obtain the PET from TerraClimate 66 , a dataset of monthly climate data at the global scale and 0.5°resolution, and calculate the cell's spatial average PET.
All climate data from the different products were downloaded through Google Earth Engine 67 .

Protected area data
In the first model, we include the percentage of protected area in a cell and year as covariate.We obtained the data on the conservation units and indigenous lands of Brazil from the Ministério de Meio Ambiente do Brasil (MMA) 68 .We calculated the percentage area of each cell and year occupied by those protected areas established on or before the year.

Data analysis
To tackle the different research questions, we use various Bayesian linear regression models and Bayesian multilevel models (M1 through to M5), whose structure and variable inclusions were guided to allow causal inference.With this aim, we first draw two Directed Acyclic Graphs (DAGs) describing our system (DAG1 and DAG2, Supplementary Figs. 1 and 10).Based on these DAGs, and depending on the causal effects of interest, we applied the backdoor criterion -a test that allows to determine whether and how we can compute a causal effect -to inform our modelling strategy and the selection of covariates to explore causal implications 42,43,[69][70][71] .
Creation of a causal inference framework.Directed Acyclic Graphs (DAGs) are graphical representations of the causal relations between the different variables of interest, along with other variables that may bias or alter the causal association between a given explanatory variable of interest and the response variables 42,43 .We build DAG1 (Supplementary Fig. 1) to represent the relations between the variables that may affect the percentage of burned natural area, and DAG2 (Supplementary Fig. 10) to depict the relations between the different variables that may affect trends in burned natural area with time.The design of both DAGs is based on the authors' expertise, on supporting literature, and exploratory data analysis of the variables that may play a role in the system and their relations.
Fire, and burned area in specific, is a complex phenomenon that involves many different drivers and actors interacting in an intricate manner, and while fully describing complex multifactorial systems such as the one studied here can be challenging, we explicitly acknowledged and accounted for what we consider to be all the major nodes and causal links relevant to the questions of the present work (details in Supplementary Methods).These DAGs therefore provided a rational basis upon which we conditioned causal inference using statistical models, while making the ecological assumptions underlying our models and causal inferences transparent.DAGs entail testable implications of relations between variables, like conditional independencies 69,70 , that we have checked based on the associations and conditional associations in the data to validate the DAG structure, wherever possible.See Supplementary Methods for an in-depth explanation of the DAGs designed.
The main interest of this study is to explore causal effects of different explanatory variables of interest to help understand the specific mechanisms by which human land occupation and the climate affect burned area, and how they interact.Thus, each model we build (see following sections) is aimed at quantifying the total causal effects of a variable of interest, notwithstanding the predictive power of the model.For each explanatory variable of interest (anthropic area percentage, climatic variables, or trends in these variables), the corresponding DAG shows the direct and indirect causal paths between this causal variable and the response variable (burned natural area percentage, or its trends), and helps identify the covariates that we need to control for to avoid confounding: that is, those variables that are necessary to close the non-causal paths between the explanatory and response variables that could introduce bias.These control variables (protected area and latitude in model M1, and latitude in model M2, see following sections) are identified using the backdoor criterion.Furthermore, this criterion also indicates those variables that must not be included to avoid opening non-causal paths.As the non-causal paths may be different for each variable of interest, the backdoor criterion may indicate different choices of covariates for each causal question and, thus, imply building different models.It is important to emphasise that, as a consequence of building our models in this manner, only the effects of the variables of interest can be interpreted causally, while the model outcomes for the control covariates have no causal interpretation as they could be confounded -and thus their specific quantification is not relevant to our research questions.
Finally, we would like to remark that, as a consequence of the explained above, our models include only those variables needed to address our research questions about causal relations.Therefore, there may be of course other variables that also affect burned area causally, or that are good predictors of this quantity, but that we have not included in our models because they are neither the variables of interest to the research questions nor needed for control.In addition, because the models are not designed to maximise their predictability, their overall explanatory power may be limited.Instead, the results from this work should be assessed in the light of the uncertainty of the model parameters, which, as we are using a Bayesian approach, stems from the posterior probability densities -which we summarise using the High Posterior Density Interval, the shortest interval in a posterior probability representing a certain confidence interval.
Effect of anthropic area on burned natural area.Model M1 is a Bayesian linear regression model with anthropic area percentage as the predictive variable of interest and the percentage of burned natural area as the response variable.We run this model for all data points, consisting of each grid cell per year.We construct the model based on DAG1 (Supplementary Fig. 1).In our system, the causal effect of anthropic area percentage on natural area burned may be confounded by the extent of protected area and by the latitude (see Supplementary Methods for further explanation), since both are common causes of the explanatory and the outcome variables (backdoor paths).Therefore, we include these two variables as controls to block the non-causal paths, so that the slope coefficient of the anthropic area covariate will be the total causal effect of this variable on burned natural area (conditional on our DAG) 69,71 .We therefore model burned natural area as: where BA represents percentage of burned natural area, α is the expected average BA when all covariates are zero (intercept), β, η, and ψ represent the burned natural area response to anthropic area percentage (A), latitude (L) and percentage of protected area (PA).It is worth stressing here that in this model only β is to be interpreted causally, while η, and ψ are only there to block non-causal paths (see above).We use parameters γ cell and θ t as varying intercepts to capture the residual variation in expected burned natural area among grid cells and years respectively.Full model specification is in Supplementary Methods.
Climate effects on burned natural area and mediation by anthropic area percentage.In models M2 and M3, we use data for each cell and year observations to test the effect of the climate on the percentage of burned natural area for different levels of anthropic area percentage.However, burned area is a stochastic process that usually follows a distribution within the family of power-law distributions 53 : small burned areas are considerably more frequent than extreme events (larger burned areas), which occur only rarely.In our data, we observe this structure for each interval of each climatic variable (Supplementary Fig. 13).We notice that it is the likelihood of the larger burned areas -the tail of the power law -that changes with the climate (Supplementary Fig. 13).Therefore, we build these two models using a subset of the data comprising the larger events per climatic interval.That is, we divide the values of a certain climatic variable into regular intervals, and for each one we subset those data points corresponding to the upper 95th percentile of burned natural area.Additionally, we observe that the cell-year observations experiencing the highest temperatures in the study period (1985-2020) do not follow the same behaviour than the rest of the data points (Supplementary Fig. 14a): all the burned natural area events for these temperatures are far smaller than the events with temperature immediately lower.When mapping the number of years each cell has experienced an average temperature equal or larger than 28.5 °C, we find that these events occurred in a small number of cells, which are clustered together and located at the border with the Caatinga biome (Supplementary Fig. 14b).This difference in behaviour could hence be related to this area's vegetation being more similar to that of the Caatinga (xeric shrubland and thorn forest), with distinct fuel types and loads.We briefly explored the climate of these cells and found that they also present some of the more extreme maximum climatological water deficit (MCWD) values in the Cerrado (Supplementary Fig. 14c), but we do not investigate this matter further as it is beyond the scope of the study.Since these cells are clustered together and jointly present a markedly different behaviour, we excluded them from this part of the analysis.
Based on DAG1 (Supplementary Fig. 1), we use two different models to quantify the causal relations of the climate on the burned natural area.When studying the effect of temperature and VPD, our model must include latitude as a predictor (model M2).This is because latitude is a common cause of temperature (and therefore VPD) and burned natural area, hence expected to confound the causal effect of temperature on burned natural area if not controlled for.To explore the causal effects of the different precipitation-related variables (same-year total precipitation, same-year dry-season precipitation, and previous-year total precipitation), there is no need to control for any other variable (model M3), as there are no alternative paths going through a common cause between the predictor and the response variable (Supplementary Fig. 1).
Models M2 and M3 have a two-level hierarchical structure where we use a linear regression to predict percentage of the burned natural area (BA) as a function of a climatic variable (C) while allowing both the slope and intercept coefficients to vary by decile of the anthropic area percentage factor, with 10 levels.That is, and In the models, [A] indicates that the corresponding coefficient is allowed to vary by decile of anthropic area, α [A] are the expected average burned natural area for each decile (intercepts), and β [A] corresponds to the response of burned natural area to the each climatic variable in each decile.Parameters γ cell and θ t are varying intercepts used to capture the residual variation in expected burned natural area among cells and years respectively.Lastly, η represents the relation between burned natural area and latitude (L).Again, only the β [A] coefficients are to be interpreted causally.The full model specifications including group-level covariance matrices, prior and hyperprior specifications can be found in the Supplementary Methods.Model M2 was run separately for T and VPD, and model M3 was run separately for same-year total precipitation, same year dry season precipitation, and previous year total precipitation.Finally, we compare the effects of the different climatic variables to explore whether any of them may have a stronger effect on burned natural area than the others.With this aim, we explore the change in burned natural area when increasing (decreasing in the case of dry season precipitation) the value of each climatic variable by 2 standard deviations from their average observed value for each decile of anthropic area.That is, we use the joint posterior of the fitted models to generate burned natural area distributions under a "mean value" and under a "mean value + (or -) 2 standard deviations" of the climatic variables.This results in two predictive posteriors for each anthropic area decile group.We then calculate the "contrasts", that is, we subtract the two group-level posteriors for each variable and decile.In this manner, we can compare the change in burned natural area when varying all climatic variables by the same amount (within each variable's distribution).
Trends in burned natural area, climate change and anthropic expansion.For each cell in the 0.2°grid over the Cerrado, we estimated the annual rate of change (temporal trends) over the period 1985-2020 of burned natural area percentage, anthropic area percentage and various climatic variables.We calculated the trends using Bayesian linear regression, where V generically represents the response variables of interest for a cell i in a certain year, t represents time measured in years, α characterises the average value V of a cell in 1985, and β is the annual rate of change of V -the parameter of interest.We ran this model separately for burned natural area, anthropic area, temperature, VPD, and precipitation using uninformative priors (Supplementary Methods).We consider that a certain cell does not show clear trends in time if its β 95% HPDI encompasses 0. Finally, we map the temporal trends obtained for each variable (Fig. 2) over the Cerrado.
We then proceed with a model (M5) aiming to test whether and how the temporal trends in anthropic area (tA) expansion and in temperature (tT) are responsible for the spatial variation in the temporal change in burned natural area (tBA).These temporal trends are the β parameters obtained in model M4.In this case, we only use those cells that show clear trends in slopes for tA, tT and tBA (95% HPDI).We built model M5 based on DAG2 (Supplementary Fig. 10), which we designed using existing knowledge of similar systems and testing the associations and conditional associations implied by the DAG among the variables to validate its statistical plausibility.In the model, we include a linear interaction term between tT and tA to allow each predictor to regulate the effect of the other.However, the extent of anthropic expansion in the study period is contingent on the percentage of anthropic area in each cell in the year 1985: cells that were mostly intact at the beginning of the period could experience any degree of expansion in anthropic area, great or small, and they consequently show a wide range of tA values.By contrast, cells that were already largely converted in 1985 could only experience small tA values.Hence, we can expect different effects of the predictor variables depending on the anthropic area extent of the cells at the start of the period.Therefore, we classify cells in five intervals of 20% of anthropic area percentage in 1985, and allow the intercept and slope coefficients of tT and tA to vary by group by building M5 as a Bayesian multilevel model, where α ½A 85 characterises the average tBA of cells with a certain level of anthropic area in 1985 [A 85 ], and β 1;½A 85 , β 2;½A 85 , and β 3;½A 85 characterise the response of tBA to tA, tT and their interaction, respectively, for each [A 85 ] group.The full model specification is in the Supplementary Methods.Finally, we explore the implications of model M5 generating predictive posterior draws for each A 85 group by fixing one of the predictors at its 90th and 10th percentile for each A 85 group, and exploring the effects of the other predictor on tBA throughout its observed value range (Fig. 3 and Supplementary Figs.11 and 12).

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

Fig. 2 |Fig. 1 |
Fig. 2 | Spatial distribution of the trends in time (slope coefficients) obtained from model M4 for each cell in a 0.2°grid over the Cerrado using the year as the predictor over the period 1985-2020.Trends for a burned natural area percentage, b anthropic area percentage, and c temperature.Grey indicates cells with no evidence of clear trends -the 95% HPDI encompasses 0.

1 ) 1 )Fig. 3 |
Fig.3| Predicted temporal trends in burned natural area (tBA) from model M5 for those cells that were mostly intact in 1985, with 0-20% of anthropic area (80-100% of natural area).These predictions are generated by fixing one predictor and exploring the effects of varying the other.a The effects of trends in temperature (tT) on tBA for two fixed values of trends in anthropic area (tA) corresponding to the 10th percentile (0.00 % year −1 ) and 90th percentile (1.36% year −1 , or 48.96% in the period) of the range of tA values observed in this group of cells.b The effects of tA after fixing tT at the 10th (0.027 °C year −1 , or 0.97 °C in the period) and 90th (0.039 °C year −1 , or 1.40 °C in the period) percentiles of its observed data values for this group of cells.Grey shaded areas correspond to the 80% HPDI.