Research ArticleECOLOGY

The carbon sink of tropical seasonal forests in southeastern Brazil can be under threat

See allHide authors and affiliations

Science Advances  18 Dec 2020:
Vol. 6, no. 51, eabd4548
DOI: 10.1126/sciadv.abd4548

Abstract

Tropical forests have played an important role as a carbon sink over time. However, the carbon dynamics of Brazilian non-Amazon tropical forests are still not well understood. Here, we used data from 32 tropical seasonal forest sites, monitored from 1987 to 2020 (mean site monitoring length, ~15 years) to investigate their long-term trends in carbon stocks and sinks. Our results highlight a long-term decline in the net carbon sink (0.13 Mg C ha−1 year−1) caused by decreasing carbon gains (2.6% by year) and increasing carbon losses (3.4% by year). The driest and warmest sites are experiencing the most severe carbon sink decline and have already moved from carbon sinks to carbon sources. Because of the importance of the terrestrial carbon sink for the global climate, policies are needed to mitigate the emission of greenhouse gases and to restore and protect tropical seasonal forests.

INTRODUCTION

Tropical forests have a key role in the global carbon dynamics by accounting for one-third of the terrestrial gross primary production and one-half of the terrestrial stored carbon (1, 2). Increasing atmospheric CO2 concentration, rising temperatures, drought events, and deforestation are expected to affect ecosystem functioning through plant physiological responses and forest loss over the coming decades (37). How tropical forests will respond to increasing atmospheric CO2 concentration and climate change are sources of uncertainty in predicting their future carbon stocks and net primary productivity (6, 8). Among tropical forests, the ones under stressful conditions, such as the seasonally dry tropical forests (SDTFs) that endure periodic droughts, can be vulnerable to these environmental changes because they are already at the edge of climate niches to sustain forest formations with high biomass (9, 10). Brazil has been the largest source of carbon emissions from gross deforestation up to 2013: In 2013 alone, 192,000 ha of Caatinga forests (where the largest continuous extent of neotropical SDTFs is located) and 24,000 ha of Atlantic forests were deforested (11). An aggravating factor is that only 6.2% of Brazilian SDTF extent is protected (12). In this context, it is crucial to advance our understanding of the carbon sink of these forests (9, 10).

Over recent decades, the terrestrial carbon sink has been increasing globally (13, 14). This phenomenon is possibly explained by the increases in atmospheric CO2 concentration (CO2 fertilization), which is expected to enhance plant growth (6, 13). Carbon dioxide has a key role in photosynthesis and can potentially increase water use efficiency (WUE) by reducing stomata conductance (15, 16). In this context, the increases in atmospheric CO2 concentration are thought to have enhanced photosynthesis more than rising temperatures have enhanced respiration (17, 18). Hence, drought-related stress on plant growth would be understated by rising atmospheric CO2 concentration. However, the mechanisms involved in the feedbacks among vegetation, atmospheric CO2, and climate are complex. For example, the increase in photosynthesis and WUE led by rising atmospheric CO2 concentration does not necessarily promote stand growth because of the effects of other co-occurring factors (19). The effects of increasing drought, rising temperatures, competition, and physiological acclimation to higher levels of CO2 can constrain tree growth and also lead to tree mortality (3, 5, 20, 21). Therefore, it is difficult to disentangle the effects of climate fluctuations and rising CO2 on carbon dynamics because these factors covary and can interact over time (5).

While recent studies have shown a long-term decline in the Amazon rainforest carbon sink, mostly driven by climate-induced tree mortality (6, 22), others have predicted that tropical rainforest carbon sink may be resilient to climate change in the next decades (18). However, it remains uncertain how Brazilian seasonal forests (which are already exposed to drought), such as deciduous forests (SDTFs) and semideciduous forests, have responded to the increasing levels of atmospheric CO2 and climate fluctuations over time (4, 9, 23). To assess how these forests are behaving over time, we use long-term seasonal forest census data from southeastern Brazil (Figs. 1 and 2) to investigate the long-term trends of carbon stocks, gains, losses, and net carbon sink. We draw our inferences from 95 census intervals nested within 32 sites, ranging between 1987 and 2020 (mean site total monitoring length, ~15 years). The spatial (number of sites and sampled area) and temporal (interval length and total monitoring time) sampling efforts varied among sites and forest types; this variability is depicted in Fig. 2, figs. S1 and S2, and Materials and Methods. The forest sites used here are in advanced successional stages, free from fire, flood, landslides, and human disturbances at least for decades before the first census of each site. The data encompass a wide environmental space and three forest types (deciduous, evergreen, and semideciduous forests) (Fig. 2 and figs. S3 and S4), allowing us to investigate whether forests under different climates have differed in their long-term trends. In addition to the rising atmospheric CO2 concentration [parts per million (ppm)] and CO2 change (ppm year−1) over time, mean annual temperature (MAT) and mean annual precipitation (MAP) have shown an unstable temporal trend in our data (Fig. 3). Therefore, we used time (year) as a proxy of the effects of rising CO2, climate fluctuations, and other unmeasured confounding effects over time. We did this because these factors can be codependent and correlated over time, being difficult to disentangle their individual effects. We fitted statistical models to assess the general long-term trends of carbon dynamics, as well as the long-term trends by forest type, and to test whether climate mediates these long-term trends. More specifically, we tested whether sites under different climate conditions have differed in their long-term trends. In this sense, we expected the long-term trends of sites under drier and warmer conditions to differ from the long-term trends of sites under wetter and colder conditions.

Fig. 1 Spatial location of the sampled sites in South America.

The 32 sampled sites belong to three forest types: deciduous forests (n = 11), evergreen forests (n = 5), and semideciduous forests (n = 16) (Sentinel-2 image).

Fig. 2 Distribution of sites and forest types over the climate space and spatial-temporal sampling efforts.

(A) Distribution of the sites and forest types within the climate space represented by mean annual temperature (MAT) and mean annual precipitation (MAP). The points are census intervals (n = 95), and their sizes are proportional to the site sampled area (mean, 1.05 ha) times interval length (mean, 5 years). (B) Frequency of site sampled areas (n = 32). The red dashed line is the mean of the sampled area among sites (1.05 ha). (C) Frequency of site total monitoring length (n = 32) (year of the last census minus year of first census). The red dashed line is the mean of the total monitoring length among sites (14.7 years).

Fig. 3 Long-term trends of the environmental variables.

(A) MAT, (B) MAP, and (C) CO2 change. The points are census intervals (n = 95). The black dashed curves were fitted with generalized additive models (GAM, including a random effect of site), and the green solid curves were fitted with LMM (including a random effect of site). Note that if the effective degree of freedom (edf) from GAM is equal to 1, then the relationship is linear.

RESULTS

Both MAT and MAP fluctuated over time (Fig. 3, A and B). MAT showed a positive trend (0.04°C by year) (Fig. 3A), while MAP showed a negative trend (−10.2 mm by year) (Fig. 3B). Moreover, CO2 change increased almost linearly with time (Fig. 3C). In general, the carbon stocks increased over time until 2013 (~0.67% by year) and then started to decline (Fig. 4A). During most of the time, the net carbon sink was above zero (positive balance between carbon gains and losses), with a slight negative trend; however, in 2013, the net carbon sink became negative (carbon losses exceeded carbon gains) (Fig. 4B), which explains the carbon stock decline after 2013. In general, the net carbon sink decreased by 0.13 Mg C ha−1 year−1. Carbon gains decreased (~2.6% by year) and carbon losses increased (~3.4% by year) almost linearly over time (Fig. 4, C and D). Predictions for the years 2013 and 2018 together with the observed mean, maximum, and minimum values of each carbon dynamics variable can be found in Table 1.

Fig. 4 Long-term trends of carbon stocks and dynamics.

(A) Carbon stocks, (B) net carbon sink, (C) carbon gains, and (D) carbon losses. In (A), censuses (n = 127) nested within sites (n = 32); in (B) to (D), census intervals (n = 95) nested within sites (n = 32). The black dashed curves were fitted with GAMs (including a random effect of site), and the orange solid curves were fitted with LMMs (including a random effect of site). Note that the slopes of the carbon stock, carbon gain, and carbon loss models were estimated in the logarithmic scale; if the edf (from GAM) is equal to 1, then the relationship is linear.

Table 1 This table shows the observed and predicted (2003 and 2018) means, minimums, and maximums of each carbon dynamics variable.

Predictions were obtained from the GAMs (including a random effect of site) (Fig. 4). Data with census intervals (n = 95) nested within sites (n = 32). Note that the minimum and maximum of the predicted net carbon sink are converged toward the mean because the estimated variance between sites was near zero.

View this table:

The long-term trends by forest type revealed that the carbon stocks of the semideciduous forests increased over time, while the carbon stocks of the deciduous and evergreen forests showed a stable trend, whereby the deciduous forests showed a slight (but nonsignificant) decrease (Fig. 5A). All forest types showed a decline in their net carbon sinks over time, with a stronger decline in deciduous forests (Fig. 5B). Carbon gains decreased and carbon losses increased in all forest types, but these trends were more pronounced in the deciduous forests than in other forest types (Fig. 5, C and D).

Fig. 5 Long-term trends of carbon stocks and dynamics by forest type.

(A) Carbon stocks, (B) net carbon sink, (C) carbon gains, and (D) carbon losses. In (A), censuses (n = 127) nested within sites (n = 32); in (B) to (D), census intervals (n = 95) nested within sites (n = 32). The curves were fitted with LMMs (including a random effect of site). Dashed curves are nonsignificant effects (significance level of 0.05).

In the final model, including climate, soil, and time, the forests under different climate conditions differed in their long-term trends (Figs. 6 to 8). Carbon stocks increased over time, except for the driest and warmest sites, in a way that the positive temporal trend of carbon stocks became weaker with decreasing MAP and increasing MAT (Figs. 6A and 7, A and B). Net carbon sink decreased over time, in a way that its negative temporal trend became weaker as MAP increases and MAT decreases (Figs. 6B and 7, C and D). Carbon gains decreased with time, whereby its negative trend became weaker with decreasing MAT and increasing MAP, and became positive under wet conditions (Figs. 6C and 8, A and B). At the same time, carbon losses increased over the years; the temporal trend of carbon losses became weaker with increasing MAP and decreasing MAT (Figs. 6D and 8, C and D).

Fig. 6 Averaged models (selected models ∆AICc ≤ 4) containing the effects of time, climate, and soil on carbon stocks and dynamics.

(A) Carbon stocks (averaged marginal R2 = 32.5%), (B) net carbon sink (averaged marginal R2 = 29.3%), (C) carbon gains (averaged marginal R2 = 48.6%), and (D) carbon losses (averaged marginal R2 = 28.7%). Dots are the conditional averaged coefficients with their 95% confidence intervals; the coefficients and confidence intervals of variables not included in the averaged model were set to zero. In (A), censuses (n = 127) nested within sites (n = 32); in (B) to (D), census intervals (n = 95) nested within sites (n = 32). The coefficients were estimated using LMMs (including a random effect of site). The coefficients of the carbon stock, carbon gain, and carbon loss models were estimated in the logarithmic scale, and all models were fitted with scaled predictors.

Fig. 7 Interaction effects between time and climate on carbon stocks and net carbon sink.

(A) Carbon stocks, time as predictor and MAP as mediating variable. (B) Carbon stocks, time as predictor and MAT as mediating variable. (C) Net carbon sink, time as predictor and MAP as mediating variable. (D) Net carbon sink, time as predictor and MAT as mediating variable (best model containing MAT − ∆AICc = 4.02). In (A) and (B), censuses (n = 127) nested within sites (n = 32); in (C) and (D), census intervals (n = 95) nested within sites (n = 32). The coefficients were estimated by averaging the LMMs (including a random effect of site) with ∆AICc ≤ 4. Note that the slope of year of net carbon sink differs between (C) and (D) because the effects showed in (D) came from the best model containing MAT. The slopes and interactions of the carbon stock models were estimated in the logarithmic scale, and all models were fitted with scaled predictors.

Fig. 8 Interaction effects between time and climate on carbon gains and carbon losses.

(A) Carbon gains, time as predictor and MAP as mediating variable. (B) Carbon gains, time as predictor and MAT as mediating variable (best model containing MAT − ∆AICc = 12.2). (C) Carbon losses, time as predictor and MAP as mediating variable. (D) Carbon losses, time as predictor and MAT as mediating variable. Data with census intervals (n = 95) nested within sites (n = 32). The coefficients were estimated by averaging the LMMs (including a random effect of site) with ∆AICc ≤ 4. Note that the effect of year on carbon gains differs between (A) and (B) because the effects showed in (B) came from the best model containing MAT. The slopes and interactions were estimated in the logarithmic scale, and all models were fitted with scaled predictors.

The effects of climate have changed over time (Figs. 6 to 9). Carbon stocks increased with MAP and decreased with MAT; these effects became stronger from past to present (Figs. 6A and 9, A and B). The effects of MAP and MAT on the net carbon sink were near zero in the past; however, from past to present, the effect of MAP became positive while the effect of MAT became negative (Figs. 6B and 9, C and D). Meanwhile, the positive effect of MAP and the negative effect of MAT on carbon gains increased over time (Figs. 6C and 9, E and F). The positive effect of MAP and the negative effect of MAT on carbon losses became weaker from past to present (Fig. 6D and 9, G and H). The only significant effect from soil variables was found for soil organic matter (SOM), which had negative effects on carbon gains and carbon losses (Fig. 6, C and D). Site area, which was used as a proxy of edge effects, has not displayed significant effects on the carbon dynamics variables (Fig. 6).

Fig. 9 Interaction effects between time and climate on carbon dynamics.

(A) Carbon stocks, MAP as predictor and time as mediating variable. (B) Carbon stocks, MAT as predictor and time as mediating variable. (C) Net carbon sink, MAP as predictor and time as mediating variable. (D) Net carbon sink, MAT as predictor and time as mediating variable (best model containing MAT − ∆AICc = 4.02). (E) Carbon gains, MAP as predictor and time as mediating variable. (F) Carbon gains, MAT as predictor and time as mediating variable (best model containing MAT − ∆AICc = 12.2). (G) Carbon losses, MAP as predictor and time as mediating variable. (H) Carbon losses, MAT as predictor and time as mediating variable. In (A) and (B), censuses (n = 127) nested within sites (n = 32); in (C) to (H), census intervals (n = 95) nested within sites (n = 32). The coefficients were estimated by averaging the LMMs (including a random effect of site) with ∆AICc ≤ 4. The slopes and interactions of the carbon stock, carbon gain, and carbon loss models were estimated in the logarithmic scale, and all models were fitted with scaled predictors.

DISCUSSION

The net carbon sink in southeastern Brazil’s seasonal forests is declining over time by decreasing carbon gains and increasing carbon losses (Figs. 4 to 8). The carbon sink became a carbon source in 2013 (Fig. 4B), which explains the decline in carbon stocks after this year (Fig. 4A). Among the forests under different climates, the driest and warmest sites are experiencing the most severe decrease in carbon gains, the most severe increase in carbon losses, and, therefore, the most severe decline in the net carbon sink (Figs. 5 to 9). The severe decrease in the carbon sink of the driest and warmest forests suggests that these forests may have reached a climatic stress threshold (9). The carbon stocks of the driest and warmest sites remain stable with a slight (but nonsignificant) decrease (Figs. 5A and 7, A and B); however, if their net carbon sink remains in a negative balance, then their carbon stocks will also decline in the near future, as observed in the general trend after 2013 (Fig. 4A).

Under wet conditions, the carbon gains increased over time (Fig. 8A), consistent with the hypothesized pantropical increase in tree growth caused by CO2 fertilization (6, 18, 24). However, the long-term decrease in carbon gains experienced by the sites under intermediate climate and by the driest and warmest sites (Fig. 8, A and B) is inconsistent with CO2 fertilization effects and with findings in the Amazon forests (6, 18, 24). Recent evidence suggests that atmospheric CO2 increases are not necessarily translated into larger amounts of sequestered carbon by old-growth forest trees (19). A large portion of this CO2 surplus can be emitted back into the atmosphere by processes such as soil respiration (19). Alternative hypotheses might also explain this effect (3, 20, 21). High levels of CO2 can increase tree-to-tree competition by enhancing the growth of some species or individuals, which, at the stand level, can decrease carbon gains by constraining the growth of the suppressed trees (5, 21). In addition, trees growing under increasing CO2 can acclimate to high CO2 availability, by reducing their photosynthetic capacity below the expected for a given CO2 level (20, 25). The potential increases in liana density caused by increasing CO2 and decreasing MAP (Fig. 3, B and C) can also decrease carbon gains by stimulating tree-liana competition for light, water, and nutrients (3, 5, 2628). Meanwhile, the potential direct effects of climate fluctuations, such as decreasing MAP and increasing MAT (Fig. 3, A and B), on carbon gains may have suppressed the effects of CO2 fertilization on photosynthesis and WUE. Because water is an important resource for photosynthesis and high temperatures enhance respiration, drought and rising temperatures can cause physiological stress (e.g., carbon starvation and hydraulic failure) and decrease tree growth (5, 6, 22, 29, 30).

The driest and warmest sites are experiencing the most severe increase in carbon losses over time (Fig. 5D and 8, C and D). This effect is possibly explained by the combined effects of increasing CO2, drought, and higher temperatures (Fig. 3) on tree mortality (5). Drought and high temperatures can directly drive tree mortality through physiological stress, which can be potentialized by the effects of CO2 on tree mortality (5, 31, 32). Increasing CO2 can accelerate the speed at which trees reach large heights, which would increase the rate at which they are exposed to dry upper canopy, lightning, and windthrow and to the physiological aspects associated with larger sizes (5, 3335). In addition, increases in liana density provoked by rising CO2 and drought can also enhance tree mortality by increasing tree-liana competition for light and water (3, 5, 2628). The decrease in carbon gains and the increase in carbon losses of the driest and warmest sites over time, which have a distinct flora, naturally associated with dry and warm conditions, suggest that these forests may have reached a stress threshold due to the effects of increasing drought, temperature, and CO2 (5, 9).

Carbon stocks and carbon gains decreased with MAT and increased with MAP over space (Fig. 9, A and B, and E and F), as expected by theory (6, 29, 36). These effects are evidence of the physiological responses to the harsh conditions imposed by high temperatures and drought, which are related to physiological stress (5, 6, 29). However, the wettest and coldest sites have both greater carbon gains and greater carbon losses (Fig. 9, E to H), consistent with the high-gain high-loss dynamic pattern (6). Therefore, the spatial effects of MAT and MAP on the net carbon sink were near zero most of time and became negative for MAT and positive for MAP from past to present (Fig. 9, C and D). These shifts in the spatial effects of climate occurred because the driest and warmest sites experienced the most severe increase in their carbon losses and the most severe decrease in their carbon gains over time. Thus, the net carbon sink of the driest and warmest sites became smaller than the net carbon sink of the wettest and coldest sites. Although MAT was more important than MAP to differentiate the forest types (fig. S3), MAP was more important than MAT to predict the net carbon sink, carbon gains, and carbon losses (Fig. 6). Soil variables were important to differentiate forest types (figs. S3 and S4); however, only SOM displayed significant effects on carbon dynamics. Forests on soils with lower levels of SOM tend to have higher carbon gains and higher carbon losses (Fig. 6, C and D). SOM is related to soil quality and higher productivity; thus, these effects are not conclusive and should be further investigated. In addition, site area, which was a proxy of edge effects, has not displayed significant effects in the carbon dynamics (Fig. 6). Therefore, the results suggest that climate is the most important predictor of the spatial patterns of carbon dynamics in our study region.

We recognize the limitations of our data, such as space-for-time (sites with only one census interval; table S1) and unbalanced spatial-temporal sampling efforts (Fig. 2 and figs. S1 and S2). However, the negative trend and the negative balance of the carbon sink were a clear pattern in our data and results. In general, these forests are shifting from carbon sinks to carbon sources. Currently, the forests under intermediate climate conditions and the forests under the driest and warmest conditions are already carbon sources, probably because they may have reached a stress threshold. Meanwhile, the carbon sink of the wettest and coldest forests is continually declining. The driest and warmest forests naturally have lower carbon stocks (Fig. 9, A and B), which will decline in the near future if their net carbon sink remains in a negative balance. These long-term trends in carbon dynamics are likely to be influenced by climate fluctuations and rising CO2 (5, 6, 22). However, because these factors are correlated and can interact over time, their mechanistic individual effects and the effects of other unmeasured drivers remain uncertain and should be further investigated. Atmospheric CO2 concentration, temperature, and drought events are expected to continue increasing in upcoming decades, implying that the ecosystem functioning of southeastern Brazil tropical seasonal forests may be under threat (5, 6, 9).

Political actions to mitigate greenhouse gas emissions, together with conservation policies to protect these ecosystems, are needed. We also argue that the driest and warmest sites (deciduous forests, SDTFs) should be further included in conservation policies and that revegetation strategies in agricultural areas can be useful to offset the decline in the carbon sink and stocks. Beyond the political implications, our findings are also useful to improve the predictions of future global carbon sink and to bring knowledge to the carbon cycle of tropical forests.

MATERIALS AND METHODS

Study design

We used tree community data from 32 repeatedly measured permanent forest sites located in southeastern Brazil (Fig. 1). The greatest distances among sites are 900 km (latitude) and 177 km (longitude) (Fig. 1); sites’ altitudes range between 450 and 1475 m above sea level (mean, 898 m); MAP is between 666.4 and 1967.3 mm (mean, 1409.2 mm) (Fig. 2A); MAT is between 16.7° and 25.2°C (mean, 21.1°C) (Fig. 2A); and site area is between 3 and 1200 ha (mean, 78 ha). Our data comprise 95 census intervals (nested within 32 sites), the first census in our dataset was in 1987, and the most recent was in 2020 (table S1). The sites encompass three forest types: deciduous forests (the driest and warmest sites; flat topography; >50% of the individuals lose their leaves in the dry season; 24 census intervals nested within 11 sites; monitored from 2002 to 2020), evergreen forests (the wettest and coldest sites; undulating topography; no deciduousness; 12 census intervals nested within 5 sites; monitored from 1995 to 2020), and semideciduous forests (intermediate climate; undulating topography; 20 to 50% of the individuals lose their leaves in the dry season; 59 census intervals nested within 16 sites; monitored from 1987 to 2020) (table S1). The soils of the deciduous forests have higher levels of phosphorus and cation exchange capacity (CEC) (figs. S3 and S4). All sites are closed canopy forests, in advanced successional stages, free from fire, flood, landslides, and human disturbances at least for the decades preceding the first census of each site. The historical record is based on personal communication with the local people, land owners, or reserve managers. In addition to the historical records, field observations also corroborate the absence of human disturbances, fire, flood, and landslides during the monitoring time. There is no record of wind disturbances, although wind effects cannot be disregarded. Because all sites are in advanced successional stages, there are no differences in the successional stages between forest types (figs. S5 to S22). The diametric structure and diversity distribution over diameter classes are very similar among sites and forest types (figs. S5 to S20) and no significant differences in the community weighted means of wood density (WD) (details on WD below) between forest types were detected (figs. S21 and S22). Sampling (sampled area) was designed to capture intrasite heterogeneity (most subplots are located in the interior of the sites, but edges were also captured by the samples) and spanned from 0.2 to 5 ha (mean of 1.05 ha) (using subplots of 400 m2, 20 m by 20 m, in most cases) in each study site, totaling 33.5 ha (Fig. 2B, fig. S1, and table S1). The sampled plots were measured at least twice in every site (table S1). The number of censuses per site varied between two and eight (mean five censuses per site), and census interval length varied among sites and within the same site over time, ranging from 1 to 9 years (mean, 5 years) (table S1). The total monitoring length per site varied from 5 to 31 years (mean, 15 years) (Fig. 2C and fig. S2). Forest inventory data are hosted in the ForestPlots.net system (www.forestplots.net).

All trees with quadratic mean diameter at reference height (1.3 m) [diameter at breast height (DBH)] ≥5 cm were measured in each census (37). This inclusion criterion allows to include tree individuals with multiple stems (stems with DBH ≤ 5 cm) if the quadratic mean diameter at reference height of the tree individual meets the criterion. The height of the point of measurement (POM) was marked in the stems and used as a reference for subsequent measures. When changing the POM was necessary for a given stem, its diameter growth was estimated using a ratio between the current and previous POMs (38). Plant identification followed Angiosperm Phylogeny Group IV (39) and was performed by experts in the field or by consulting herbaria, and species names were standardized using Flora do Brasil 2020 (http://reflora.jbrj.gov.br/reflora/floradobrasil). To obtain WD values, the species were matched to the global WD database (40). When WD was not available at the species level, the average WD of other species from the same genus or family was used.

The aboveground woody biomass (AGWB) of each tree was obtained with the modified pantropical allometric equation of Chave et al. (41) (Eq. 1) through the BIOMASS R package (42). We used the equation with the parameters for tree diameter, WD, and environmental stress (E) because tree height data are absent in our censuses. Equation 1 is broadly applied to tropical forests worldwide for deriving from a large dataset encompassing the pantropical region of the globe. But another equation (Eq. 2), proposed by Sampaio and Silva (43), was developed specifically to estimate tree aboveground biomass of deciduous tropical forests in the Brazilian Caatinga domain, which is the case for the deciduous forests of our dataset. Therefore, to reduce the uncertainty in the biomass estimates of the deciduous forests, we combined (averaged) the estimated values of Eqs. 1 and 2 for deciduous forest trees. The correlation between the estimated values of the two equations is 0.81 (Pearson correlation), but the mean of the estimated values from Eq. 1 was 2.2% greater than the estimated values from Eq. 2. The AGWB was converted to aboveground carbon (AGC) assuming that AGC is 45.6% of the AGWB (44). We estimated the carbon stock (Mg C ha−1) of each site in each census as the sum of each living tree’s AGC, scaled to hectares. Carbon gains (Mg C ha−1 year−1) were estimated as the sum of the growth of the surviving trees (tree AGC at the end of the interval minus the tree AGC at the start of the interval) plus the sum of the AGC of the recruited trees (tree individuals or stems that reached DBH ≥ 5 cm at the end of the interval, assuming DBH = 0 at the start of the interval). The sum of the growth of the survivors with the AGC of the recruited trees was divided by the interval length (in years) and then scaled to hectares. Carbon losses (Mg C ha−1 year−1) were estimated as the sum of the AGC of the individuals and stems that died during a census interval, divided by the census interval length (in years), and then scaled to hectares. As carbon gains and losses are biased by census interval length, we corrected their values using the equation CIC1 proposed by Talbot et al. (38). We calculated net carbon sink (Mg C ha−1 year−1) as carbon gains minus carbon losses. Note that the length of each census interval was calculated using the rounded years (i.e., year of the end of the interval minus the year of the start of the interval).AGB (Mg)=exp(2.0240.896×E+0.920×ln(WD)+2.795×ln(DBH)0.0461×(ln(DBH)2))(1)andAGB (kg)=0.1730×DBH2.295(2)where DBH is the diameter at breast height (cm), WD is the wood density (g cm−3), and E is the environmental stress.

Climate data

We obtained monthly mean temperature (°C), monthly average daily maximum temperature (°C), monthly precipitation (mm), and daily potential evapotranspiration (PET; mm) from the Climatic Research Unit (CRU TS version 4.04, released 24 April 2020; https://crudata.uea.ac.uk/cru/data/hrg) (45). We used data from WorldClim 2.1 (46) to downscale the CRU temperature and precipitation data to 1 km2 and then applied the monthly correction for all months in each census interval. Delta spatial downscaling method was used for downscaling; see Peng et al. (47) for details. We calculated annual maximum climatological water deficit (CWD) by summing the differences between monthly downscaled precipitation (mP) and monthly PET (PET minus mP) only when this difference was positive (i.e., evapotranspiration exceeded precipitation: water deficit). We used the average values of MAT, maximum temperature (max temp), CWD, and MAP of the years within each census interval. For the sites measured in 2020, for which climate data were not available in CRU TS version 4.04, we used climate data ranging from the start of the interval to the end of 2019. The climate space and the temporal trends can be found in Figs. 2A and 3 (A and B). Note that for the first census of the census level data (carbon stocks), we used the averaged climate values from 5 years before the first census of each site. Max temp is strongly correlated with MAT (r = 0.99) and MAP is strongly correlated with CWD (r = −0.96); thus, we opted to use MAT and MAP because these variables are easier to obtain and are usually more accessible than max temp and CWD.

Atmospheric CO2 concentration

We used annual mean values of atmospheric CO2 (ppm) concentration from the Mauna Loa record (48). We also calculated the rate of CO2 change (ppm year−1) by subtracting the CO2 concentration at the end of the interval by the CO2 concentration at the start of the interval and then dividing this difference by the census interval length. Both CO2 concentration and CO2 change increased over time (Fig. 3C and figs. S23 and S24). The Pearson correlation between CO2 concentration and year is 0.99 while between CO2 change and year the correlation is 0.85. Therefore, we opted to use time (year) in the models instead of CO2 variables. We did this because the potential effects of CO2 are not distinguishable from the effects of other environmental variables over time, such as climate fluctuations and other unmeasured potential predictors.

Soil data

To control the potential effects of soil, we collected soil surface samples (20 cm of depth) in the first census of each site. Soil surface samples were collected from three points in each subplot and later combined into one composite sample. Following the protocol by the Empresa Brasileira de Pesquisa Agropecuária (49), the following soil attributes were obtained: pH in water (pH), phosphorus (mg cm−3) (P), potassium (mg × cm−3) (K), calcium (cmolc × dm−3) (Ca), magnesium (cmolc × dm−3) (Mg), aluminum (cmolc × dm−3) (Al), sum of bases (SB) (cmolc × dm−3), cation exchange capacity (cmolc × dm−3) (CEC), soil organic matter (dag × kg−1) (SOM), and clay percentage (%) (clay). We calculated the mean of each soil variable for each site.

Statistical analysis

All analyses were carried out in the R environment version 3.6.1 (50). Because of the nonindependence between observations within the same site (temporal autocorrelation; repeated measurements in the same sites over time), our modeling and statistical framework were based on tools that control for the nonindependence between observations. We used generalized additive models (GAMs) with a random effect of site and linear mixed-effects models (LMMs), which allow the intercepts to vary randomly among sites. We used the mgcv package (51) to fit the GAMs and the lme4 package (52) to fit the LMMs. We did not include random slopes because some sites have few census intervals, preventing the models to achieve convergence. All models were fitted under normality assumptions (Gaussian family). We assessed residual normality and homoscedasticity through residual inspection. We applied logarithmic transformation on carbon stocks, carbon gains, and carbon losses to meet normality assumptions. As plots vary in spatial and temporal sampling efforts (Fig. 2), we weighted the observations (prior weights) by sampled area and census interval length. For carbon stock models, we weighted by the cubic root of the sampled area (36), and for the other response variables, we weighted by the cubic root of census interval length plus the fourth root of sampled area minus one (6).

The modeling framework can be divided into three parts: (i) the general long-term trends of carbon dynamics (carbon stocks, net carbon sink, carbon gains, and carbon losses), (ii) the long-term trends of carbon dynamics by forest type, and (iii) a final model including time (as a proxy of the effects of rising CO2, climate fluctuations, and other unmeasured confounding effects over time), climate and soil, and the interactions between time and climate. The interactions between climate and time allowed to evaluate whether sites under different climate conditions differ in their long-term trends and whether the effects of climate have changed over time. Note that we have not included climate and forest type in the same model because forest type is correlated with climate and soil (Fig. 2 and figs. S3 and S4), leading to high variance inflation factor (VIF) values when they are in the same model. We used GAMs only in part 1 because we do not have sufficient data to fit the models of parts 2 and 3 with smoothing splines (overfitting). In all models, we used the year corresponding to the midpoint of each census interval (except for carbon stock model, to which we used census year).

To estimate the long-term trends of carbon dynamics variables (bivariate relationships) (part 1), we used GAMs and LMMs. More specifically, we regressed the carbon dynamics variables as function of time (year). We estimated the long-term trends of each forest type using LMM (part 2). Because of the strong multicollinearity between the environmental variables, before building the global model of part 3, we removed collinear variables to avoid redundant and circular interpretations. In general, we removed variables with Pearson correlation (r ≥ |0.75|), with the following exceptions: Although strongly correlated (r ≥ |0.75|), we maintained MAP and MAT (r = −0.89), and CEC as a proxy of soil fertility (which is strongly correlated with climate) (see correlation matrix at fig. S24 for details). We also controlled for the potential effects of site area (proxy of edge effects) by including the area of the site (site area) in the global model. Thus, the final set of variables included in the global model was year, MAP, MAT, P, CEC, SOM, clay, and site area. The global model was built in a way that the independent effect of each explanatory variable was accounted, as well as the interactions between year with MAP and MAT (Eq. 3). All predictors were scaled and centered to zero mean and unit variance. In addition, we checked for spatial autocorrelation in global models’ residuals using Moran’s I test, implemented in package ncf (53).Carbon dynamicsyear×(MAT+MAP)+P+CEC+SOM+clay+site area(3)where carbon dynamics is carbon stocks, net carbon sink, carbon gains, and carbon losses. Additive and interaction effects are represented, respectively, by “+” and “×.” Note that the independent effects of year, MAT, and MAP are included in the interaction term.

An information theoretical approach based on the Akaike Information Criterion of second order (AICc) was used for model selection (54). From the global model (Eq. 1) of each response variable, we obtained the set of best models (those with ΔAICc ≤ 4) (54). To avoid multicollinearity issues, the selected models were constrained to have explanatory variables with r ≥ |0.6| (55), ensuring low VIF values (VIF ≤ 4). To avoid overfitting, we also limited the best models’ degrees of freedom, ensuring at least 15 observations per term. Using multimodel inference, we averaged the coefficients of the selected models and used the conditional averaged coefficients as a final result (54). The relative importance of the predictors was not considered, given that some variables were not contained in the same number of models because of collinearity, which could bias the sum of Akaike weights (54). Our conclusions relied on the significant conditional averaged coefficients (significance level of 0.05). The MuMIn package (56) was used for model selection, model averaging, and to obtain the marginal R2 (variance explained by the fixed effects, we used the average between the marginal R2 of the best model and the marginal R2 of the global model) (57). Graphics were obtained through the packages ggplot2, viridis, and corrplot (5860).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/51/eabd4548/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank the Federal University of Lavras (UFLA) and all students who collected the data used in this study over time. Funding: This study was funded by the State of Minas Gerais Research Foundation (FAPEMIG), the National Council for Scientific and Technological Development (CNPq), and the Coordination for the Improvement of Higher Education Personnel (CAPES). Author contributions: V.A.M. and R.M.d.S. conceived the ideas and designed methodology. A.B.M.S., C.R.d.S., T.A.V., L.C.A.d.S., R.T.P., N.C.A.F., G.G.P.d.P., P.A.C., P.O.G., M.A.L.F., F.M.G., F.d.O., D.T.G., J.D.M., P.F.S., F.d.C.A., N.d.A.-C., C.L.F., G.C.d.O.M., M.A.L.F., W.B.d.S., R.M.d.S., and V.A.M. collected the data. V.A.M. analyzed the data. V.A.M., A.B.M.S., and N.d.A.-C. led the writing of the manuscript. C.L.F. and C.R.d.S. curated the data. P.A.C., N.d.A.-C., M.C.F.d.O., L.C.A.d.S., and J.D.M. edited the manuscript. All authors contributed critically to the drafts and gave final approval for publication. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Forest inventory raw data are hosted in the ForestPlots.net system and bound to data-use restrictions defined on Forestplots.net (www.forestplots.net) (see area codes in table S1). All data needed to evaluate the conclusions in the paper are present in the Supplementary Materials (data files S1 and S2). The R codes are available upon request from the corresponding author.
View Abstract

Stay Connected to Science Advances

Navigate This Article