We used data of forest dynamics from 95 plots encompassing 432 censuses from a published dataset plus data obtained via ForestPlots.net (Hubau et al., 2020; ForestPlots.net et al., 2021; Figure 2). Censuses took place between 1981 and 2017, including drought and non-drought periods. Plot size ranged from 0.25 to 9 ha, with a mean size of 1.14 ha and total plot area of 108 ha. Census intervals ranged from 4 months to 5.6 years, and plots have been monitored for an average of 15 years. For each census interval, we calculated Aboveground Biomass (AGB) gains caused by growth and recruitment, AGB losses from mortality and AGB net change as the difference between AGB gains and AGB losses. To calculate AGB, we used a community weighted mean wood density dataset available on ForestPlots.net that is derived from a global wood density database (Lopez-Gonzalez et al., 2011). This same database was used in the analyses which include community weighted mean wood density. AGB rates were transformed into Aboveground Carbon (AGC) stocks using the ratio 2:1 (IPCC, 2006). We excluded censuses with an interval length longer than twice the mean value of all censuses (>5.77 years) from our analysis. We exclude these censuses because long intervals are less likely to detect any effects of drought. All plots are in mature lowland forest (< 1,500 m above sea level), in non-ADE sites, without known recent fire disturbances or selective logging (i.e., from 1981) and represent a climatic and edaphic variation gradient across the Amazon Basin.

Soil samples were collected from the same plots as the forest inventory data (Figure 2) at depths of 0–5, 5–10, 10–20, and 20–30 cm as part of past studies (Koele et al., 2017; Quesada et al., 2020) and this study with soil samples taken between 2017 and 2019. In plots ≥1 ha we sampled a minimum of five different locations across the plot, and in plots <1 ha we sampled three locations. Then, we combined the samples in each plot by depth to result in one sample per plot per depth. For each sample, we measured the total organic carbon (TOC), pyrogenic carbon (PyC), and non-PyC fraction of the organic carbon (OC). The hydrogen pyrolysis technique (HyPy) (Ascough et al., 2009) was used to quantify the PyC fraction of TOC, which represents stable polycyclic aromatic carbon with a ring number of >7 (Meredith et al., 2013). However, there is potentially more PyC that is not quantified by this technique since PyC produced by forests fires includes a variable fraction of PyC with an aromatic ring size of <7 that is not quantified by this technique. The non-PyC fraction of organic carbon (OC) is the difference between TOC and PyC and will include low temperature PyC not quantified by HyPy (if present).

We calculated the weighted mean for all depths to estimate the PyC and non-PyC OC concentrations in the soil surface (0–30 cm). The statistical analyses use continuous values; for Figure 5 we assigned each plot to one of two PyC classes based on the PyC median distribution: “Low”– values smaller than 0.048% and “High”– values equal or higher than 0.048%.

The Maximum Cumulative Water Deficit (MCWD) is a well-established metric to indicate water stress in forests (Aragão et al., 2007). MCWD is calculated from the cumulative difference in precipitation and evapotranspiration using an average evapotranspiration (E) for lowland moist tropical forests of approximately 100 mm per month (da Rocha et al., 2004; von Randow et al., 2004). Consequently, if the monthly precipitation (P) is lower than 100 mm, the forest enters a state of water deficit (WDef), which can accumulate over time (CWD). To calculate CWD, we used monthly precipitation data from the Climate Hazards Group InfraRed Precipitation with Stations (CHIRPS) from January 1981 until December 2017 (Funk et al., 2015). This dataset integrates satellite data with in situ rainfall gauge station measurements to produce gridded rainfall time series, and performs well across several regions of the world (Duan et al., 2016; Katsanos et al., 2016; Verdin et al., 2016; Perdigón-Morales et al., 2018). In the Brazilian Amazon, the CHIRPS dataset explained 73% of the observed rainfall in field gauge stations with a root mean square error below 15 mm per month (Anderson et al., 2018). The following rule was applied to calculate the CWD per plot (p) for each month (n), using the evapotranspiration fixed at 100 mm month^–1:

The most negative CWD value among all months in each calendar year was taken as the MCWD. Subsequently, we calculated MCWD for each census interval, which refers to the most negative CWD during the census interval. The mean annual MCWD from 1981 to 2017 for each plot was used as the baseline to compute MCWD anomalies (Z-score) for each census interval, normalized by the standard deviation (σ). Censuses with anomalies ≤−1.65 were classified as severe drought events, based on a confidence level of 90% (Aragão et al., 2018).

Soil physicochemical data derived from plot samples were provided by ForestPlots.net (ForestPlots.net et al., 2021). We used data for surface (0–30 cm) soil including soil texture (silt, clay, sand fractions), total P, pH and exchangeable cations (K, Mg, Ca). The soil samples were all analyzed following the protocol described in Quesada et al. (2010).

All soil fertility variables were centered, scaled and condensed into one axis using a principal component analysis (PCA) using the R package FactoMineR (Husson et al., 2016). The first eigenvalue explained 73.3% of the variance and was used as a single soil fertility variable in the analysis, visualized using the R package factoextra (Kassambara and Mundt, 2017) (Supplementary Figure 1 and Supplementary Table 1).

To evaluate the relationship between soil PyC and soil fertility, OC, wood density and soil texture we used Pearson’s correlation tests. We applied a natural log transformation to soil PyC, soil fertility and OC to meet assumptions for data normality.

We applied a generalized linear mixed effects model (GLMM) to evaluate the relationship between soil PyC and AGC gains, losses and net change. To account for census replication per plot and for regional differences across the Amazon Basin, we used the plot code nested in plot clusters as a random intercept effect in the models, with a fixed slope. The plot clusters refer to plot codes that are located within the same local group and are represented by the same three first letters of the plot codes (Lopez-Gonzalez et al., 2011).

To evaluate the relationship of soil PyC on the change in AGC dynamics with drought severity, we filtered our database to only include censuses classified as severe drought (MCWD anomalies ≤−1.65σ) (Feldpausch et al., 2016). In this analysis, we also applied a GLMM using the same random effect cited above. Soil PyC, severe MCWD anomalies, OC, soil fertility, wood density and interaction between variables were used as predictors to evaluate differences in AGC dynamics according to drought severity. We applied a square-root transformation of AGC loss, and a log transformation of soil PyC and MCWD anomalies to ensure a normal distribution of fitted residuals. All variables were standardized to allow comparisons between model outputs. We applied a backwards stepwise selection procedure and selected the best models according to the lowest Akaike Information Criterion corrected for sample size (AICc) values using the MuMIn package (Barton, 2020). To evaluate spatial autocorrelation that could bias the models’ coefficients, we used the models’ residual to plot spline correlograms with 95% confidence interval based on 1,000 bootstrap resamples using the spline.correlog function in the package ncf (Bjørnstad and Falck, 2001). In order to understand whether the effects of soil PyC were driven by non-PyC organic carbon soil amelioration, we re-ran the models with OC instead of soil PyC. All GLMM models were performed using the lme4 package (Bates et al., 2018). All statistical analyses were conducted using R statistical software version 4.0.3 (R Core Team, 2020).

Soil PyC was significantly correlated with all measured variables. We found highly significant positive correlations between soil PyC with soil fertility (PCA axis-1) (r = 0.49, p < 0.001), OC (r = 0.30, p = 0.003), clay (r = 0.52, p < 0.001) and silt (r = 0.34, p = 0.002), and a negative correlation between soil PyC with wood density (r = −0.33, p < 0.001) and sand (r = −0.6, p < 0.001) (Figure 3).

Our results showed no significant relationship between soil PyC and AGC gains, AGC losses or AGC net change (Figure 4 and Supplementary Table 2). There was no significant autocorrelation in the model residuals and thus no spatial bias in the model results (Supplementary Figure 2).

To investigate the effect of PyC on AGC dynamics during drought events, we selected only censuses that had experienced a severe drought (i.e., MCWD anomalies ≥−1.65σ). Soil PyC, MCWD anomalies, OC, soil fertility, wood density, and interactions between variables were used as predictors to evaluate AGC dynamics with drought severity (Supplementary Table 3). The best models selected to explain AGC dynamics based on AICc values are reported in Table 1. Spline correlograms of the best models showed no significant spatial autocorrelation in the model residuals (Supplementary Figure 3). Wood density and the interaction of soil PyC with MCWD anomalies were significant predictors of AGC gains (i.e., woody productivity) during severe droughts. Our results showed lower rates of AGC gains in forests with greater wood density. We also found a significant interaction between soil PyC and MCWD anomalies on AGC gains (p = 0.041), with the slope between AGC gains and MCWD anomalies becoming less negative with increasing soil PyC (Table 1).

We found the difference in AGC gains between forests with PyC concentrations lower and those with PyC concentrations greater than the median (0.048%) to increase with drought severity, as AGC gains significantly declined at a greater rate in low soil PyC forests compared with high soil PyC forests (Figure 5). We compared events around the mean MCWD anomaly (∼−2σ, p < 0.046) and the most severe MCWD anomaly (∼−4.3σ, p < 0.001) for forests with the first quartile of observed soil PyC concentration (0.038%) to forests with the third quartile of observed soil PyC concentration (0.056%), which was ∼1.5 times greater. Our models suggest that forests with the third quartile of observed soil PyC concentrations have 3.8% greater AGC gains under mean drought conditions, but 33.7% greater AGC gains under the most extreme drought conditions compared with forests with the first quartile of observed soil PyC concentrations. These results may translate to a difference in the gain of up to 0.68 Mg C ha^–1yr^–1 during a most extreme drought.

The best model to explain AGC losses also included WD and the interaction between soil PyC and MCWD anomalies. However, the effect size did not significantly differ from zero. We could not explain any of the variance in AGC net change during severe droughts with the tested variables (Table 1). We also found the non-PyC OC had no significant effect on AGC gain and loss with drought severity when tested instead of soil PyC (Supplementary Table 4).

Our findings indicate that soil PyC is associated with increased soil fertility (Figure 3A). Studies of ADE and biochar (PyC) addition indicate that PyC can improve soil fertility via several mechanisms (Glaser et al., 2001; Czimczik and Masiello, 2007). PyC has a polycyclic aromatic structure and can persist in the environment for centuries to millennia (Bird et al., 2015). Therefore, PyC can reduce aluminum toxicity by increasing pH, increase cation exchange capacity, and improve water holding capacity (Glaser et al., 2001; Czimczik and Masiello, 2007). The ADE are historical indigenous lands where fire was frequently used as a land-management tool. Covering only 3% of Amazonia (McMichael et al., 2014), these small-scale areas (e.g., 0.5–300 ha in size) (Paz-Rivera and Putz, 2009) are more fertile and richer in soil PyC than adjacent areas (Glaser et al., 2002; Liang et al., 2006; Glaser, 2007). Our study highlights the impact of soil PyC-associated changes in soil fertility for non-ADE soils, representing the remaining 97% of the Amazon Basin, where we found a positive correlation between PyC and soil fertility.

We quantified total resistant PyC, representing PyC produced in situ by local fires and PyC that may have been deposited by aeolian transport from remote fires. The aerosol-derived PyC is derived semi-continuously from ancient to modern fires and fluxes have recently been estimated (fossil fuel + biomass burning) at approximately 6 kg km^–2 yr^–1 in Amazonia (Coppola et al., 2019). Given this rate, it would still take approx. 27,000 years to accumulate the estimated store of 1.1 Pg (0–30 cm) of PyC in the Amazon Basin from aerosol deposition alone (Koele et al., 2017). This suggests that the majority of the PyC analyzed derives from ancient local fires. Atmospheric transport, over centennial-scales, may also have been important, and these processes combined may have improved soil fertility over time.

Soil texture is a key determinant of soil fertility. There is a large range of physical and chemical soil properties across the Amazon Basin that vary according to gradients of pedogenic development, where nutrient pools are lowest in the most weathered soils (Quesada et al., 2010). We found the same pattern for soil PyC, with soil PyC concentration negatively correlated with soil sand content (Figure 3D), and occurring in highly weathered soils. Because of macropore predominance, sandy soils are more aerated and, consequently, destabilize PyC through high rates of oxidation and potentially greater rates of vertical translocation (Major et al., 2010). On the other hand, soils that experienced less weathering and, consequently, have more clay and silt, showed a positive association with PyC (Figures 3E, F) because of lower PyC oxidation rates. Our results corroborate other studies showing that soils with >50% clay content have significantly more PyC (Reisser et al., 2016). Moreover, clay-rich soils present more opportunities for organo-mineral interactions, helping to stabilize PyC (Sørensen, 1972; Six et al., 2002; Reisser et al., 2016).

We also found a negative correlation between wood density and soil PyC (Figure 3C). A previous large-scale fire experiment has already shown that tropical soft-wood trees are at greater risk of death from fire, in part because tree species with low wood density are less likely to close wounds post-fire (Balch et al., 2015). Moreover, wood density had a significant positive relationship with tree survival rates in areas burned repeatedly, where increasing wood density by 0.8 g cm^–3 can enhance tree-level survival probability by 15% (Brando et al., 2012). Note that only one previous work has investigated the relationship between soil PyC and wood density and did not find a significant tendency (Massi et al., 2017). However, it included fewer than half the number of plots used in our study.

It is well established that wood density is negatively correlated with soil fertility (Baker et al., 2004; Quesada et al., 2012). We found a strong correlation between soil PyC and soil fertility (Figure 3A), and other previous studies have shown soil PyC can increase soil fertility as discussed previously in this section. Despite not being able to separate the origin of soil PyC (aerosol deposition or locally produced) in our analyses, we hypothesize that forests with high soil PyC concentrations added by aerosol deposition are more likely to have low wood density. This may be because soil PyC can increase soil fertility, leading to forests with fast stem turnover (Quesada et al., 2012). Wood density is also known to be related to successional stage and thus disturbance history. Meanwhile, forests with high concentrations of locally produced soil PyC may have experienced ancient fire disturbances in situ. Low wood density in these forests with high concentrations of soil PyC may be an indicator of disturbance and an earlier successional state, rather than merely because of changes in soil fertility. We attribute these differences to ancient fires because there was no evidence of recent disturbance events, although it should be noted that modern records of disturbance only go back <100 years.

Whilst we found correlations between soil PyC and soil fertility, soil texture and tree wood density, we were unable to directly test causal relationships in this study and these results should therefore be interpreted with caution. It should be noted that the highest concentrations of soil PyC were found in the West Amazon (Figure 2) that also has the most fertile soils, the lowest sand content and the lowest stand-level wood density (Quesada et al., 2012).

Our results do not show a relationship between soil PyC and AGC dynamics when analyzed across all census intervals. This result may be related to low concentrations of soil PyC found in the majority of the plots analyzed in this study. A study from eastern and southern Amazonia showed that biomass and forest composition had legacy effects in forests on or near ADE soils with high soil PyC and a history of ancient fires, which in general presented higher aboveground biomass (de Oliveira et al., 2020). It is likely that soil PyC has stronger relationships at high concentrations since the effect of PyC on soil fertility and water holding capacity is likely to be small at the concentrations found in this study (Glaser et al., 2002; Glaser, 2007). Moreover, disturbances in ADE are likely to be much greater (de Oliveira et al., 2020) and therefore they may still have differences in species composition and traits that would not have persisted in forests that were less disturbed.

Drought events can occur because of natural processes that are related to changes in sea surface temperature, e.g., Atlantic Multidecadal Oscillations, El Niño Southern Oscillations and Pacific Decadal Oscillations (Marengo and Espinoza, 2016; Aragão et al., 2018). However, anthropogenic actions such as land-use and land-cover changes and greenhouse gas emissions may exacerbate the intensity of droughts by changing patterns of large-scale atmospheric circulation (Spracklen and Garcia-Carreras, 2015; Llopart et al., 2018). Drought events in the Amazon are becoming more frequent, prolonged and intense as a result of both globally and locally driven climate change (Malhi et al., 2008; Dubreuil et al., 2012; Bonini et al., 2014). Amazonian forests are vulnerable to drought, which causes biomass losses and reduced forest productivity, enough to temporarily reverse a large multi-decadal aboveground carbon sink into mature forest biomass (Phillips et al., 2009; Feldpausch et al., 2016; Anderson et al., 2018; Gatti et al., 2021).

Our results show that during extreme drought events, forests with greater soil PyC have significantly higher rates of AGC gain compared to forests with lower soil PyC, with this difference increasing with drought severity (33.7% greater gains under the most extreme drought conditions) (Figure 5). These results may be driven by the capacity of soil PyC to hold more water in the soil (Glaser et al., 2002; de Melo Carvalho et al., 2014), alleviating the effects of severe droughts, and allowing trees to continue to grow under drought conditions. This significant effect on AGC gains during drought events in forests with higher soil PyC compared to those with lower soil PyC, may also be related to soil fertility, since forests with high concentrations of soil PyC are also more fertile soils. The drought immunization response is associated with the PyC fraction of TOC and not the non-PyC fraction of the organic carbon, since we found no significant relationship between OC and changes in AGC gains and losses with drought severity. Our models also suggest that PyC rather than soil fertility drives the patterns we observed. This indicates that the responses we found are potentially related to the effects of ancient fire on soil and vegetation (Quesada et al., 2012).

Some studies have shown that during extreme droughts, old-growth forests can lose more AGC than forests disturbed by past fires (Brando et al., 2014; Berenguer et al., 2021). Since old-growth forests have trees with higher wood density than those in disturbed forests, the death of few large trees in an old-growth forest can cause a greater amount of carbon loss than in disturbed forests (Brando et al., 2014; Berenguer et al., 2021). Therefore, our results may demonstrate that forests with a stronger history of fires, here identified by high concentrations of soil PyC, can be more resistant to current drought events. These forests may have lower overall biomass resulting from an establishment of long-term successional species, and consequently do not show a significant reduction in biomass growth during drought events.

In contrast to finding a significant relationship between PyC and AGC gains, we were unable to explain any of the variance in AGC losses (Table 1). Berenguer et al. (2021) found that forests increase the carbon loss rates for up to 3 years after a drought, and consequently, forest inventories taken in the first year after a drought may not detect the full impact of drought on tropical vegetation. Since many of our census intervals did not capture 3 years of post-drought dynamics, we may be failing to capture any effect of PyC on long-term carbon losses. Moreover, high spatial and temporal heterogeneity in tree mortality may prevent a signal from being detected in this study (Johnson et al., 2016; Pugh et al., 2020). Many of our plots were small (≤1 ha) or had short census intervals (<2 years), making it difficult to estimate AGC losses precisely. Much larger datasets may be needed to separate a trend from the large natural variance in tree mortality. Despite an increase in AGC gains, we are not able to detect an increase in AGC net change. Since spatial variation in carbon stocks depends more on losses than gains (Johnson et al., 2016; Hubau et al., 2020; Pugh et al., 2020), the caveats outlined above are also likely to hold true for the AGC net change data.

It is critical to understand the impacts of ancient fire on contemporary forest dynamics. Our results highlight the importance of understanding ancient fire regimes when predicting how species composition and carbon storage will change as a result of drought. Our results suggest ancient fires may influence current forest dynamics by altering soil fertility and favoring tree species capable of continued growth and recruitment during droughts. A major shift in the frequency and intensity of fires and droughts has occurred this century, bringing large uncertainties for future predictions of the carbon cycle (Aragão et al., 2007; Aragão et al., 2018; Silva Junior et al., 2019). This increase in fire represents a new fire regime for the Amazon (Aragão et al., 2018; Silveira et al., 2020). These fires bring several impacts to tropical rainforests such as changes in forest structure, species composition and carbon dynamics (Sato et al., 2016; Prestes et al., 2020; Silva et al., 2020; Pontes-Lopes et al., 2021). Forests that burned 30 years ago still have ∼25% less aboveground biomass than unburned forests (Silva et al., 2018), showing these impacts can persist over decades. From our analysis using soil PyC as a proxy and legacy of ancient fires, it is possible that forests which burned centuries ago may still not have recovered to their pristine state but instead continue to exhibit some attributes associated with recently disturbed or secondary forests (Berenguer et al., 2018; de Oliveira et al., 2020; Heinrich et al., 2021). Our study, however, suggests that mature forests that experienced fires centuries or millennia ago have greater resistance to short-term droughts, as a consequence of ancient fires. Further studies and experiments are needed to identify whether the effects of PyC on fertility and/or water holding capacity drive the patterns we observe, or whether PyC simply acts as a proxy of ancient fires.