Soil–Atmosphere GHG Fluxes in Cacao Agroecosystems on São Tomé Island, Central Africa: Toward Climate-Smart Practices

1. Introduction

Cocoa (Theobroma cacao L.) is one of the most important tropical crops, predominantly cultivated in West and Central African countries, which together accounted for approximately 71.3% of global production—around 3.12 million tons of cocoa—during the 2023–2024 period [1]. However, the carbon footprint of cocoa production remains considerable, with average emissions estimated at 1.47 kg CO₂e per kilogram of cocoa produced. Deforestation is identified as the main driver of these emissions, while tree biomass and the adoption of good management practices can improve the carbon balance, aligning profitability with low climate impact [2,3]. Although cocoa is typically grown in agroforestry systems known for their agro-environmental benefits, drivers such as land-use change, unsustainable conventional practices, and climate change can adversely impact the net balance of greenhouse gas (GHG) emissions, including carbon dioxide (CO₂), methane (CH₄), and nitrous oxide (N₂O) [2,4,5,6].

According to the 2019 annual report on GHG emissions monitoring and verification [7], São Tomé's estimated emissions of the main GHGs in 2018 reached 204.98 Gg CO₂-eq (excluding FOLU), representing a 31% increase compared to the last inventory conducted in 2012, while carbon sequestration increased by 23%. Emissions from the agriculture/livestock/forestry and other land use (AFOLU) sector were estimated at 24.41 Gg CO₂-eq. Among its GHG reduction targets, São Tomé projected a goal of reducing total GHG emissions by 109 KtCO₂-eq, nearly doubling the initial target submitted in 2015. This target corresponds to an approximate 27% cut in national GHG emissions compared to the business-as-usual scenario by 2030, as outlined in its Nationally Determined Contributions (NDCs) [8].

The accumulation of GHGs in the atmosphere is driving a sustained increase in global temperatures [9,10]. This phenomenon is widely recognized as a major driver of climate change, leading to alterations in precipitation patterns, sea-level rise, and other critical ecological processes [11,12,13]. Agriculture is a major contributor to global GHG emissions, mainly due to land-use change, deforestation, fertilizer application, enteric fermentation, and manure management [14,15].

CO₂ emissions from agriculture are primarily associated with fossil fuel combustion, soil degradation, and forest conversion for cropland expansion [11,16]. In addition, soils naturally release CO₂ through microbial and root respiration. These processes are highly sensitive to soil temperature, moisture, nutrient availability, and organic matter quality [17,18,19]. These emissions can also be modulated by soil texture, vegetation cover, and topographic position [20,21,22,23].

Methane (CH₄), the second most important biogenic GHG, is mainly produced under anaerobic conditions through enteric fermentation and organic matter decomposition in poorly aerated soils [11]. However, upland soils can act as CH₄ sinks through methanotrophic activity, especially in well-aerated environments such as tropical uplands [24,25]. In contrast, periodically saturated or waterlogged soils, such as those in lowland agroforestry systems, can become net sources of CH₄ [26]. Recent studies have highlighted how CH₄ fluxes vary according to land cover, soil physical conditions, and vegetation structure [12,27,28].

Although present in lower concentrations, nitrous oxide (N₂O) has a global warming potential 298 times greater than CO₂ [9] and is strongly linked to microbial nitrification and denitrification processes [29]. In tropical agricultural systems, N₂O emissions are often exacerbated by excessive nitrogen fertilizer inputs and poor nutrient use efficiency [11].

Recent research in tropical forests has shown that both the magnitude and direction of GHG fluxes can shift depending on topographic position, water saturation, and vegetation characteristics [23,26]. Similarly, in tropical agroecosystems, factors such as fertilizer application, tree density, canopy cover, and crop age have been shown to differentially influence CO₂, CH₄, and N₂O fluxes [12,13].

In this context, management practices grounded in agroecological principles—such as species diversification (agrobiodiversity), the use of leguminous cover crops, organic fertilization, and shade management—can help mitigate and stabilize net GHG emissions over the long term. These practices enhance carbon sequestration, improve soil biogeochemical processes, reduce nitrogen losses, and regulate the soil–atmosphere water balance [12,13,26,30]. Particularly in cocoa agroforestry systems, practices such as the integration of perennial and food-producing trees, organic fertilization, permanent soil cover, residue management, balanced canopy (shade) regulation, and appropriate irrigation and drainage systems are considered sustainable strategies that can help reduce CO₂, CH₄, and N₂O emissions, depending on the specific environmental conditions of each region [2,3,4,31,32].

In this context, in island countries such as São Tomé and Príncipe, where cocoa cultivation is primarily carried out by smallholder farmers under diverse environmental conditions and typically conventional management, the dynamics of GHG fluxes remain complex and poorly understood. Therefore, understanding how these fluxes vary in response to biophysical and management conditions is essential to developing climate-smart practices to address climate change in cocoa agroecosystems [4].

In this research, we tested the following hypotheses: (1) more biodiverse cacao systems produce lower GHG fluxes (CO₂, CH₄, and N₂O) compared to conventional systems, with differential responses depending on the type of gas; (2) GHG fluxes are influenced by the crop development stage, generally showing higher emissions during the growing stage compared to the productive stage; and (3) GHG fluxes are driven by different climatic and edaphic interactions, with gas-specific responses. In contrast, topographic factors are expected to have minimal or no direct influence on the regulation of these fluxes under the same management conditions. The objective of this study was to assess changes in CO₂, CH₄, and N₂O fluxes in cacao agroecosystems under contrasting management and crop development conditions in two climatic zones of São Tomé Island (Central Africa). Our findings contribute to identifying emission patterns and drivers across cacao agroecosystems, providing insights to guide the development of climate-smart practices and support policy guidelines aimed at achieving the NDC targets in São Tomé Island.

2. Materials and Methods

2.1. Study Sites

This study was conducted on São Tomé Island, part of the Democratic Republic of São Tomé and Príncipe, located in the Gulf of Guinea off the western equatorial coast of Central Africa (Figure 1). The island is of volcanic origin and forms part of the Cameroon Volcanic Line [33]. São Tomé presents a rugged topography with altitudinal gradients ranging from sea level to its highest peak, Pico de São Tomé, at 2,024 m above sea level. The island falls within the São Tomé, Príncipe, and Annobón moist lowland forest ecoregion and is home to a mosaic of forest and agricultural land uses [34]. A large portion of the remaining old-growth forest is protected under the Ôbo Natural Park, which spans approximately 235 km² and covers 27% of the island’s surface [35]. Outside this protected area, the landscape includes secondary forests, shaded cocoa plantations, oil palm plantations, mixed subsistence agriculture, and urban settlements, each with varying degrees of human intervention. Secondary forests are typically the result of natural regrowth after plantation abandonment, while shaded plantations combine cocoa crops with a canopy of native and exotic trees [36].

São Tomé has a humid tropical climate characterized by two distinct rainy seasons, from March to May and from October to November, and two dry periods, from June to September and from December to February [37]. Average annual temperatures range from approximately 27°C at sea level to 21°C in the highlands, with relative humidity remaining high throughout the year. Precipitation is highly variable across the island, ranging from about 900 mm in the northeastern lowlands to nearly 6,000 mm in the southern mountainous region [38,39]. The soils are derived mainly from basaltic and phonolitic volcanic rocks and are classified as deep, fertile Andosols and Nitisols, with good drainage and high organic matter content [40]. Field measurements were carried out in two climatically contrasting regions based on a regional precipitation map of São Tomé Island [39]: a dry zone (northern region, Lobata District; 0°23′–0°25′ N, 6°35′–6°37′ E; 750–1,370 mm year^-1) and a wet zone (southeastern region, Cantagalo District; 0°13′–0°15′ N, 6°39′–6°41′ E; 1,380–2,000 mm year^-1) covering a gradient of elevation (50–343 m a.s.l.) and slope (4.4°–20.6°).

2.2. Sampling Design

The campaign was carried out between April and May, 2025, coinciding with the transition from the short rainy season to the onset of the dry period on São Tomé Island. We selected eight cacao agroecosystem types defined by the combination of three factors: climatic zone (wet or dry), system type (biodiverse or conventional), and cacao development stage (growing or productive) (Figure 1). At each combination, three replicate sampling plots (20 m × 30 m) were established, resulting in a total of 24 plots. Plots were purposively selected (spaced at least 10 m apart), based on pre-existing agroecological conditions, and aimed at internal homogeneity in slope, soil type, canopy structure, and management history.

Cacao crops ages ranged from 2 to 4 years in growing plots and from 14 to 22 years in productive ones. Overall, dominant vegetation included T. cacao, Musa acuminata, Manihot esculenta, and various fruit and timber species such as Cedrela odorata, Psidium guajava, and Persea americana. Cacao trees in growing plantations exhibited diameter at 30 cm height (DH₃₀) values between 1.3 and 3.2 cm and heights from 1.4 to 2.0 m, while productive plantations showed diameter at breast height (DBH) values from 13.9 to 15.1 cm and heights from 4.2 to 4.5 m. Soil organic carbon (SOC) stocks ranged from 55.7 to 75.9 Mg ha^-1 and aboveground biomass carbon (AGBc) from 11.8 to 28.7 Mg C ha^-1, reflecting the variability in vegetation structure and organic matter across sites.

Regarding cacao agroecosystem management, biodiverse systems were characterized by the integration of more sustainable practices, including the diversification of shade trees, fruit trees, and staple food crops (>7 species), canopy management to regulate shading, the use of leguminous ground covers, and the application of organic soil amendments. In contrast, conventional systems exhibited less sustainable practices, with low species diversity (3–5 species), reduced shading levels, limited ground cover, and no implementation of organic soil management practices..

This diversity of biophysical and structural conditions provided a robust framework to evaluate the influence of climate, management, and phenology on soil GHG fluxes in tropical cacao agroecosystems (Table 1; Supplementary Figure S1).

2.3. Flux Measurements

In each sampling plot, a 5-m-radius circular subplot was established, within which three static polypropylene chambers (20 cm diameter × 20 cm height) were installed in an equilateral triangle configuration with 1 m spacing between them, adapted from Daniel et al. [26]. Chamber bases were inserted 5 cm into the soil at least 24 hours prior to gas sampling to minimize disturbance. Gas samples were collected between 9:00 a.m. and 3:00 p.m. to reduce diurnal variability in fluxes [23,41]. Four gas samples were taken from each chamber at 0, 15, 30, and 45 minutes after closure, using 12 mL pre-evacuated glass vials (Labco Ltd., UK). A 15 mL syringe was used to extract gas from the chamber headspace through a rubber septum. Chamber temperature was recorded at each interval using a digital needle thermometer (Full-Scale Traceable® Thermometer, Traceable®, USA), and soil temperature (ST) (°C) and soil water content (SWC) (%) were measured with a portable sensor (TEROS 11, Decagon Devices, USA) placed at the center of the chamber triangle. One blank vial per plot was included for ambient background control.

Samples were carefully packaged and sent to the laboratory, where CO₂, CH₄, and N₂O concentrations were analyzed using a gas chromatograph (GC-2014, Shimadzu, Japan) equipped with a flame ionization detector (FID) for CO₂ and CH₄ (with methanizer) and an electron capture detector (ECD) for N₂O. Analytical conditions were: column temperature 80 °C; FID and ECD temperatures 250 °C and 325 °C, respectively; methanizer temperature 380 °C; injection volume 2 mL (automatic loop); and nitrogen carrier gas flow of 30.83 mL/min. The limits of detection (LOD) and quantification (LOQ) were 0.0211 and 0.0704 ppm for CH₄, 37 and 125 ppm for CO₂, and 0.004 and 0.012 ppm for N₂O, respectively.

Soil GHG fluxes were calculated in R using the gasfluxes package v.0.7 [42]. The minimum detectable flux (MDF) for each gas was estimated following the procedure described by Parkin et al. [43]. For each time series of gas concentrations, four models were applied: linear regression (LR), robust linear regression (RLR), Hutchinson–Mosier (HMR), and non-steady-state diffusion estimator (NDFE) [26,44,45]. These approaches allowed assessing the sensitivity of flux estimations under varying statistical assumptions and concentration patterns. Overall, 91.67% of CO₂ fluxes were estimated using RLR and 8.33% using LR. For CH₄, 95.83% were calculated with RLR and 4.17% with LR. In the case of N₂O, 90.28% of fluxes were estimated using RLR and 9.72% with LR. Gas concentrations (ppm) were converted to mass using the ideal gas law and normalized by the surface area of each chamber. Fluxes below the detection limit were excluded from analysis: 5.56% for CO₂, 12.5% for CH₄, and 15.28% for N₂O [43].

The global warming potential (GWP) at factorial combination was estimated by projecting the soil CO₂, CH₄, and N₂O fluxes measured during the sampling campaign to annual values. GWP was expressed in t CO₂-eq ha^-1 year^-1 using the global warming potential conversion factors of 27.2 for CH₄ and 273 for N₂O, as reported by the IPCC [46].

2.4. Environmental Variables Measurements

After the final GHG measurements, nine soil subsamples (one per chamber across three plots) were collected from the chamber bases at a depth of 0–20 cm. These were then homogenized to obtain three composite samples per factorial combination, following the approach of Pang et al. [47]. Soil samples were processed in the laboratory for subsequent physicochemical analyses. Soil texture (sand, clay, and silt) (%) was determined using the Bouyoucos hydrometer method. Bulk density (BD, g cm⁻³) was measured with an Eijkelkamp hand auger. Soil pH and electrical conductivity (EC, dS m⁻¹) were assessed using the saturation paste and the conductometric method [48]. Cation exchange capacity (CEC, meq 100 g⁻¹) was determined by titration with 1 M NaOH, and base saturation (BS, %) was calculated as the sum of exchangeable base cations (Ca²⁺, Mg²⁺, K⁺, and Na⁺), divided by CEC [48]. Soil organic carbon (SOC, Mg ha⁻¹) was estimated from the oxidizable organic carbon (COOx) content, obtained using the potassium dichromate colorimetric method (NTC 5403, Walkley–Black), and converted following the methods of Nelson and Somers [49] and Batjes [50]. Available phosphorus (P, mg kg⁻¹) was determined using the Bray II method, and total nitrogen (N, %) content via the Kjeldahl method [48].

Climatic and topographic variables were processed per plot using R v.4.4.3 [51] and RStudio v.2025.05.04 [52]. Elevation (m a.s.l.) was derived from AWS Terrain Tiles via the elevatr package [53]. A digital elevation model (DEM, ~30 m resolution) of São Tomé was downloaded using elevatr, from which slope (degrees) was derived with the terrain() function from the terra package [54], and values were extracted per plot using sf and terra. Daily weather data were then extracted from the NASA POWER platform using the get_power() function from the nasapower package [55], based on plot coordinates and sampling dates. Variables included air temperature at 2 m (Temp, °C), relative humidity at 2 m (RH, %), precipitation (Prec, mm day^-1), 5-day cumulative precipitation (Prec_5d, mm), and solar radiation (RS, MJ m⁻² day⁻¹). Finally, reference evapotranspiration (ET₀, mm day^-1), was calculated using a simplified version of the Hargreaves equation [56].

2.5. Statistical Analysis

Linear mixed-effects models (LMMs) were fitted to assess the fixed effects of zone climatic, system type, development stage and their interactions on GHG fluxes using the lme function from the nlme package [57] from R 4.4.3 [51], implemented through the InfoStat v.2020 interface [58]. Plot was included as a random effect (1|Plot). Normality and homoscedasticity were evaluated thorough exploratory analysis of the model residuals. CO₂ flux data were log-transformed, whereas CH₄ and N₂O fluxes were transformed using the formula log(flux − min(flux) + 1) in order to meet the assumption of normality [23]. To address heteroscedasticity across system types, a variance structure (varIdent) was applied [57]. Post hoc comparisons of means across fixed effects were performed using Fisher’s LSD test (α=0.05).

To explore the strength of association between gas fluxes and environmental factors (i.e., soil properties and aboveground conditions), pairwise Spearman correlation matrices were computed using the rcorr function from the Hmisc R package [59]. Based on significant correlations, dependency relationships between gas fluxes and environmental variables were further examined through simple linear regression using the lm function from the stats package [60] in R. Model evaluation included the coefficients of determination (R²) and residual diagnostics to assess the quality of fit. Resulting models were visualized using the ggplot2 R package [61].

Finally, structural equation models (SEMs) were used to examine the direct and indirect effects of environmental variables on each gas flux. The best SEMs were selected based on the Chi-square test (χ²), degrees of freedom (df), Comparative Fit Index (CFI), and Standardized Root Mean Square Residual (SRMR). SEMs were performed using the sem function from the lavaan package [62], and visualized with the grViz function from the DiagrammeR package [63] in R. All statistical analyses in R were conducted using the RStudio v.2025.05.0 interface [52].

3. Results

3.1. Soil CO₂, CH₄, and N₂O Fluxes

Significant effects on soil CO₂ fluxes were observed for the highest-order interaction among climatic zone, system type, and cacao development stage (p < 0.05) (Figure 2A,B; Supplementary Table S1), as well as for the interaction between system type and development stage (Supplementary Figure S2). Overall, the highest CO₂ fluxes were recorded in conventional growing systems located in the wet zone (125.5 mg C m⁻² h⁻¹), whereas the lowest emissions were observed in biodiverse productive systems within the same zone (28.3 mg C m⁻² h⁻¹) (Figure 2A,B; Supplementary Table S2). Averaging across climatic zones, CO₂ fluxes were significantly higher in conventional systems at the growing stage (100.77 ± 8.07 mg C m⁻² h⁻¹) than in biodiverse systems at the productive stage (37.09 ± 8.07 mg C m⁻² h⁻¹) (Supplementary Figure S2A). A similar pattern was observed between the two system types during the productive stage, with higher CO₂ fluxes in conventional systems (53.90 ± 4.00 mg C m⁻² h⁻¹) compared to biodiverse systems (37.09 ± 8.07 mg C m⁻² h⁻¹) (p < 0.05). During the growing stage, CO₂ fluxes were also significantly higher in conventional systems compared to biodiverse systems.

For CH₄, no significant effects on soil CH₄ fluxes were observed for the highest-order interaction (Figure 2C,D; Supplementary Table S1). However, a significant interaction between system type and development stage was evidenced (p < 0.01) (Supplementary Figure S2B). Biodiverse–growing systems exhibited the highest mean CH₄ emissions (5.08 ± 1.50 μg C m⁻² h⁻¹), whereas conventional–growing systems showed the strongest CH₄ uptake (−12.58 ± 2.35 μg C m⁻² h⁻¹) (Supplementary Figure S2B; Supplementary Table S2). During the productive stage, no significant differences were observed between the two system types (p > 0.05).

Regarding N₂O fluxes, the highest-order interaction showed no significant effect (Figure 2E,F; Supplementary Table S1). Only the interactions between system type and development stage, as well as between climatic zone and development stage, showed significant effects on mean N₂O fluxes (both, p < 0.05) (Supplementary Figures S2C and S3). The highest N₂O fluxes were recorded in conventional growing systems (32.28 ± 5.76 μg N m⁻² h⁻¹), whereas biodiverse growing systems exhibited the lowest mean values (−12.59 ± 4.63 μg N m⁻² h⁻¹) (Supplementary Figure S2C). This pattern was further supported by the interaction between climatic zone and development stage, with the highest N₂O emissions observed during the growing stage in the dry zone (22.00 ± 5.22 μg N m⁻² h⁻¹), compared to the productive stage in the same zone (−6.25 ± 5.22 μg N m⁻² h⁻¹) (Supplementary Figure S3). During the growing stage, no significant differences were observed between the two zones (p > 0.05).

Figure 3 summarizes the differences and sums of means for each GHG across main effects. CO₂ flux showed the greatest contrast among the difference values with higher emissions at the growing stage (+37.3 mg C m⁻² h⁻¹) and in the dry zone (+11.3), and lower in biodiverse systems (-26.4) (Figure 3A). CH₄ emissions were slightly higher in biodiverse systems (+8.5) but lower in the dry zone and growing stage (−4.8 and −4.7, respectively). For N₂O, fluxes were higher in the growing stage (+13.4), dry zone (+9.5), and lower in biodiverse systems (−26.9).

Regarding the sum of fluxes (Figure 3B), CO₂ showed a greater contribution of positive fluxes from the growing stage (64.5%) compared to the productive stage. Likewise, conventional systems contributed more (60.3%) than biodiverse systems, and the dry zone contributed more (54.4%) than the wet zone. For CH₄, the largest contributions of negative fluxes corresponded to the growing stage (79.2%), conventional systems (66.2%), and the dry zone (78.5%). In the case of N₂O, although more balanced, positive fluxes remained dominant in the growing stage (73.5%), conventional systems (61.6%), and the dry zone (83.2%).

3.2. Soil Temperature and Water Content

During the field campaign, significant differences in ST were observed for the highest-order interaction (Figure 4A,B; Supplementary Table S1). Overall, the highest ST values were recorded in biodiverse systems at the productive stage in the wet zone (31.09 °C), whereas the lowest values were observed in biodiverse systems at the same stage in the dry zone (26.60 °C).

For SWC, no significant differences were observed for any of the studied effects or their interactions (Figure 4C,D; Supplementary Table S1). However, some trends were observed, with higher mean values in conventional productive systems in the wet zone (35.94%) and lower values in biodiverse growing systems from the same zone (22.29%) (Figure 4C,D)

3.3. Environmental Drivers of Soil GHG Fluxes

The analysis of climatic, topographic, and soil variables revealed contrasting patterns among climatic zones, system types, and developmental stages (Supplementary Table S3). Climatically, Prec_5d and Prec were consistently higher in the wet zone than in the dry zone for both stages, with particularly high values in conventional systems during the growing stage (Prec_5d: 23.14 ± 5.89 mm; Prec: 7.19 ± 2.14 mm). ET₀ and SR were greater in the dry zone, especially during the growing stage (ET₀: 402.06 ± 0.01 mm; SR: 22.28 ± 0.01 MJ m⁻² d⁻¹). RH reached higher values in the wet zone (>87%), whereas Temp showed minimal variation (27.2–27.5 °C). Slope showed notable variation between zones and systems. In the growing stage, conventional systems in the wet zone recorded the steepest slopes (20.59 ± 0.01°), whereas the lowest slopes were observed in the dry zone (4.44–6.36°). In the productive stage, slopes were moderate (10.95–12.06°). Elevation presented a marked contrast: systems in the dry zone were located at low altitudes (50–214 m a.s.l.), while those in the wet zone reached higher elevations, up to 343 m a.s.l. in conventional systems during the growing stage. Among soil properties, Clay was higher in the wet zone (up to 21%), whereas Sand predominated in the dry zone (>60% in some cases). Biodiverse systems in the dry zone had higher SOC values (up to 75.87 ± 0.75 Mg C ha⁻¹) and N than conventional systems. BD was greater in conventional systems (up to 1.19 g cm⁻³). Soil pH ranged from 5.35 to 5.77, with slightly higher values in the dry zone. CEC was higher in biodiverse systems in the wet zone (12.10 ± 1.01 meq 100 g⁻¹).

Spearman’s correlation analysis revealed multiple significant relationships between CO₂, CH₄, and N₂O fluxes and various climatic and soil variables (Figure 5).

CO₂ flux showed significant positive correlations with Prec, Temp, Prec_5d, and ST, and significant negative correlations with RH, SR, and Clay. Simple linear regression analysis of CO₂ flux against each of these variables indicated that only the relationship with Clay was not significant (Supplementary Figure S4). CH₄ flux exhibited significant negative correlations with Prec_5d, Prec, and Temp, and positive correlations with RH, SR, ET₀, pH, and SWC. However, the dependencies of CH₄ flux on SR and ET₀ were not significant (Supplementary Figure S5). N₂O flux was positively correlated with Prec, BD, N, SOC, CEC, EC, and BS, with no significant negative correlations detected with environmental variables. The dependencies of N₂O flux on CEC and BS were not significant (Supplementary Figure S6).

The SEM results explained 63% of the variance in CO₂ flux, 32% in SOC, and 14% in ST (Figure 6A). The overall model fit was very good (χ² = 4.72, df = 6, p > 0.05; CFI = 1.000; SRMR = 0.080), indicating that the specified structure was consistent with the observed data. ST had a positive and significant effect on CO₂ flux (β = 0.59, p < 0.001). SR showed a significant negative effect on CO₂ flux (β = –0.47, p < 0.01). BD had a positive and significant impact on SOC (β = 0.52, p < 0.05).

The SEM explained 58% of the variance in CH₄ flux, 24% in ST, and 22% in soil pH (Figure 6B). The model fit was good (χ² = 8.15, df = 7, p > 0.05; CFI = 0.954; SRMR = 0.040). Regarding direct effects, CH₄ flux was positively influenced by soil pH (β = 0.38, p < 0.05). In contrast, Prec had a significant negative effect on CH₄ flux (β = –0.47, p < 0.05), as did Temp (β = –0.39, p < 0.05). Temp indirectly affected CH₄ flux by increasing soil pH (β = 0.49, p < 0.05). Additionally, RH had a positive effect on soil pH (β = 0.64, p < 0.05) but a negative effect on ST (β = –1.63, p < 0.05). Conversely, ST was significantly increased by SR (β = 1.02, p < 0.05) and decreased by RH (β = –1.63, p < 0.05).

The SEM for N₂O flux showed a very good fit (χ² = 1.736, df = 2, p > 0.05; CFI = 1.00; SRMR = 0.033) (Figure 6C). The model explained 69% of the variance in N₂O flux, 37% in EC, and 17% in ST. SOC was the strongest and most significant predictor of N₂O flux (β = 0.67, p < 0.001). ST also had a significant positive effect (β = 0.30, p < 0.05). A negative effect of pH on N₂O emissions was observed (β = –0.36, p < 0.01). P had a significant positive effect on EC (β = 0.58, p < 0.001). Higher Temp values were positively associated with increased ST (β = 0.44, p < 0.05).

3.4. Global Warming Potential (GWP)

The global warming potential (GWP) of soil GHG fluxes varied markedly across system types and development stages (Table 2). In both climatic zones, higher GWP values were observed in conventional systems compared to biodiverse systems. Notably, in the wet zone, GWP reached 40.90 Mg CO₂ eq ha^-1 year^-1 in conventional systems during the growing stage, over three times higher than in the biodiverse counterpart (12.98 Mg CO₂ eq ha^-1 year^-1). Conversely, the lowest GWP values were observed in biodiverse systems during the productive stage in the wet zone (9.05 Mg CO₂ eq ha^-1 year^-1). In the dry zone, differences were less pronounced, though conventional systems still tended to emit more than biodiverse systems at the same development stage. Overall, the relative contribution of CH₄ and N₂O fluxes to total GWP was low, accounting for approximately 2.58%, with CO₂ being the dominant component of soil-derived climate impact.

4. Discussion

4.1. CO₂ Fluxes are Higher in Conventional Systems and During the Growing Stage

Our results show that CO₂ fluxes in cacao agroecosystems are strongly influenced by the interaction between climatic zone, system type, and development stage (Figure 2 and Figure 3). The highest soil CO₂ emissions were observed in conventional systems during the growing stage in the wet zone, reaching 125.5 mg C m^-2 h^-1, whereas the lowest values were evidenced in biodiverse systems during the productive stage within the same zone (28.3 mg C m⁻² h⁻¹). The highest CO₂ flux recorded in the conventional system was more than three times higher than that observed in the same system type during the productive stage (37.09 mg C m⁻² h⁻¹), highlighting the influence of management intensity and phenological stage on soil respiration. These findings align with the seasonal trends reported by Pang et al. [12], where soil CO₂ emissions in temperate forests were consistently higher during the growing season due to increased root activity and microbial metabolism, both strongly driven by temperature and substrate availability [64]. These results align with previous studies reporting that the more intensive management practices typical of conventional systems, combined with lower canopy complexity, increase soil CO₂ emissions by raising soil temperature and altering soil moisture dynamics [6,65,66].

The correlation and pathway patterns observed in our study suggest that CO₂ fluxes are strongly regulated by both climatic and edaphic factors, particularly precipitation, relative humidity, solar radiation, air temperature, clay content, and soil temperature (Figure 5 and Figure 6). This increase in CO₂ fluxes is due to greater soil respiration under warm and humid conditions that favor microbial and root activity [67,68]. The negative correlation with clay content could be associated with reduced aeration and CO₂ diffusion in finer-textured soils [69]. Similar to the observations of Pang et al. [12], our results suggest that finer textures and higher soil moisture levels may constrain gas diffusivity and alter the balance between aerobic and anaerobic processes. The negative relationship observed between solar radiation and soil CO₂ flux is explained by the greater presence of shade trees in cacao agroforestry systems, which substantially reduce the incident radiation on the soil. This attenuation of radiation creates more favorable microclimatic conditions—such as higher soil moisture and thermal stability—that stimulate soil respiration and, consequently, increase CO₂ flux [70]. In our study, neither elevation nor slope had a significant association with CO₂ fluxes. These results align with Daniel et al. [26], where the topographic gradient did not affect CO₂ emissions, but contrast with Pang et al. [12], who argue that site elevation can be used as a predictor of regional CO₂ emissions.

4.2. Net CH₄ Emissions in Biodiverse Systems and Greater Uptake in Conventional Systems

CH₄ fluxes in cacao agroecosystems were strongly influenced by the interaction between system type and development stage (Figure 2 and Figure 3). Biodiverse growing systems exhibited net CH₄ emissions (5.08 µg C m⁻² h⁻¹), whereas conventional growing systems acted as stronger sinks (−12.58 µg C m⁻² h⁻¹). These results indicate that CH₄ flux direction may shift depending on vegetation structure, soil moisture, and organic inputs, as elevated rainfall and litter-derived labile carbon can enhance methanogenesis under temporarily saturated conditions [47,71]. Although biodiverse cacao systems in the growth stage exhibited net CH₄ emissions, their greater capacity to store carbon in biomass and soil (Table 1) offsets the impact of these fluxes on the system’s net carbon balance. Systems with higher tree diversity and canopy cover not only reduce incident radiation and stabilize soil temperature and moisture but also enhance long-term carbon stocks [72,73]. A recent study reported that diversified agroforestry systems have a significantly lower carbon footprint (932.1 ± 251.6 kg CO₂e ha⁻¹) compared to traditional monoculture systems (1914.4 kg CO₂e ha⁻¹) [72].

In terms of environmental drivers of CH₄ fluxes (Figure 5 and Figure 6), soil pH had a positive effect, indicating that less acidic conditions favor methane production, as most methanogens have an optimal pH range close to neutrality [74]. Conversely, precipitation and air temperature showed negative effects, suggesting that warmer and wetter conditions may reduce methane emissions by limiting anaerobic conditions in the soil. In well-aerated soils, increases in surface moisture associated with recent rainfall can stimulate methanotrophic oxidation, thereby lowering net emissions. This is consistent with the positive correlation observed between CH₄ flux and soil water content (Figure 5). These patterns reflect the dual nature of CH₄ as both a product of anaerobic processes and a substrate for aerobic oxidation [75,76].

Air temperature also indirectly influenced CH₄ fluxes by increasing soil pH, possibly through effects on base availability or organic matter decomposition dynamics [76]. Relative humidity had a positive effect on soil pH but a negative effect on soil temperature, reinforcing its role as a microclimatic regulator and suggesting an indirect positive influence on CH₄ emissions. In our study, none of the topographic variables significantly affected methane fluxes. This contrasts with Pang et al. [12], who reported higher methane emissions at sites with greater elevation, attributed to increased rainfall.

Overall, our findings support the view that CH₄ dynamics in tropical cacao agroecosystems arise from complex interactions between climatic factors (air temperature, relative humidity, and precipitation) and edaphic properties (particularly soil pH), modulated by management practices such as tree diversity.

4.3. Higher N₂O Emissions in Conventional Systems and During Growing Stages

For N₂O fluxes, the significant interactions observed between system type and development stage, as well as between climatic zone and stage, highlight the complexity of the controlling factors (Figure 2 and Figure 3). Growth-stage conventional systems recorded the highest fluxes (32.28 µg N m⁻² h⁻¹), whereas biodiverse systems functioned as net sinks (−12.59 µg N m⁻² h⁻¹). In the dry zone, emissions were markedly higher during the growth stage compared to the productive stage. Higher emissions in conventional systems may be associated with more intense nitrification–denitrification cycles resulting from soil disturbance and shorter fallow periods. In particular, moisture events following dry spells trigger substantial emission peaks, suggesting that the absence of vegetative cover and conventional management practices contribute significantly to nitrous oxide emissions [77,78]. Similarly, a recent study conducted in cacao plantations in West Africa found that shaded agroforestry systems (i.e., combining cocoa with shade or other fruit or trees and food crops) exhibit lower N₂O fluxes compared to conventional open-sun systems, owing to their greater tree cover, complex microbiotic interactions, and improved soil quality [3].

Regarding environmental drivers, soil organic carbon had the strongest and most significant direct positive effect on N₂O flux (Figure 5 and Figure 6). This finding is consistent with previous studies showing that soils with high organic matter content promote heterotrophic microbial activity and denitrification, a key process in nitrous oxide production [79]. Conversely, Pang et al. [12], reported lower N₂O emissions in sites with high dissolved organic carbon (DOC) and ammonium availability, where microbial communities may facilitate the complete reduction of N₂O under moist, low-oxygen conditions.

Soil temperature also had a direct positive effect, whereas air temperature exerted an indirect positive effect on N₂O fluxes, indicating that higher temperatures enhance microclimatic energy flows and, consequently, respiration rates and microbial processes such as nitrification and denitrification [80,81]. Soil pH had a direct negative effect on N₂O emissions, suggesting that acidity constrains the final reduction of N₂O to N₂ during denitrification, thereby favoring N₂O accumulation and release [82].

In our study, elevation and slope showed no significant impacts on N₂O fluxes, suggesting that topographic effects may occur but are context-dependent [83]. However, Pang et al. [12], observed increased soil N₂O fluxes at higher elevations, driven by greater precipitation and more pronounced freeze–thaw cycles during winter.

4.4. Global Warming Potential (GWP) and the Role of Climate-Smart Practices

The global warming potential (GWP) of soil GHG fluxes varied markedly across zone climatic, system type and cacao development stage, with the highest values observed in conventional growing systems in the wet zone, reaching 40.90 Mg CO₂ eq ha^-1 year^-1. In contrast, the lowest GWP was recorded in biodiverse systems at the productive stage in the wet zone (9.05 Mg CO₂ eq ha^-1 year^-1), indicating that both climate and management factors interact to influence the net climate impact of cacao agroecosystems.

On average, CO₂ contributed 97.42% of total GWP, with CH₄ and N₂O together accounting for only 2.58%. This dominance of CO₂ aligns with the findings of Pang et al. [12], who also reported that soil respiration is the primary driver of GWP in forested systems, especially under conditions of high organic matter turnover and biological activity. The limited contribution of CH₄ and N₂O reinforces the idea that, although these gases have much higher warming potentials per unit mass, their relatively low flux rates in upland tropical soils result in marginal climate impacts when compared to CO₂.

The strong differences observed between conventional and biodiverse systems, especially in the wet zone, are likely associated with management-induced changes in soil structure, organic matter inputs, and microbial activity. As previously shown in tropical and subtropical systems, management practices that enhance litter cover, reduce disturbance, and promote root diversity can stabilize soil respiration and reduce net GHG emissions [12,26]. Moreover, increased shading and surface cover in biodiverse plots may moderate soil temperature and moisture fluctuations, further buffering CO₂ release [3,68].

The cacao developmental stage had a notable influence on GWP values, with growing-phase crops consistently exhibiting higher emissions than productive-phase counterparts. This trend likely reflects the physiological and microbial dynamics associated with early vegetative growth, including intensified root respiration, greater microbial decomposition of recently incorporated organic matter, and favorable moisture conditions [76]. Similar observations have been reported in tropical agroforestry systems, where initial stages of plant establishment often coincide with peaks in biological activity and CO₂ efflux [12,23]. In cacao systems, these conditions are further amplified by management practices such as organic fertilization, pruning, or replanting, which may temporarily stimulate soil respiration rates [80]. Moreover, less-developed canopies in growing plots allow greater radiation penetration and soil warming, potentially enhancing microbial metabolism and CO₂ release [3,65]. These findings emphasize the importance of considering phenological stage when assessing the climate impact of agroecosystems, as GHG emissions may vary substantially throughout the crop cycle. Incorporating stage-specific dynamics into monitoring and mitigation frameworks could therefore improve the accuracy of emission estimates and inform the timing of targeted interventions in tropical cacao landscapes.

5. Conclusions

This study revealed that the highest CO₂ fluxes occurred in conventional cacao systems during the growing stage under locally wetter conditions (wet zone), whereas biodiverse systems consistently exhibited lower emissions. CH₄ fluxes generally acted as a sink, particularly in conventional systems located in the dry zone. However, net CH₄ emissions were recorded in biodiverse systems in the wet zone, suggesting that local management practices and site-specific conditions can reverse the natural sink function of soils. N₂O fluxes were generally low or negative and showed high variability across cacao agroecosystems, with the highest emissions observed in conventional growing systems in the dry zone.

Our results demonstrates that GHG dynamics in tropical cacao agroecosystems are primarily governed by intricate interactions between climatic factors—such as air temperature, relative humidity, and precipitation—and key edaphic properties, including soil pH, soil organic carbon, soil temperature, and clay content, while local topographic variables (elevation and slope) showed no significant influence.

In this study, CO₂ fluxes were the dominant driver of the GWP, accounting for approximately 97.6% of total emissions expressed as t CO₂-eq ha⁻¹ year⁻¹. The highest annual GWP was recorded in conventional cacao systems during the growing stage in wet environments (40.90 t CO₂-eq ha⁻¹ year⁻¹), whereas the lowest occurred in biodiverse productive systems under the same local conditions (9.05 t CO₂-eq ha⁻¹ year⁻¹). Although CH₄ and N₂O contributed less than 3% to the total GWP, their individual dynamics remain relevant due to their high global warming potential per unit mass.

Overall, our findings highlight the substantial climate mitigation potential of biodiverse cacao agroecosystems (i.e., biodiversity-based agroforestry), particularly when integrated with other agroecological management practices (e.g. canopy regulation, organic fertilization, and the use of legume cover crops). The promotion of such climate-smart practices could significantly reduce net GHG emissions from tropical cacao plantations, especially in wet environments and during early crop developmental stages. Future research should focus on year-round monitoring of these emissions and the integration of soil microbiological assessments to improve the accuracy of GHG balance estimates and guide the development of effective mitigation strategies in tropical cacao agroforestry systems.