The study was conducted from June 10, 2024, to January 27, 2025, at The Starbucks Coffee Company’s Alsacia Farm in Alajuela, Costa Rica (10°06’01.3” N, 84°11’41.8” W). The soil at the experimental plot, was classified as an Andisol with an ustic moisture regime, characterized by high organic matter content (> 5%) and low fertility (CEC< 952 mg kg^-1).

Over the crop cycle, temperatures ranged from 13°C to 19°C, with a mean of 16.8°C, and total accumulated precipitation reached 3900 mm. Meteorological data were recorded by an automatic weather station located in Laguna Fraijanes, 17 km from the experimental site, and provided by the National Meteorological Institute under contract IMN-DIM-CM-179-2025. Between August 16, 2024, and January 27, 2025, measurements of Water-Filled Pore Space (WFPS) in the soil were taken using an MP406 sensor connected to an ICT Bluetooth interface around the gas sampling chamber, during each gas flux sampling day. The WFPS ranged from 35% to 55%, with no significant differences in mean values observed among the three fertilization treatments or sectors (Supplementary Table 1).

To characterize the physical and chemical properties of the soil, samples were sent to the Soil Laboratory of the Agronomic Research Center (CIA), University of Costa Rica. Soil texture was determined using the Bouyoucos hydrometer method. Bulk soil density was measured from undisturbed samples collected at a depth of 10 cm, using 250 cm³ soil cores at six locations (one per block).

A composite soil sample from the 20 cm topsoil was collected for fertility analysis from Sector A and Sector B (Figure 1A). Samples were taken from 45 points per sector on July 13, 2024, and at the end of the experiment on January 27, 2025. The Total carbon (C) and nitrogen (N) contents were analyzed through dry combustion, while soil pH was determined in a 1:2.5 soil-to-water suspension. The concentrations of calcium (Ca), magnesium (Mg), potassium (K), and total phosphorus (P) were obtained using a modified KCl-Olsen extraction method. Trace elements including zinc (Zn), copper (Cu), iron (Fe), and manganese (Mn) were measured by atomic absorption spectroscopy (novAA 400p, Analytik Jena, Germany).

The experiment was conducted using a randomized block design with three treatments and six repetitions. The blocks were distributed into two sectors based on topography (Sector A and Sector B) (Figure 1A), each with distinct soil textural class (Table 1). A total of 18 plots, each measuring 15 m × 16 m, were established. The plots were planted with the “Caturra” coffee variety at a spacing of 0.90 m between plants and 1.80 m between rows, resulting in a planting density of 6, 216 plants ha^-1. Adjacent plots were separated by four coffee rows to minimize edge effects. Additionally, the initial foliar chemical status was assessed using composite samples consisting of 90 leaves per sector, collected on July 20, 2024 (Supplementary Table 2), to evaluate the nutritional status of the plants in both sectors and to ensure comparable baseline conditions.

Details of the experimental treatments, the total N applied in each fertilization strategy, and the corresponding application schedule are presented in Table 2. The total N rate for F and Y treatments (346 kg N ha^-1) was divided as follows: flowering stage (100 g plant^-1), grain filling (115 g plant^-1) and ripening (112 g plant^-1). In the case of U treatment, 414 kg ha^-1 were divided into two stages floral differentiation (100 g plant^-1) and grain filling (115 g plant^-1).

Coffee yield was evaluated across four harvest dates (November 19, December 6, December 26, 2024; and January 28, 2025) within the 100 m² productive plot of each block. On each date, only ripe cherries were harvested. Bean classification was performed according to Starbucks quality standards, excluding green beans, floaters, dried beans, and those damaged by the coffee berry borer. For each experimental unit, the total fresh weight of harvested coffee cherries was recorded and normalized by the number of productive plants.

Insect pests and disease control were carried out through regular field monitoring. To support floral differentiation and improve plant health, boric acid, magnesium sulfate, potassium nitrate, and a phytoprotectant containing calcium, silicon, and copper were applied.

Weed management was conducted using chemical control, complemented by the application of a selective systemic herbicide for both pre- and post-emergence control (a.i. clomazone). Areas infested with weeds that were not chemically treated were managed mechanically using a brushcutter, while climbing plants were removed manually.

Management of coffee leaf rust (Hemileia vastatrix Berk. & Br.) involved the application of fungicides with different modes of action to prevent resistance development. These included quinone outside inhibitors (a.i. strobilurins), ergosterol biosynthesis inhibitors (a.i. triazoles), spore germination inhibitors (a.i. pyraclostrobin), fungicides that disrupt neurotransmitter synthesis (a.i. dithiocarbamates), and copper-based compounds. Control of the coffee berry borer (Hypothenemus hampei) was carried out using pyrethroids, including the active ingredient dimethoate, which inhibits acetylcholinesterase.

Emissions of CO₂, CH₄, and N₂O were monitored using the non-steady-state chamber method. Frames were set 15 cm below the soil surface, two weeks prior to the start of the sampling period. Two bases were installed in each experimental unit to allow alternating chamber placement across sampling dates. Each chamber, with a surface area of 0.16 m² and a height of 0.1 m, was covered with insulating material to minimize temperature fluctuations. A small internal fan ensured gas homogenization throughout the measurement period. A flexible gas bag was connected to each chamber to maintain equilibrium between the inside and outside pressure.

A total of 33 sampling rounds were conducted. Sampling was performed on day 0 (prior to fertilizer application) and on days 1, 2, 4, 7, 8, 11, and 15 after each fertilization. Thereafter, sampling was carried out every two weeks until the next fertilization event or harvest. Gas samples were collected at 0, 20, and 40 min after chamber closure and injected into vacuum vials, for a total of 54 samples per day. In addition, two air samples, two recovery controls, and one blank were analyzed on each sampling day for quality data control.

Gas samples were analyzed using a gas chromatograph (Agilent 7890A, USA) equipped with a methanizer, a flame ionization detector (FID), and an electron capture detector (ECD), coupled to an autosampler (Agilent 7697A, USA). The instrument was calibrated for each target analyte using a four-point calibration curve, performed at the start of the analysis, after every 24 h of continuous operation, and upon completion of the injection sequence. Standard gases prepared with air as the balance exhibited an uncertainty of approximately 5%. To verify equipment performance, a certified gas mixture containing CO₂, CH₄, and N₂O was analyzed after every set of 14 vials. The observed deviations from the certified concentrations did not exceed 20%.

The hourly concentration changes for CH₄ and CO₂ fluxes was obtained from a linear regression while a quadratic regression was used for N₂O. The flux of each GHG (f_i) was determined according to Equation 1.

where atmospheric pressure (P, in Pa) and headspace temperature (T, in °C) were recorded during each measurement round using a Kestrel 4000 weather meters (Loftopia LLC., MI, USA). The molar mass (MMi) was set at 12 μg C μmol^-¹ for CO₂ and CH₄, and 28 μg N μmol^-¹ for N₂O. The ideal gas constant (R) was taken as 8.314 J K^-¹ mol^-¹ (Dawar et al., 2021). Chamber height (h, in meters) was adjusted by accounting for the distance between the frame and the soil surface.

Data analysis of fluxes was conducted under the following criteria: CO₂, CH₄ and N₂O fluxes were rejected if the coefficient of determination (r²) of the linear regression for the change in CO₂ concentration was less than 0.88, values below this threshold were considered indicative of potential leakage. For N₂O and CH₄, the r² acceptance thresholds were 0.85 and 0.70, respectively. To ensure data quality, 4% of the CO₂ and 9% of the N₂O flux readings were removed from the dataset, while 26% of the CH₄ fluxes were excluded due to low r² values. To estimate the missing values, the average slope of the replicates for each treatment was used according to the corresponding measurement day.

Cumulative CH₄ and N₂O emissions were calculated using the trapezoidal method (Pérez-Castillo et al., 2021), applying the median flux for each specific treatment and sampling day when a flux value was excluded based on the criteria described in the previous section.

Nevertheless, since fluxes could not be measured for an entire year, the cumulative emissions of N₂O (g N₂O kg of N_applied^-1 ha^-1) were normalized by total amount of N applied for each fertilization strategy for comparative purposes. The yield-scaled emissions (g CO₂e kg_coffee^-¹) were estimated by converting cumulative CH₄ and N₂O fluxes into CO₂ equivalents, using conversion factors of 27 for CH₄ and 273 for N₂O, according to IPCC guidelines (IPCC, 2021). These values were then divided by the total kilograms of coffee produced per hectare and per treatment.

Normality and homogeneity of variances were assessed using the Shapiro–Wilk test and Bartlett’s test, respectively. A one-way ANOVA was conducted, considering blocks A and B as an experimental error factor, to identify significant differences among treatments. When significant differences were detected by ANOVA, Tukey’s Honest Significant Difference (HSD) test was applied for multiple comparisons of treatment means. Statistical significance was set at p < 0.05. All statistical analyses were performed using RStudio software, version 4.5.1. To identify the variables most strongly associated with coffee productivity, we applied an information-theoretic model selection approach based on the Akaike Information Criterion (AIC). The drivers of yield variation were constructed using key soil properties and gas fluxes and management factors. For each model, AIC and Akaike weights were used to assess the relative likelihood of each model. Full model comparisons are provided in Supplementary Table 3.

Coffee yield did not show statistically significant differences between the evaluated N fertilization strategies. In contrast, significant differences in productivity were observed between sectors (p = 0.0026). Sector A exhibited higher mean yields, ranging from 5, 060 to 7, 337 kg ha^-¹, while in Sector B yields ranged from 2, 530 to 3, 036 kg ha^-¹ (Figure 1B).

Soil texture in Sector A contained 50% more clay and 47% more silt than in Sector B, resulting in a different textural class between the two sectors. The mean soil bulk density across sectors was 0.74 g cm^-³ (Table 1). The initial soil chemical conditions were similar between sectors, except for iron (Fe), element with higher content in Sector A (993 ± 23 mg kg^-¹) than in Sector B (883 ± 40 mg kg^-¹) (Table 3). At the end of the crop cycle, available phosphorus (P) in Sector A increased significantly by nearly threefold compared to its initial value (p = 0.026), whereas P levels in Sector B remained unchanged. Final Fe content decreased in Sector B compared to its initial condition (p = 0.022), whereas in Sector A it remained unchanged.

Daily fluxes of CO₂ and N₂O increased following the three N fertilization events (June 10, August 12, and October 14) and their associated management strategies (Figures 2A–C). Daily CO₂ fluxes ranged from 42 to 123 mg CO₂ m^-² h^-¹. CO₂ fluxes were influenced by the phenological development of the coffee crop (Figure 2A), they were highest during the flowering stage of the coffee plants and decrease during the grain filling and ripening stages. CH₄ fluxes consistently show negative values, regardless of the N fertilization strategy applied (Figure 2B).

Daily N₂O emissions showed similar temporal patterns throughout the entire measurement period for treatments F and Y, with no significant differences between the two strategies. Several notable N₂O emission peaks were observed following both the first and second fertilization events, closely associated with episodes of intense rainfall followed by short drying periods (Figures 2C–D, Supplementary Table 4). The highest N₂O emission peak across all treatments occurred after the second fertilization on August 12, coinciding with two substantial rainfall events on August 3 (55 mm) and August 12 (71 mm). A secondary emission peak was recorded on August 14, two days after fertilization, aligning with a short dry spell (August 14: 0.22 mm).

Cumulative CO₂ fluxes showed no significant differences among fertilization events or across N fertilization strategies. In the case of CH₄, all treatments displayed negative cumulative fluxes (Figure 3B), which became more pronounced on January 26.

Cumulative N₂O fluxes differed among fertilization treatments. When urea was applied (first and second fertilization events), it generated the highest emissions (Figure 3C). During the first event, Y produced N₂O fluxes 3.8 times higher than those of F, whereas during the second and third events, emissions were 2.3 and 1.8 times higher in F than in Y, respectively. After the second fertilization, when all treatments received 122 kg N ha^-1, N₂O emissions from Y and F were 38% and 42% lower than those from U over the period August 12–September 26. Overall, cumulative N₂O emissions were significantly higher in treatment U than in F and Y, with no significant differences between the latter two (Figure 3C).

The comparison of N₂O cumulative emission was carried out between treatments F and Y, since both received the same total N rate over the entire crop cycle. No statistically significant differences in the N₂O emission were observed between F and Y (Table 4). Statistically significant differences were found in the yield-scaled emissions for treatment F compared to treatment U at p = 0.007. The yield-scaled GEI emissions expressed as total kg CO₂e per kg of coffee, was lower under fertilization strategy F compared to the conventional urea treatment, although not significantly different from strategy Y.

The yields obtained from treatments F, U, and Y in Sector A were comparable to the national average projected for the 2024–2025 period (5, 591 kg ha^-¹) (ICAFE, 2024). In contrast, yields in Sector B were about 50% lower (Figure 1B). Either among treatments or between field sectors, the N fertilization strategies evaluated in this study did not result in significant differences in coffee yield, similar to findings reported in Brazil, where the use of urea, urea + NBPT, and ammonium nitrate, did not affect short-term coffee productivity (Sarkis et al., 2023). The absence of yield differences does not necessarily imply that these fertilization strategies have equivalent effects on soil N dynamics, or that N gains and losses are comparable across treatments (Sadeghian-Khalajabadi et al., 2022), as additional factors (particularly soil properties) should be considered, as significant variability can occur even at the microtopographic scale within the same farm in Andisol soils (Francisco et al., 2023).

The sharp contrast between Sectors A and B could be influenced by differences in soil properties, primarily soil texture and extractable P levels (Tables 1, 3). Higher clay and silt content may affect key factors such as nutrient retention and availability (Le et al., 2023). The lower coffee yield observed in Sector B may be associated with its sandier texture (sandy loam), which likely reduced nutrient retention compared to the medium-textured Andisols considered optimal for coffee cultivation in Costa Rica (ICAFE, 2008), as observed in sector A (Supplementary Table 3). Despite the relevance of soil texture for coffee productivity, diagnostic analyses of this parameter are not commonly practiced by local producers. Moreover, the proportion of Andisols with sandy loam texture in Costa Rica is not reported in the national datasets published by the Costa Rican Coffee Institute (ICAFE). Also, typical characteristics of Andisols were observed, such as low bulk density and a high phosphorus (P) fixation capacity, the last one evidenced by extractable P levels below a critical threshold of 100 mg kg^-¹ in both sectors (Cabalceta and Cordero, 1994; Hifnalisa et al., 2020). Andisols are known to contain minerals such as allophane, imogolite, and halloysite, which contribute to P fixation. The higher total iron (Fe) content observed in Sector A may be attributed to a greater abundance of alumino-ferrous minerals, which could also explain the lower initial levels of extractable Fe in Sector B, likely related to its lower clay content (Table 1) (Knabner and Amelung, 2014; Hifnalisa et al., 2020).

Phosphorus plays a key role in supporting the growth and yield of coffee plants, as it not only contributes directly to plant development but also enhances the uptake of other essential nutrients. P deficiency has been reported to reduce crop productivity by 30–40% (Tang and Riley, 2021; Hifnalisa et al., 2020; Nurcholis et al., 2024). In this study, the final P content in Sector B confirmed its status as the most limiting nutrient affecting productivity (Table 3). Additionally, both the initial and final soil pH values in the two sectors were within the “very strongly acid” to “strongly acidic” range, conditions under which inorganic phosphorus tends to precipitate with Fe and aluminum (Al) or becomes adsorbed onto clay minerals and oxide surfaces, thereby reducing its availability to plants. This effect was more pronounced in Sector B, which exhibited a lower buffering capacity (Msimbira and Smith, 2020; Quispe et al., 2025).

Following the implementation of N fertilization strategies, available P levels increased in Sector A (Table 3), whereas no significant changes were observed in Sector B. Notably, available P in Sector A closely tripled relative to its initial value, suggesting that improvements in soil fertility may partly explain the higher coffee yields recorded in this sector. These findings are consistent with those of Le et al. (2023), who reported that increases in available soil P were positively and proportionally correlated with coffee yield. Furthermore, their study indicated that a higher proportion of fine soil fractions (clay and silt) was associated with enhanced coffee productivity, supporting the previously noted textural differences between the two sectors.

Daily CO₂ fluxes during the evaluated phenological stages are comparable to the average daily fluxes reported for coffee cultivation in Brazil, which ranged from 82.1 to 206.7 mg CO₂ m^-² h^-¹ (Vitória et al., 2019).

In this study, CO₂ emissions increased consistently after each fertilization event, reflecting a rapid stimulation of soil microbial activity in response to nutrient inputs. The rise in mineral N likely enhanced microbial N cycling, which is tightly linked to carbon metabolism, thereby intensifying soil respiration. These findings suggest that fertilization acts as a short-term driver of CO₂ fluxes, underscoring the importance of nutrient management strategies that minimize unintended carbon losses (Sosulski et al., 2020).

During the maturation stage, coffee plants undergo root senescence and a reduction in vegetative growth, while up to 70% of the annual dry matter is allocated to fruit development. At this stage, nutrient partitioning toward reproductive organs can reach levels four times greater than those directed to new vegetative tissues, with fruits immobilizing up to 95% of the absorbed K, P, and N. This shift in nutrient allocation likely reduces nutrient availability in the rhizosphere, which may partly explain the lower CO₂ fluxes observed after the third fertilization event, corresponding to the late fruit development phase (Burgos et al., 2024; Magalhães et al., 2025).

The observed negative CH₄ fluxes indicate that the coffee production system acted as a methane sink, regardless of the nitrogen fertilization strategy applied. This behavior is consistent with previous findings showing that agroforestry coffee systems often capture CH₄ emissions compared to open-field systems, mainly due to differences in microclimate, soil organic matter dynamics, and carbon availability (Berhanu et al., 2023). The low bulk density of Andisols also promotes a high proportion of air-filled pore space, creating aerobic conditions that inhibit the activity of methanogenic microorganisms (Filho et al., 2025).

The processes of organic matter formation in Andisols generate carbon forms that are not readily available or mineralizable by soil microbes. These soils typically contain large amounts of humus, high concentrations of active aluminum complexed with organic matter, and low pH, which together enhance the microbial stability of organo–metallic associations (Aran et al., 2001). The resulting slow mineralization rates may favor methanotrophic bacteria capable of oxidizing methane diffusing from the atmosphere, further contributing to the negative CH₄ fluxes observed. Consequently, methane production is likely suppressed in these soils, contrasting with the low but positive emissions often reported for other soil orders under coffee cultivation (Berhanu et al., 2023; Le Mer and Roger, 2001).

Increment of daily fluxes of N₂O was likely more pronounced after the two first fertilization due to the stimulation of soil microbial activity driven by greater N availability under favorable soil moisture and temperature conditions (Montenegro, 2020; Nurcholis et al., 2024; Filho et al., 2025). However, these responses are inherently complex and may also be influenced by additional factors, including interactions between autotrophic root respiration, heterotrophic microbial activity, and various soil properties such as texture, pH, temperature, moisture content and organic matter levels (Chinchilla-Soto et al., 2021; Quiñones-Huatangari et al., 2022).

Daily fluxes of N₂O patterns suggest that the combination of N addition, heavy rainfall, and brief drying periods plays a key role in triggering N₂O emissions. Distinct peaks in N₂O emissions were observed following the first and second fertilizer applications, coinciding with periods of heavy rainfall and subsequent brief drying phases (Figures 2C–D; Supplementary Table 4). The post-rainfall emission patterns are well documented in agricultural soils and reflect the dynamic response of nitrification and denitrification processes to changing moisture conditions (San Martin Ruiz et al., 2021), whose activity is influenced by rainfall through its effect on soil moisture. In this study, soil moisture remained within the 40–60% WFPS range (Supplementary Table 1), likely supporting the coexistence of oxic and anoxic microsites and enabling both processes to occur simultaneously (Quiñones-Huatangari et al., 2022; Sarkis et al., 2023). When heavy rains occur, waterlogged conditions drive high N₂O productions by denitrification which remains dissolved in topsoil layer until it is likely released of N₂O when soil dries (Giraldo-Sanclemente et al., 2025).

In this study, the soil water-filled pore space (WFPS) remained within the 40–60% range, which does not promote anaerobic conditions and therefore limits CH₄ production (Filho et al., 2025). The negative cumulative CH₄ fluxes observed are consistent with previous findings in Costa Rican coffee systems on Andisol soils, which also showed net soil methane sequestration that was not associated with nitrogen fertilization (Hergoualc’h et al., 2008). Although the exclusion rate of CH₄ data appears relatively high, it did not influence the overall patterns of CH₄ flux. To confirm this, comparative analyses were performed both including and excluding the omitted data. In both cases, the coffee soil system consistently functioned as a CH₄ sink, independent of soil sector or N fertilization treatment. The excluded data often showed negative fluxes with small and inconsistent slopes, suggesting that they likely resulted from ambient air fluctuations rather than actual CH₄ uptake by the soil. Considering that each treatment included six replicates, this adjustment was implemented as part of our quality-control procedure without compromising the robustness of the dataset. Additionally, the reported exclusion rate includes the 4% of CO₂ data initially removed under the first-tier gas quality criterion; thus, the CH₄-specific exclusion rate was 22%.

The N₂O emissions observed across the treatments, encompassing the three main fertilization periods of the coffee crop, fell within the range reported by Quiñones-Huatangari et al. (2022) (0.2–12.8 kg N ha^-¹ yr^-¹) in their review of fertilization-related N₂O emissions from coffee cultivation in Costa Rica, Nicaragua, and Ecuador. Cumulative N₂O fluxes did not differ significantly between treatments F and Y, mainly due to the temporal variability of daily emission peaks observed during the first four days following each of the three fertilization events (Figures 2C, 3C). These fluctuations were closely associated with rainfall events followed by short drying periods. This alternating emission pattern reflects the dynamic response of nitrification and denitrification processes to changing soil conditions, ultimately resulting in similar cumulative emissions between both treatments by the end of the cycle. In contrast, the highest cumulative N₂O emissions were observed in treatment U, where a higher nitrogen rate was applied as conventional urea during the first fertilization event (285 kg N ha^-¹ in treatment U vs. 106 kg N ha^-¹ in treatments F and Y, respectively) (Figures 2C, 3C). Finally, in the second fertilization event, when all treatments received the same N rate, the reduction of N₂O emissions from treatments Y and F compared to treatment U aligns with findings observed in coffee systems in Brazil, where the use of urea treated with NBPT and ammonium nitrate reduced N₂O emissions by 50.6% and 78.5%, respectively, relative to conventional urea (Sarkis et al., 2023).

Under the experimental conditions at the Alsacia farm in the Central Region of Costa Rica, the three fertilization strategies implemented by producers did not result in statistically significant differences in coffee yield. Nevertheless, at least biannual monitoring is necessary to generate robust evidence on the long-term effects of these strategies on crop productivity.

The findings of this study emphasize the pivotal role of soil properties in regulating coffee yield, highlighting the need for site-specific nutrient management to simultaneously support productivity and environmental sustainability. In Andisols, soil texture and P availability represent key determinants that can guide more efficient fertilization practices as enhancement of the P application could maximize the use of N (Zhang et al., 2017).

In addition, the selection of N fertilization strategies should integrate agrometeorological forecasts. In tropical regions such as Costa Rica, where the onset of the rainy season can lead to rapid and pronounced increases in N₂O emissions, the timing of N application emerges as a critical management variable.

Collectively, these results reinforce the importance of adapting nutrient management strategies to local soil and climatic conditions to foster productive, resilient, and climate-smart coffee agroecosystems.