Denitrification as the dominant process in nitrous oxide production in the water column of two eutrophic reservoirs

Reservoirs are important sites for nitrogen cycling and a significant global source of the potent greenhouse gas nitrous oxide (N₂O). They receive nitrogen inputs from agriculture and urban sources, fueling N₂O production via nitrification, denitrification, and photochemodenitrification. However, existing estimates of N₂O production in reservoirs remain uncertain because most studies have focused on N₂O in rivers or lake sediments, often overlooking the water column of lentic systems. Here, we present the first integrated assessment of N₂O production pathways in reservoir water columns using stable isotope tracer incubations alongside analyses of in situ natural abundance of nitrogen pools and functional genes involved in nitrification (amoA) and denitrification (nirS), across two eutrophic reservoirs with contrasting morphometries. We used ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ tracers to quantify rates of N₂O production, nitrification, and nitrate reduction at the beginning and the end of the stratification period. Notably, nitrate concentration decreased by up to 49 % over the two months. N₂O production from ammonium ranged from 0.02 to 48.6 nmol N L⁻¹ d⁻¹, while N₂O production from nitrate varied from 0.2 to 61.0 $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$. High rates of nitrification, nitrate reduction to nitrite, and rapid nitrite turnover were observed, with total N₂O production significantly correlated with nirS gene abundance. A strong positive correlation was found between δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ and both N₂O concentration and nirS abundance. These findings reveal that denitrification and nitrite dynamics play a central role in N₂O formation within reservoir water columns, advancing understanding of nitrogen loss and greenhouse gas emissions from lentic systems.

1 Introduction

Reservoirs created by damming rivers are an important global source of the greenhouse gas nitrous oxide (N₂O) to the atmosphere (Li et al., 2024; Wang et al., 2023). N₂O is about 273 times as potent as carbon dioxide for atmospheric warming on a 100-year time horizon (IPCC, 2021), and is the main driver of stratospheric ozone depletion (Ravishankara et al., 2009). Reservoirs receive substantial nitrogen (N) loading from agriculture and urban areas in their watersheds, processing it throughout different microbial and abiotic pathways, and then emitting back a fraction to the atmosphere as dinitrogen gas (N₂) and, significantly, as N₂O (Leon-Palmero et al., 2025; León-Palmero, 2023). Reservoirs accounted for 50 % (i.e., 0.44 Tg N yr⁻¹) of the total increase in N₂O emissions from inland waters between 1900 and 2010 (i.e., 0.89 Tg N yr⁻¹) (Wang et al., 2023). This rapid rise in N₂O emissions from reservoirs is linked to the growing number of reservoirs worldwide (Lehner et al., 2011), as well as an increase in N₂O production within these systems (Wang et al., 2023). Nevertheless, current estimates of N₂O emissions remain highly uncertain because they rely on limited datasets, and direct measurements of N₂O production rates in these reservoirs are scarce. Compared to other inland waters such as lakes and rivers, reservoirs have received far less attention, despite processing a disproportionately large fraction of N (Harrison et al., 2009), leading to elevated N₂O production rates and substantial emissions (Beaulieu et al., 2015; León-Palmero et al., 2020b, 2023; Rodríguez-Velasco et al., 2024). In fact, in Mediterranean reservoirs, N₂O emissions can occasionally surpass the combined climatic forcing of CO₂ and CH₄ (e.g., Iznájar reservoir, León-Palmero et al., 2020b). A recent study even estimated that N₂O accounted for more than 80 % of the total GHG emissions from hydroelectric reservoirs in China in 2020 (Chen et al., 2025). Therefore, it is crucial to quantify these production rates and understand the factors controlling N₂O production in reservoirs, especially considering the global increase in reservoir construction (Zarfl et al., 2015).

Microbial transformations that lead to the production and consumption of N₂O include ammonia oxidation, nitrifier denitrification, and denitrification, and they are affected by the availability of N-substrates, dissolved oxygen (DO), and phosphorus availability (Beaulieu et al., 2015; Codispoti, 2010; Ji et al., 2018; León-Palmero et al., 2023). N₂O is a byproduct of ammonia oxidation to nitrite (i.e., first step of nitrification), which is performed by ammonia-oxidizing bacteria (AOB) and ammonia-oxidizing archaea (AOA) in oxygenated waters (Könneke et al., 2005; Kowalchuk and Stephen, 2001), with the latter dominating in Mediterranean reservoirs (León-Palmero et al., 2023). At low oxygen concentrations, nitrifiers increase the yield of N₂O production, relative to the ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) oxidized, by nitrifier denitrification (via AOB), hybrid formation (AOA), or hydroxylamine oxidation (AOA), although some details of the reactions remain unresolved (Stein, 2019; Wan et al., 2023; Ward, 2013). Lastly, denitrification is the reduction of nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) to nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$), nitric oxide (NO), N₂O, and N₂, coupled to organic matter oxidation. Hence, denitrification can act as a source or sink of N₂O depending on the rate of N₂O reduction to N₂, which is catalyzed by the enzyme N₂O reductase. Denitrification is an anaerobic pathway, and oxygen regulates the activity of the denitrifying enzymes, especially the N₂O reductase (Bonin et al., 1989; Zumft, 1997). Therefore, at low but non-zero oxygen concentrations, N₂O reductase might be inhibited, promoting partial denitrification and resulting in net N₂O production. Moreover, many bacteria can denitrify in both oxic and anoxic conditions (Hochstein et al., 1984; Lloyd et al., 1987), and the presence of denitrifying bacteria has been demonstrated in the oxic and anoxic water column of lakes (Junier et al., 2008; Kim et al., 2011; Pajares et al., 2017) and reservoirs (León-Palmero et al., 2023).

Moreover, other specific factors may influence the production, accumulation, and emission of N₂O in reservoirs, such as morphometry (i.e., depth and shape) and water residence time (Hayes et al., 2017; Liang et al., 2019). The morphometry of a reservoir and water residence time affect thermal and oxygen stratification, as well as N₂O storage in the water column. Deep reservoirs can produce and accumulate large concentrations of N₂O in the hypolimnion during thermal stratification, particularly under anoxic conditions and high N concentrations. In contrast, denitrification can be a sink of N₂O in the anoxic hypolimnion when N concentration is low (Beaulieu et al., 2015; León-Palmero et al., 2023). Shallow systems tend to emit N₂O continuously due to weak thermal stratification and less capacity to accumulate N₂O. Further studies on N₂O production in the water column of reservoirs with different morphometries are required to improve our knowledge of N₂O emissions. To address this gap, we present the first integrated assessment of N₂O production pathways in reservoir water columns, combining stable isotope tracer incubations with analyses of in situ natural abundances of the N pools and functional genes involved in N₂O cycling to quantify N₂O production rates and trace the origin of the N₂O in the water column of two reservoirs. We used ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ to quantify the rates of N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and ammonia oxidation to nitrite and nitrate; and ¹⁵N-NO${}_{\mathrm{3}}^{-}$ to trace the formation of N₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ from NO${}_{\mathrm{3}}^{-}$ reduction. Incubations were performed at three depths at the beginning and end of summer stratification. We selected a shallow and a deep reservoir (Cubillas and Iznájar, respectively) located in watersheds with high N inputs, both of them monomictic with significant emissions and concentrations of N₂O (León-Palmero et al. 2020b, 2023), providing an ideal setting to explore N₂O cycling.

2 Material and Methods

2.1 Study reservoirs, morphometry, and watersheds

This study was conducted in southeastern Spain (Fig. S1) in two monomictic reservoirs with contrasting morphometries. Cubillas (37.27° N, 3.68° W) is a small and shallow reservoir with a surface area of 1.94 km² and a total capacity of 19 hm³ (mean depth = 9.66 m). Iznájar (37.26° N, 4.33° W) is a big and deep reservoir with a surface area of 26 km² and a total capacity of 981 hm³ (mean depth = 37.55 m) (open database IDEAndalucia; https://www.ideandalucia.es/portal/, last access: 13 March 2024). Both reservoirs are impacted by large agricultural and urban areas in their watersheds, which results in large inputs of N and phosphorus (León-Palmero et al., 2020b, 2023). More information about the watersheds, morphometry, and water column characterization is provided in previous studies (e.g., León-Palmero et al., 2020a, b).

We sampled the water column of these reservoirs at the beginning (4 and 9 July) and the end (5 and 7 September) of the summer stratification in 2018. During the study period, intense human usage caused a decline in the volume and water level in both reservoirs, although this decline was more evident in the smaller reservoir (i.e., Cubillas). Cubillas reservoir decreased in volume from 17 hm³ in July to 11 hm³ in September and experienced a 3.4 m reduction in the water level. The hydraulic residence time during the study period was 83 d. Iznájar reservoir decreased in volume from 575 hm³ in July to 480 hm³ in September, with a 5.4 m reduction in the water level. The hydraulic residence time was 255 d during this period. The reservoir volumes and water levels on specific dates were obtained from the Confederación Hidrográfica del Guadalquivir open database (CHG; https://www.chguadalquivir.es/saih/, last access: 13 March 2024).

2.2 Vertical profiles and biogeochemical characterization

We sampled the water column near the dam, in the open water of the reservoir, at the same location during both the July and September campaigns. First, we conducted a vertical profile of the water column using a Sea-Bird 19plus CTD profiler, obtaining continuous measurements of temperature (°C), dissolved oxygen (DO, µmol L⁻¹), and conductivity (µS cm⁻¹) in the reservoirs. Based on the temperature and DO profiles, we sampled three depths representing the epilimnion, oxycline, and hypolimnion or bottom waters. Water was collected at these three depths using a 5 L UWITEC bottle for further analyses and incubation experiments.

Samples for dissolved N₂O analysis were taken in 250 mL air-tight Winkler bottles in duplicate, preserved with a solution of HgCl₂ (final concentration 1 mmol L⁻¹) to inhibit biological activity, and sealed with Apiezon^® grease to prevent gas exchange. Samples were stored in the dark at a controlled temperature (25 °C) for less than six months until analysis at the University of Cádiz. Dissolved N₂O concentration was measured using headspace equilibration in a 50 mL air-tight glass syringe in triplicate in each bottle from each sample. N₂O concentration was quantified using a daily calibrated gas chromatograph (Bruker^® GC-450) as detailed in a previous study (León-Palmero et al., 2023).

Water samples for chemical and biological analysis were maintained at 4 °C until arrival at the laboratory. Particulate material from 500 to 1000 mL of water was filtered through pre-combusted (450 °C for 3 h) Whatman GF/F glass-fiber filters with a nominal pore size of 0.7 µm. Chlorophyll a (Chl a) was extracted from the filtered material and measured following the standard method (APHA, 1992). To obtain the cumulative Chl a (a proxy for fresh organic matter exported to the water column) in the whole water column (mg Chl a m⁻²) from the discrete depths, we summed the concentration of Chl a of each stratum using the trapezoidal rule (León-Palmero et al., 2020a). Dissolved organic carbon (DOC), ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ were assayed in the filtered water. Samples for DOC determination were acidified with phosphoric acid (final pH <2) and measured by high–temperature catalytic oxidation using a Shimadzu total organic carbon analyzer (Model TOC-V CSH) (Álvarez-Salgado and Miller, 1998). ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration was assayed using the UV spectrophotometric method at the wavelength of 220 nm and correcting for DOC absorbance at 275 nm (APHA, 1992). ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations were measured by Inductively Coupled Plasma Optical Emission Spectrometry at the Centro de Instrumentración Científica of the Universidad de Granada.

2.3 Natural abundance of stable isotopes (δ¹⁵N and δ¹⁸O)

Two 60 mL glass serum bottles per depth were collected after overflow without headspace and poisoned with HgCl₂ to analyze the natural isotopic composition (δ¹⁵N) of the ambient pools of N₂O, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Samples were maintained in darkness at room temperature for under six months before shipment to Princeton University for analysis. A 3 mL headspace was created with He before measuring the N₂O, including standards with a known amount of N₂O gas and internal standards for ¹⁵N-N₂O. The total N₂O in each bottle was extracted by purging with helium for 35 min at 38 mL min⁻¹. Then, N₂O was trapped by liquid nitrogen and isolated from interference by gas chromatography (Frey et al., 2020; Ji et al., 2015). We detected the nitrogen masses 44 (i.e., ⁴⁴N₂O representing ¹⁴N¹⁴N¹⁶O), 45 (i.e., ⁴⁵N₂O representing ¹⁴N¹⁵N¹⁶O or ¹⁵N¹⁴N¹⁶O), and 46 (i.e., ⁴⁶N₂O representing ¹⁵N¹⁵N¹⁶O), and the isotope ratios $\mathrm{45}/\mathrm{44}$, $\mathrm{46}/\mathrm{44}$ with a GC-IRMS system (Delta V Plus, Thermo). Standards in 20 mL glass vials with a known amount of N₂O gas were measured every two to three samples to calibrate for the N₂O concentration. The total N₂O concentration and ${}^{\mathrm{45}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}{/}^{\mathrm{44}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ and ${}^{\mathrm{46}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}{/}^{\mathrm{44}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ ratios were converted to moles of ⁴⁴N₂O, ⁴⁵N₂O and ⁴⁶N₂O. Both δ¹⁵N-N₂O (‰) vs. Air-N₂ and δ¹⁸O-N₂O (‰) vs. Vienna Standard Mean Ocean Water (VSMOW) were determined. Isotope measurements were linearity and offset corrected using an internal N₂O reference gas with known isotopic composition. The N₂O reference had the following isotopic composition: δ¹⁵N = $-\mathrm{0.65}\pm \mathrm{0.08}$ ‰ and δ¹⁸O = 37.37±0.27 ‰ present in ⁴⁵N₂O and ⁴⁶N₂O. Ideally, two known N₂O reference gases would have been used for correction; however, due to this limitation, natural abundance isotope data were used to analyze trends in the sample dataset, rather than making comparison with previous studies.

The natural isotopic composition of the ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pools (i.e., δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ and δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$) were determined by converting these compounds to N₂O and analyzing the isotopic composition of the resulting N₂O. ${\mathrm{NO}}_{\mathrm{2}}^{-}$ was converted to N₂O by using the azide method (McIlvin and Altabet, 2005). Sample size was adjusted to contain 10 nmol of ${\mathrm{NO}}_{\mathrm{2}}^{-}$, transferred into 20 mL glass vials, and purged with He for 10 min. The ${\mathrm{NO}}_{\mathrm{2}}^{-}$ was then converted to N₂O using sodium azide in acetic acid. During this reaction, one N from azide is transferred into the N₂O molecule; hence the resulting values were corrected by multiplying by 0.5. The ¹⁵N-N₂O generated was measured on a Delta V Plus (Thermo) as described above.

We used the denitrifier method to convert ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to N₂O (Granger and Sigman, 2009; Sigman et al., 2001; Weigand et al., 2016). The method is based on the isotopic analysis of the N₂O generated from the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ by denitrifying bacteria that lack N₂O-reductase activity (i.e., Pseudomonas chlororaphis). Sample size was adjusted to 20 nmol nitrate ${\mathrm{NO}}_{\mathrm{3}}^{-}$. The ¹⁵N-N₂O generated was measured on a Delta V Plus (Thermo) as described above. We included known ${\mathrm{NO}}_{\mathrm{3}}^{-}$ isotope international standards (USGS34 and IAEA N3) and converted them to N₂O using the denitrifier method to correct δ¹⁵N-N₂O values.

2.4 Functional genes

The abundance of unique functional genes involved in N₂O cycling was quantified using quantitative PCR (qPCR), similarly to a previous study (León-Palmero et al., 2023). We pre-filtered water samples through 3 µm pore size filters, and concentrated the samples by centrifugation, then extracted DNA following Boström et al. (2004), and applied PCR and qPCR to assess presence, and abundance of target genes. We used standard reaction mix recipes, thermocycling conditions, and primer requirements specified by the manufacturer. Specific primers were selected from studies performed in natural freshwater samples when available. DNA from pure cultures was used as positive controls and for qPCR standard preparation. We targeted ammonia oxidizers using the archaeal amoA gene, as AOA dominated over AOB in these reservoirs (León-Palmero et al., 2023). Comammox amoA genes were targeted in PCR assays using degenerate PCR primers for clades A and B (Pjevac et al., 2017), but no positive control was available in this case. The nirS gene abundance was used as a proxy for denitrifiers, while nosZ gene (Clade I) abundance was assessed only at the deepest layer, assayed only bacteria reducing N₂O to N₂. More details on the DNA extraction method, qPCR quantification, primers, specific conditions, standards, and positive controls are provided in the Supplement (Sect. S1).

2.5 Experimental setup of ¹⁵N tracer incubations

Reservoir water from the three depths was drawn from the sampling bottle into 60 mL glass serum bottles after overflow. Once in the lab, samples from oxic water depths (refer to Table 1) were purged uncapped for 2 min to remove excess N₂O, and a 3 mL headspace with ambient air was maintained after being exposed to ambient air for 30 min. Samples from anoxic waters were sealed with butyl rubber septa and crimped with aluminum seals immediately after filling. In these samples, a 3 mL helium headspace was retained after purging for 4 min. The serum bottles were weighed before and after filling them to account for the exact water volume in each sample. Table 1 compiles the incubation setup, conditions, and concentration of inorganic nitrogen added in each treatment. In the first treatment, we injected nine bottles from the same depth with ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ tracer (¹⁵NH₄Cl ≥ 98 atom % 15N, Sigma Aldrich) to a final concentration of 0.5 µmol L⁻¹, obtaining a fraction labeled of the substrate pools between 0.1 and 1.0. In this treatment, we also added ¹⁴N-${\mathrm{NO}}_{\mathrm{3}}^{-}$, equivalent to 0.10 of the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pool. In the second treatment, ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ tracer (K¹⁵NO₃, 98 atom % 15N, Sigma Aldrich) was injected to obtain a fraction labeled of the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pool about 0.10. We also added ¹⁴N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ to a final concentration of 0.5 µmol L⁻¹. Samples were incubated in the dark at the in situ temperatures from 13 to 26 °C (Table 1).

The first treatment (${}^{\mathrm{15}}\mathrm{N}-{\mathrm{NH}}_{\mathrm{4}}^{+}+{}^{\mathrm{14}}\mathrm{N}-{\mathrm{NO}}_{\mathrm{3}}^{-}$) was performed at all the depths (n=12), but the second treatment (¹⁵N-NO${}_{\mathrm{3}}^{-}$ + ¹⁴N-${\mathrm{NH}}_{\mathrm{4}}^{+}$) was performed only at the oxycline and hypolimnion (n=7, Table 1). Incubations were terminated by adding 0.1 mL saturated mercuric chloride (HgCl₂) to two bottles at t₀ (≈0.25 h), two at t₁ (≈1–3 h), two at t₂ (≈12 h), and three at t₃ (≈24 h). All samples were stored at room temperature in the dark for less than six months and shipped to the laboratory at Princeton University for further analysis.

Table 1 Incubation conditions and concentration of inorganic nitrogen compounds added in each treatment. Concentrations are measured in µmol N L⁻¹. More details are provided in the main text.

2.6 ¹⁵N-N₂O production rates from ¹⁵NH${}_{\mathrm{4}}^{+}$ and ¹⁵NO${}_{\mathrm{3}}^{-}$

The total N₂O in each incubation bottle was extracted by purging with helium and measured with a GC-IRMS system (Delta V Plus, Thermo) as explained above. We included standards in 20 mL glass vials with a known amount of N₂O gas every two to three samples to calibrate for the N₂O concentration. The total N₂O concentration and ${}^{\mathrm{45}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}{/}^{\mathrm{44}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ and ${}^{\mathrm{46}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}{/}^{\mathrm{44}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ ratios were converted to moles of ⁴⁴N₂O, ⁴⁵N₂O and ⁴⁶N₂O. N₂O production rates for each treatment were calculated from the slope of the increase in mass 45 and 46 during the linear phase over the four timepoints. The N₂O production ($R_{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$, $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) was calculated according to the following Eq. (1) (Santoro et al., 2020):

$$\begin{array}{}\text{(1)}& R_{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}={\left(F_{\mathrm{N}}\right)}^{-\mathrm{1}}\left({\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{45}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}{\mathrm{\Delta }t}}+\mathrm{2}{\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{46}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}{\mathrm{\Delta }t}}\times {\left(F_{\mathrm{N}}\right)}^{-\mathrm{1}}\right)\end{array}$$

where Δ⁴⁵N₂O and Δ⁴⁶N₂O represent the variation in the concentration of ⁴⁵N₂O and ⁴⁶N₂O over the incubation time (Δt), and the F_N represents the fraction of ¹⁵N in the initial substrate pool (${\mathrm{NH}}_{\mathrm{4}}^{+}$ or ${\mathrm{NO}}_{\mathrm{3}}^{-}$), which is assumed to be constant over the incubation time. The equation includes an extra factor of (F_N)⁻¹ to account for the probability of ⁴⁶N₂O production, which is proportional to (F_N)⁻². Natural abundance 1000 ppm N₂O carrier gas (50 µL in He) was injected before measurement to trap the produced labeled N₂O and to ensure a sufficient mass for isotope analysis.

2.7 ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ production

After N₂O analysis, we analyzed the samples incubated with ¹⁵NH${}_{\mathrm{4}}^{+}$ and ¹⁵NO${}_{\mathrm{3}}^{-}$ for ¹⁵NO${}_{\mathrm{2}}^{-}$ production to determine the rates of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation to NO${}_{\mathrm{2}}^{-}$ (ammonia oxidation), and NO${}_{\mathrm{3}}^{-}$ reduction to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (first step of denitrification). The ${\mathrm{NO}}_{\mathrm{2}}^{-}$ was converted to N₂O by using the azide method (McIlvin and Altabet, 2005), and the ¹⁵N-N₂O generated was measured on a Delta V Plus (Thermo) following the procedure and corrections described earlier. The rates of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation to NO${}_{\mathrm{2}}^{-}$ ($R_{{\mathrm{NO}}_{\mathrm{2}}^{-}\mathrm{from}{\mathrm{NH}}_{\mathrm{4}}^{+}}$, $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) and first step in denitrification ($R_{{\mathrm{NO}}_{\mathrm{2}}\mathrm{from}{\mathrm{NO}}_{\mathrm{3}}^{-}}^{-}$, $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) were calculated following Eqs. (2), (3):

$$\begin{array}{}\text{(2)}& {\displaystyle }{R}_{{\mathrm{NO}}_{\mathrm{2}}^{-}\mathrm{from}{\mathrm{NH}}_{\mathrm{4}}^{+}}={\left(F_{{\mathrm{NH}}_{\mathrm{4}}^{+}}\right)}^{-\mathrm{1}}{\displaystyle \frac{\mathrm{\Delta }\left[{}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{\mathrm{\Delta }t}}\text{(3)}& {\displaystyle }{R}_{{\mathrm{NO}}_{\mathrm{2}}^{-}\mathrm{from}{\mathrm{NO}}_{\mathrm{3}}^{-}}={\left(F_{{\mathrm{NO}}_{\mathrm{3}}^{-}}\right)}^{-\mathrm{1}}{\displaystyle \frac{\mathrm{\Delta }\left[{}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{\mathrm{\Delta }t}}\end{array}$$

where Δ[¹⁵NO${}_{\mathrm{2}}^{-}$] represents the variation in the concentration of ¹⁵NO${}_{\mathrm{2}}^{-}$, $F_{{\mathrm{NH}}_{\mathrm{4}}^{+}}$ represents the fraction of ¹⁵NH${}_{\mathrm{4}}^{+}$ in the initial substrate pool, $F_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$ represents the fraction of ¹⁵NO${}_{\mathrm{3}}^{-}$ in the initial substrate pool, and Δt is the incubation time. Each rate was calculated from the first two time points, and two or three replicates per time point. Additionally, we also calculated the turnover time of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (${\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}^{-}}$, days), which represents the average time required to replace the nitrite pool given the measured production rate following Eq. (4):

$$\begin{array}{}\text{(4)}& {\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}^{-}}={\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{R_{{\mathrm{NO}}_{\mathrm{2}}^{-}\mathrm{from}{\mathrm{NO}}_{\mathrm{3}}^{-}}}}\end{array}$$

where [${\mathrm{NO}}_{\mathrm{2}}^{-}$] represents the concentration of NO${}_{\mathrm{2}}^{-}$ (nmol N L⁻¹), and $R_{{\mathrm{NO}}_{\mathrm{2}}^{-}\mathrm{from}{\mathrm{NO}}_{\mathrm{3}}^{-}}$ represents the production rates of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (nmol N L⁻¹ d⁻¹).

2.8 ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ production

¹⁵NO${}_{\mathrm{3}}^{-}$ production rate was measured by the increase in ¹⁵NO${}_{\mathrm{3}}^{-}$ in the samples incubated with ¹⁵NH${}_{\mathrm{4}}^{+}$. We converted ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to N₂O using the denitrifier method (Granger and Sigman, 2009; Sigman et al., 2001; Weigand et al., 2016), and analyzed the resulting ¹⁵N-N₂O following the previously outlined procedure and corrections. Net production of ¹⁵NO${}_{\mathrm{3}}^{-}$ ($R_{{\mathrm{NO}}_{\mathrm{3}}^{-}\mathrm{from}{\mathrm{NH}}_{\mathrm{4}}^{+}}$, $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) is referred to here as nitrification (i.e., it includes the two-step process of oxidizing ammonium to nitrite to nitrate) and was calculated following Eq. (5):

$$\begin{array}{}\text{(5)}& R_{{\mathrm{NO}}_{\mathrm{3}}^{-}\mathrm{from}{\mathrm{NH}}_{\mathrm{4}}^{+}}={\left(F_{{\mathrm{NH}}_{\mathrm{4}}^{+}}\right)}^{-\mathrm{1}}{\displaystyle \frac{\mathrm{\Delta }\left[{}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}^{-}\right]}{\mathrm{\Delta }t}}\end{array}$$

where Δ[¹⁵NO${}_{\mathrm{3}}^{-}$] represents the variation in the concentration of ¹⁵NO${}_{\mathrm{3}}^{-}$, $F_{{\mathrm{NH}}_{\mathrm{4}}^{+}}$ represents the fraction of ¹⁵NH${}_{\mathrm{4}}^{+}$ in the initial substrate pool, and Δt is the incubation time. Each rate was calculated from the first two time points, and two or three replicates per time point.

2.9 Determination of N₂O yields

The N₂O yield during ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (N₂O−yield_Amox, %) was defined as the percent of the total N transformed to N₂O during the incubation with ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ (Eq. 6):

$$\begin{array}{}\text{(6)}& {\mathrm{N}}_{\mathrm{2}}\mathrm{O}-{\mathrm{yield}}_{\mathrm{Amox}}={\displaystyle \frac{R_{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}{R_{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}+R_{{\mathrm{NO}}_{\mathrm{2}}^{-}}}}\times \mathrm{100}\end{array}$$

The N₂O yield during nitrification (i.e., ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation to ${\mathrm{NO}}_{\mathrm{3}}^{-}$) (N₂O−yield_Nit, %) was defined as the percent of the total ${\mathrm{NH}}_{\mathrm{4}}^{+}$ transformed to N₂O during the incubation with ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ (Eq. 7):

$$\begin{array}{}\text{(7)}& {\mathrm{N}}_{\mathrm{2}}\mathrm{O}-{\mathrm{yield}}_{\mathrm{Nit}}={\displaystyle \frac{R_{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}{R_{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}+R_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}}\times \mathrm{100}\end{array}$$

The N₂O yield during denitrification (N₂O−yield_Denit, %) was calculated as follows (Eq. 8):

$$\begin{array}{}\text{(8)}& {\mathrm{N}}_{\mathrm{2}}\mathrm{O}-{\mathrm{yield}}_{\mathrm{Denit}}={\displaystyle \frac{R_{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}{R_{{\mathrm{NO}}_{\mathrm{2}}^{-}}+R_{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}}\times \mathrm{100}\end{array}$$

2.10 Data analysis

Statistical analyses were conducted in R (R Core Team, 2014) version 4.4.0. Data visualization was also performed in R, with final figure adjustments made using Inkscape (Inkscape Project, 2017). We assessed normality using the Shapiro-Wilk test of normality analysis and homogeneity of variances across groups using Levene's test. For normally distributed data with equal variances, we applied one-way ANOVA (F). When normality was met but variances were unequal, we used Welch's t test. For data that violated normality assumptions, we employed the Kruskal–Wallis rank-sum test (K–W) or the Wilcoxon test (W). Outliers were identified using the Grubbs test (G). Statistical significance was set at p<0.05. Linear regressions were used throughout the study to evaluate the rates and drivers of N₂O concentration and production. Model assumptions were assessed, and the model performance evaluated using adjusted R² values and predictor significance was determined using p values (α=0.05). Each sample was assigned a unique identifier (#1–12), which is shown in Table 1 and in the figures to facilitate data interpretation and highlight observed trends.

3 Results

3.1 Dissolved N₂O and other biogeochemical variables in the vertical profiles

The water column of Cubillas reservoir was thermally stratified in July (16.5–25.9 °C), such that DO varied dramatically with depth, with a DO peak at the top of the thermocline (400 µmol L⁻¹, 5.6 m) and decreasing concentrations until anoxia at 8 m (Fig. 1a). Dissolved N₂O concentration increased from 0.11 in the epilimnion to 6.38 µmol N L⁻¹ at the bottom of the reservoir. The decrease in the water level during the summer months due to human management presumably caused the mixing of the water column at the end of the summer, as evidenced in the homogenization of the temperature and DO profiles (Fig. 1a). Dissolved N₂O distribution remained mostly homogeneous in September, ranging from 0.22 to 0.42 µmol N L⁻¹ (Fig. 1a, Table S1 in the Supplement). The water column was always supersaturated in N₂O. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration decreased significantly from July to September (Fig. 1a, Table S1). The average ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration was reduced by half, from 321.2 µmol N L⁻¹ in July to 162.4 µmol N L⁻¹ in September. ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration varied from 13.8 to 33.0 µmol N L⁻¹ (mean = 22.0 µmol N L⁻¹). ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration was below detection level at some depths, peaking at 4.3 and 6.9 µmol N L⁻¹ in bottom waters. DOC concentrations varied from 217.6 to 247.7 µmol C L⁻¹ (Table S1), and Chl a concentrations ranged from 5.4 to 18.1 µg L⁻¹ (Fig. 1, Table S1).

Iznájar reservoir's water level decreased by over 5 m in summer, but thermal and oxygen stratification persisted due to its greater depth relative to Cubillas (Fig. 1b). The water column was always supersaturated in N₂O (Table S1). Dissolved N₂O increased with depth and over time, ranging from 0.05 to 0.26 µmol N L⁻¹ in July, up to 3.60 µmol N L⁻¹ in September, with the larger increase in the hypolimnion (Fig. 1b, Table S1). ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration also decreased from July to September, from 373.7 to 329.3 µmol N L⁻¹ (average values, Fig. 1b), with the lowest values at the oxycline, where ${\mathrm{NO}}_{\mathrm{2}}^{-}$ peaked. ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was only detected in the oxycline in July and in the hypolimnion in September, with values of 5.7 and 8.7 µmol N L⁻¹, respectively. The DOC concentrations varied from 186.0 to 228.0 µmol C L⁻¹, and the Chl a concentrations from 3.8 to 12.4 µg L⁻¹ (Fig. 1, Table S1).

Figure 1 Physico-chemical profiles of Cubillas (a) and Iznájar (b) reservoirs. The color scheme for all data is the same for both reservoirs: July (orange) and September (purple). N₂O concentration (µmol N L⁻¹, mean ± standard error) and natural abundance (δ¹⁵N-N₂O, ‰), water temperature (°C), DO concentration (µmol L⁻¹), Chl a concentration (µg L⁻¹), and the concentrations (µmol N L⁻¹) and natural abundances (δ¹⁵N, ‰) of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$. The dashed lines represent the suboxic zone (DO <10 µmol L⁻¹).

3.2 Changes in concentration and isotopic composition of N₂O and inorganic nitrogen

Figure 1 and Table S2 illustrate depth distributions of DIN concentrations and isotopic compositions. Relationships between DIN concentrations and isotopic compositions are shown in Fig. 2. The natural abundance δ¹⁵N-N₂O in the Cubillas reservoir ranged from −2.1 ‰ in the bottom waters in July to 3.6 ‰ in the epilimnion in September, while the δ¹⁵N-N₂O in the Iznájar reservoir ranged from −8.7 ‰ in the hypolimnion in July to −2.3 ‰ in the hypolimnion in September (Figs. 1 and 2). The δ¹⁸O-N₂O ranged from 41.6 ‰ in the bottom waters of the Cubillas reservoir in July to 64.4 ‰ in the bottom waters of the Cubillas reservoir in September (Fig. 2b, c). δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ was consistently positive (i.e., ¹⁵N enriched pool) in all the samples analyzed, and it varied from 8.9 ‰ to 13.4 ‰ (Fig. 2e). In the Iznájar reservoir, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration also decreased from July to September, along with an increase in δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ (e.g., Fig. 2e, #7–9). In the study reservoirs, δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ varied more than δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$. In general, δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ increased with depth, showing changes in a few meters, from ¹⁵N-depleted to ¹⁵N-enriched values, except for the Iznájar reservoir in the July sampling (Fig. 1b).

Figure 2 Relationships between the concentrations of the dissolved N₂O, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (µmol N L⁻¹), and their natural isotopic compositions. Note the logarithmic scales in the N₂O concentration axis. The lines represent the trends over depth or time mentioned in the Discussion. The slope ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}:{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}=\mathrm{2.5}$ in (c) is indicative of active N₂O reduction (Ostrom et al., 2007). Correspondence between numbers and samples is shown in Table 1 and Figs. 3 and 4. In panel (a), the red line represents the trend associated with N₂O consumption, whereas the black lines represent trends associated with N₂O production.

3.3 Distribution of N₂O production and nitrification rates from ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$

N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ ranged from 0.06 to 48.57 $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in the Cubillas reservoir (Fig. 3), and from 0.02 to 3.72 $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in the Iznájar reservoir (Fig. 4) (n=12, Table S3). Ammonia oxidation rates (i.e., ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production from NH${}_{\mathrm{4}}^{+}$, $R_{{\mathrm{NO}}_{\mathrm{2}}^{-}\mathrm{from}{\mathrm{NH}}_{\mathrm{4}}^{+}}$) were only significant in Iznájar's hypolimnion in September, reaching 215.8±38.0 $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ (${\mathrm{N}}_{\mathrm{2}}\mathrm{O}-{\mathrm{yield}}_{\mathrm{Amox}}=\mathrm{0.041}$ %) (Table S3). In contrast, significant nitrification rates (i.e., ${\mathrm{NO}}_{\mathrm{3}}^{-}$ production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$, $R_{{\mathrm{NO}}_{\mathrm{3}}^{-}\mathrm{from}{\mathrm{NH}}_{\mathrm{4}}^{+}}$) were detected at all study depths except in the hypolimnion of Iznájar in September (Figs. 3 and 4, Table S3). Nitrification rates varied from 6.1 to 56.1 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in Cubillas, and from 0.0 to 36.7 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in the Iznájar reservoir. The nitrification rates were significantly higher in July (mean ± SD = 24.6±19.4 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) than in September (7.3±6.7 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$), and in Cubillas (mean ± SD = 22.2±17.9 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$), than in the Iznájar reservoir (9.6±13.6 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) (p<0.05, in both cases). The N₂O yields during nitrification (N₂O−yield_Nit) varied from 0.000 % to 0.086 %, with the maximum yield observed in the bottom waters of Cubillas in July (Table S3). The production of N₂O from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was significantly correlated with the in situ ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration except in the hypolimnion of both reservoirs in September (n=10, adj R²=0.44, p<0.05) (Fig. 5a). These two samples, which were excluded from this analysis, contained the highest ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations (>6 µmol L⁻¹). The N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was an exponential function of the nitrification rates (Fig. 5b, adj R²=0.60, p<0.01).

3.4 Distribution of N₂O production and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction rates from ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$

N₂O production from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ varied from 0.2 to 18.1 $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in the Cubillas reservoir, and from 0.4 to 61.0 $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in the Iznájar reservoir (Figs. 3 and 4, Table S3). The highest rates were detected in the oxyclines. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (i.e., first step of denitrification, $R_{{\mathrm{NO}}_{\mathrm{2}}^{-}\mathrm{from}{\mathrm{NO}}_{\mathrm{3}}^{-}}$) varied from 13.7 to 33.2 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in Cubillas, and from 10.1 to 28.6 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in the Iznájar reservoir. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction rates were significantly higher in July (27.5±7.0 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) than in September (12.2±1.9 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) (p<0.05). This decrease in the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction rates was accompanied by a decrease in the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration from July to September in both reservoirs. Among all the samples, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ turnover varied from 0.2 d in the hypolimnion to 4.1 d in the oxycline of Iznájar in September (Table S3). The N₂O yield of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction (N₂O−yield_Denit) varied from 0.001 % to 0.132 % in the Cubillas reservoir, and from 0.003 % to 0.603 % in the Iznájar reservoir. The maximum yields occurred in the oxycline of Iznájar reservoir in September and the oxycline-bottom waters of Cubillas in September. N₂O production from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was not significantly correlated to the in situ ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration (p>0.05).

3.5 In situ abundance of functional genes

The in situ abundance of the functional genes (archaeal amoA, nirS and nosZ) varied with depth, time, reservoirs, and with the N transformation rates (Figs. 3 and 4, Table S4). Archaeal amoA abundance ranged from 0 to 2.7×10³ copies mL⁻¹ (n=12). In the Cubillas reservoir in July, the archaeal amoA gene was detected only in the oxycline, where ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration was maximal and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ minimal. We detected the archaeal amoA gene at all three depths in September, and its abundance decreased with depth. In the Iznájar reservoir, the archaeal amoA gene was detected at all depths, with the minimum abundance in the oxycline in July. Archaeal amoA abundance wasn't correlated with the N₂O concentration (p>0.05), the N₂O production rates from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (p>0.05), or the nitrification rates (p>0.05).

The nirS abundance ranged from 4.5×10⁴ to 5.3×10⁵ copies mL⁻¹ in Cubillas, and from 8.1×10⁴ to 4.7×10⁶ copies mL⁻¹ in Iznájar (n=12). nirS was present in all the samples, and its abundance increased with depth and over time in Iznájar. The nosZ gene was only quantified in the deepest layers (n=4), where it ranged from 800 to 2.1×10³ copies mL⁻¹ and was higher in September than in July in both reservoirs. N₂O production from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was not significantly correlated with the in situ nirS gene abundance (p>0.05).

Figure 3 Vertical profiles of the N₂O concentration, production rates, marker genes (colored bars), and other relevant biogeochemical variables in the Cubillas reservoir in July (a) and September (b). Dissolved N₂O (µmol N L⁻¹, mean ± standard error), and DO concentration (µmol L⁻¹); Chl a concentration (µg L⁻¹), and DIN concentration (µmol N L⁻¹); N₂O production (nmol N L⁻¹ d⁻¹) and nitrification (${\mathrm{NO}}_{\mathrm{3}}^{-}$ production, $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) from ${\mathrm{NH}}_{\mathrm{4}}^{+}$; N₂O production (nmol N L⁻¹ d⁻¹) and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production ($\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) from ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and the abundance of the target genes (×10³ copies mL⁻¹, mean ± standard deviation). Numbers next to N₂O concentrations refer to the sample ID in Table 1. The light gray area represents the suboxic zone (DO <10 µmol L⁻¹) and the dark grey the sediment. < LOD means below level of detection. Note the logarithmic scales for some panels. nosZ gene abundance was only determined in the deepest layers. N₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production were only determined in the oxycline and hypolimnion.

Figure 4 Vertical profiles of the N₂O concentration, production rates, marker genes (colored bars), and other relevant biogeochemical variables in the Iznájar reservoir in July (a) and September (b). Dissolved N₂O (µmol N L⁻¹, mean ± standard error), and DO concentration (µmol L⁻¹); Chl a concentration (µg L⁻¹), and DIN concentration (µmol N L⁻¹); N₂O production (nmol N L⁻¹ d⁻¹) and nitrification (${\mathrm{NO}}_{\mathrm{3}}^{-}$ production, $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) from ${\mathrm{NH}}_{\mathrm{4}}^{+}$; N₂O production (nmol N L⁻¹ d⁻¹) and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production ($\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) from ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and the abundance of the target genes (×10³ copies mL⁻¹, mean ± standard deviation). Numbers next to N₂O concentrations refer to the sample ID in Table 1. The light gray area represents the suboxic zone (DO <10 µmol L⁻¹). Note the logarithmic scales for some panels. nosZ gene abundance was only determined in the deepest layers. N₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production were only determined in the oxycline and hypolimnion.

Figure 5 Drivers of N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$. (a) Exponential relationship between in situ ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration (µmol N L⁻¹) and N₂O production rates (nmol N L⁻¹ d⁻¹), (b) relationship between nitrification rates (nmol N L⁻¹ d⁻¹) and N₂O production. ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations >6 µmol N L⁻¹ are shown in open symbols but excluded from the analysis in (a).

3.6 Relationships between N₂O concentration, production, and biogeochemical markers

In both reservoirs, the higher N₂O concentrations were found in the deepest layers under suboxic conditions (i.e., DO <10 µmol L⁻¹) (León-Palmero et al., 2023; Pinti, 2014), and coincided with the highest cumulative Chl a concentration (mg Chl a m⁻²), and the highest abundances of nirS gene (Figs. 1, 3 and 4). N₂O concentration decreased exponentially as DO concentration increased (Fig. 6a), but it increased in a power function correlated with cumulative Chl a concentration (Fig. 6b). N₂O concentration was also a power function of the nirS abundance (Fig. 6c). It is thus consistent that nirS abundance showed a negative correlation with DO concentration (Fig. 6d) and a positive correlation with cumulative Chl a concentration (Fig. 6e). Total production of N₂O, calculated as the sum of the production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and from ${\mathrm{NO}}_{\mathrm{3}}^{-}$, was significantly positively correlated with the nirS gene abundance (Fig. 6f, n=11). Sample #12 was excluded of this analysis.

Additionally, there was a positive correlation between δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ and the δ¹⁵N-N₂O (Fig. 6g). We also detected a strong correlation between δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O concentration (Fig. 6h). The abundance of the archaeal amoA gene was not correlated to δ¹⁵N-NO${}_{\mathrm{2}}^{-}$ (p>0.05). In contrast, δ¹⁵N-NO${}_{\mathrm{2}}^{-}$ was significantly correlated with the nirS abundance (Fig. 6i, n=12, adj R²=0.28, p<0.05). Particularly, the nirS gene abundance explained up to 94 % of the variance in δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the Iznájar reservoir (Fig. 6i, n=6, adj R²=0.94, p<0.001).

Figure 6 Drivers of dissolved N₂O concentration and production. Dissolved N₂O concentration (µmol N L⁻¹) as a function of (a) DO (µmol L⁻¹); (b) cumulative Chl a concentration (mg Chl a m⁻²), and (c) nirS gene abundance (copies mL⁻¹). nirS abundance as function of the (d) DO, and (e) cumulative Chl a concentration. (f) Total production of N₂O (nmol N L⁻¹ d⁻¹) is a function of the nirS abundance. Note that sample #12 (Hypolimnion of Iznájar in September) in (f) is an outlier, and it was not included in the analysis. (g) δ¹⁵N-N₂O as function of the δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ (‰), (h) dissolved N₂O as function of the δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ (‰), and (i) δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ as function of nirS gene abundance. A second dashed trend line and equation have been drawn in (i) only for the Iznájar samples (n=6). Note the logarithmic scales in the x and y axes. Correspondence between numbers and samples is shown in Table 1 and Figs. 3 and 4.

4 Discussion

N loading from the surrounding watershed significantly impacts the studied reservoirs, resulting in ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations exceeding 300 µmol N L⁻¹. The water columns of reservoirs have the capacity to process and remove significant amounts of N, as shown here through changes in DIN and N₂O concentrations (Fig. 1), detection of N removal processes in ¹⁵N isotope tracer experiments, presence of functional genes encoding the loss pathways (Figs. 3 and 4), and interpretation of patterns in natural abundance of N and O isotopes in the DIN and N₂O pools (Figs. 2, 6). ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration decreased by 49 % and 12 % in Cubillas and Iznájar, respectively, in just two months, which represents a substantial net N loss. This net loss in the water column likely reflects a combination of processes, including denitrification, algal assimilation followed by sedimentation of organic matter, and other biogeochemical transformations. N removal processes also drive the production of the potent greenhouse gas N₂O. The studied reservoirs had large accumulations of N₂O in their deep waters, up to 6.38 µmol N L⁻¹ in Cubillas reservoir in July, and up to 3.60 µmol N L⁻¹ in Iznájar reservoir in September. During the study period, this accumulation of N₂O in the water column of Cubillas and Iznájar reservoirs was affected by the water column depth and thermal stratification. Many reservoirs in the Mediterranean region are subject to significant evaporation during the summer, as well as intense human management, resulting in substantial fluctuations in water level. Although both reservoirs experienced a decrease in water depth, this change affected the water column biogeochemistry only in the Cubillas reservoir, likely due to its smaller size. Use of the Cubillas reservoir caused a water-level drawdown from July to September, which reduced the hydrostatic pressure and altered the water column stratification. Unstratified conditions exposed the high N₂O deep waters to the reservoir surface, which likely led to a massive release of N₂O both directly from the reservoir and, particularly, by degassing at the dam outflow or further downstream. The dam outflow is typically located at the oxycline-hypolimnion level, where the highest concentrations of greenhouse gases are found. Unfortunately, we were unable to quantify these N₂O fluxes, but the concentration detected in bottom waters in July (6.38 µmol N L⁻¹, depth = 9.5 m) versus September (0.42 µmol N L⁻¹, depth = 6.2 m) suggests a massive release of N₂O to the atmosphere during the summer. In contrast, the Iznájar reservoir did not lose thermal stratification from July to September and developed a steep oxygen gradient and an anoxic hypolimnion throughout the summer. N₂O concentration increased throughout the water column during the summer, with the most significant increase occurring in the hypolimnion (1400 % in the hypolimnion vs. ∼300 % increase in the epilimnion and oxycline), which implies that N₂O likely remains stored in that layer, and may be emitted during the fall mixing. These hydrological patterns imply dynamic N biogeochemistry during the summer stratification, which were detected explicitly by our suite of biogeochemical measurements.

4.1 Active N₂O production indicated by ¹⁵N tracer incubations and functional genes

We detected significant production of N₂O from both ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$. The rates of N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ reported in this study are larger than those found in Lake Lugano (Frame et al., 2017) and closer to those detected in the Chesapeake Bay (Tang et al., 2022). These rates are also larger than the rates found in the eastern tropical South Pacific oxygen minimum zone (Frey et al., 2020; Ji et al., 2015). N₂O production rates were significantly correlated with the availability of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and with nitrification rates, but not with archaeal amoA gene abundance. Despite the hypolimnion of Iznájar in September (#12) being apparently anoxic, we detected a significant production of N₂O from ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ammonia oxidation, and the presence of archaeal amoA genes. This combination of processes and gene detection suggests that trace amounts of oxygen may have been present at levels below the detection limit of our oxygen sensor. Similarly, the presence of trace levels of oxygen may explain the production of N₂O from ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and the nitrification rates in the anoxic waters of Cubillas, although in that case we did not detect the presence of archaeal amoA genes. The highest amoA abundance was measured in the oxycline of Cubillas in July (i.e., 2.7×10³ copies mL⁻¹), but amoA was not detected in the surface and bottom waters within the same profile, precisely where the highest N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ occurred. The absence of detectable archaeal amoA genes in samples with high N₂O production may reflect primer bias rather than true absence of ammonia-oxidizing archaea. Previous work in San Francisco Bay revealed that dominant AOA clades were not amplified by commonly used primers, including those employed in this study (Rasmussen and Francis, 2022). It is therefore possible that important AOA lineages present in these reservoirs were missed, leading to an underestimation of amoA abundance. We did not measure the bacterial amoA gene abundance, because AOA had previously been identified as the dominant ammonia-oxidizers in the study reservoirs (León-Palmero et al., 2023). Therefore, we cannot assess the potential contribution of AOB. We tested for Comammox using specific primers and did not detect them in any sample. Additionally, sample water was pre-filtered before DNA extraction (pore size = 3 µm), which may have excluded microbes attached to particles or suspended sediment, potentially including AOA or Comammox groups.

Significant nitrification rates were detected in 11 out of 12 samples, with values similar to those found in another eutrophic freshwater system, Lake Mendota (Hall, 1986), and several orders of magnitude higher than reported open ocean nitrification rates (e.g., 0.4–10 $\mathrm{nmol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) (Small et al., 2013, and references therein). The detection of high nitrification rates, but no significant ammonia oxidation, might suggest that comammox is occurring at these depths. However, our PCR analysis showed no evidence of the presence of comammox bacteria (Fig. S2), although, because no positive control was available, we cannot completely exclude their presence. Therefore, we consider the possibility that complete ammonia oxidation could contribute to the observed nitrification rates. Alternatively, we hypothesize that the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production by ammonia oxidation was tightly coupled to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ consumption by ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidizers, such that it could not be detected in the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ pool. ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production from ammonia oxidation was only detected in one sample in which we did not detect a significant nitrification rate (i.e., hypolimnion of Iznájar reservoir in September, #12), suggesting that ${\mathrm{NO}}_{\mathrm{2}}^{-}$ could accumulate due to a decoupling of ammonia oxidation and nitrite oxidation in this sample. Ammonia oxidation is the rate-limiting step for nitrification in most systems, which is why ${\mathrm{NO}}_{\mathrm{2}}^{-}$ rarely accumulates in the environment and could explain our observed mismatch between ammonia oxidation rates and total nitrification rates (Kowalchuk and Stephen, 2001). The rates of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ production detected here were often sufficient to account for a complete turnover of the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ pool during the incubation, consistent with the idea that ${\mathrm{NO}}_{\mathrm{2}}^{-}$ did not accumulate, even though the in situ concentrations were substantial.

The production of N₂O from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was generally higher than from ammonium, suggesting that ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is the main substrate for N₂O production. The highest rates occurred in oxycline samples, where ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration was often lowest, and the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration peaked. However, the N₂O production from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was not significantly correlated with the in situ concentration of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, probably because N₂O production rates are not limited by ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability. These rates were higher than the rates found in ocean waters (Ji et al., 2015), and in the Chesapeake Bay (Tang et al., 2022), but similar to those found in the eastern tropical South Pacific oxygen minimum zone (Frey et al., 2020). Similarly, these previous studies in oxygen minimum zones found the highest rates of N₂O production close to the oxic-anoxic interface (Frey et al., 2020; Ji et al., 2015).

Denitrification is the main microbial process leading to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ removal in aquatic systems. Denitrifying bacteria (as represented by the nirS gene) were consistently found throughout the reservoir water columns and reached their highest abundances in the suboxic waters. Their abundance was not significantly correlated with the N₂O production from ${\mathrm{NO}}_{\mathrm{3}}^{-}$, likely because of the small sample size (n=7). Frey et al. (2020) found that the nirS gene was not significantly correlated with N₂O production from ${\mathrm{NO}}_{\mathrm{3}}^{-}$, but was correlated with ${\mathrm{NO}}_{\mathrm{2}}^{-}$. The total N₂O production, calculated as the sum of the production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Table S3), was significantly correlated with nirS gene abundance (Fig. 6f), highlighting the importance of denitrification in the overall production of N₂O. This is consistent with the higher production obtained from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ than from ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and with the evidence from natural abundance isotopes, discussed below. The rates of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in this study were up to 1000 times higher than those in the ocean (Füssel et al., 2012; Ji et al., 2015) and in the Chesapeake Bay (Tang et al., 2022). These eutrophic reservoirs exhibit high productivity, with elevated concentrations of NO${}_{\mathrm{3}}^{-}$ and organic matter fueling intense denitrification and N₂O production. This rapid processing activity may reflect a system-level response to external nutrient loading, whereby a portion of the nitrogen input is redirected toward atmospheric release (León-Palmero, 2023).

4.2 Natural abundance stable isotopes support the role of denitrification

In general, N₂O production by denitrification, nitrifier denitrification and bacterial nitrification produces a significant isotopic fractionation of ¹⁵N, meaning that the lighter ¹⁴N is preferentially used in N₂O production, resulting in a N₂O pool depleted in ¹⁵N relative to the respective substrate and a higher δ¹⁵N value in the substrate left behind (Wenk et al., 2013 and references therein). In contrast, AOA produce N₂O that is enriched in ¹⁵N relative to the substrate, increasing δ¹⁵N-N₂O, with an isotopic fractionation value of $\approx -\mathrm{6}$ ‰ (Santoro et al., 2011; Stieglmeier et al., 2014). At the same time, the consumption of N₂O by denitrifiers increases the proportion of ¹⁵N and ¹⁸O in the remaining N₂O pool, increasing δ¹⁵N-N₂O and δ¹⁸O-N₂O values (Wenk et al., 2016).

To identify trends over depth or time, and interpret them in relation to the processes that leave their signatures in the isotopes, each sample is identified on the cross plots with a unique number (Table 1 and Figs. 2, 3, 4, 6). The trends that we observed in the natural isotopic composition of the N species suggested that denitrification was a significant process in the water column, in agreement with the rate data. In general, the increase in the N₂O concentration with depth was coupled to the δ¹⁵N-N₂O decrease (e.g., #1–3, #5–6 or #7–9 in Fig. 1 and black trend lines in Fig. 2a), which indicates net production of N₂O by water column denitrification, nitrifier denitrification and/or bacterial nitrification. In contrast, the opposite trend occurred in Iznájar in September (#10–12, Fig. 1b and red trend line in Fig. 2a), which suggests that N₂O may be a mix of consumption by denitrifiers and production by AOA in the hypolimnion at the end of the summer. There was also an increase in the δ¹⁸O-N₂O with depth in each profile, accompanied by an increase in N₂O concentration, which suggests a parallel production and consumption of N₂O at the deeper layers. That trend was not observed in Cubillas reservoir in July, but rather a noticeable increase in the δ¹⁸O-N₂O in bottom waters from July to September along with N₂O concentration decrease (Fig. 2b, #3 and #6), indicating active N₂O reduction. Besides, many samples are located along the slope ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}:{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}=\mathrm{2.5}$ in Fig. 2c, which is indicative of active N₂O reduction (Ostrom et al., 2007). We detected the nosZ gene, which encodes the reduction of N₂O during denitrification, in hypolimnetic waters with higher abundances in September. N₂O consumption can occur in the anoxic hypolimnion of Mediterranean reservoirs and result in undersaturations up to 27 % in those with low N availability (León-Palmero et al., 2023). However, in the investigated reservoirs, the N₂O reduction by nosZ-carrying denitrifiers did not cause an undersaturation of N₂O in the investigated time frame, which is consistent with previous findings in eutrophic reservoirs with high N availability (León-Palmero et al., 2023).

In the Iznájar reservoir, the decrease in ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration coincided with the increase in δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$, suggesting that denitrification is consuming the lighter ${\mathrm{NO}}_{\mathrm{3}}^{-}$ during these months (Fig. 2e, #7–9). We detected that δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ was correlated with δ¹⁵N-N₂O (Fig. 6g), which is indicative of denitrification. Over time, as more N₂O is produced from ${\mathrm{NO}}_{\mathrm{3}}^{-}$, the NO${}_{\mathrm{3}}^{-}$ pool may get substantially enriched in ¹⁵N, and δ¹⁵N-N₂O values may also increase, creating a trend line where higher δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ corresponds to higher δ¹⁵N-N₂O values. In general, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction enriches ¹⁵N in the remaining ${\mathrm{NO}}_{\mathrm{2}}^{-}$ pool, while the production of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ may decrease its δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$. In the study reservoirs, the production of N₂O by denitrification may have enriched ¹⁵N in the remaining ${\mathrm{NO}}_{\mathrm{2}}^{-}$ pool, as evidenced by the tight coupling between N₂O concentration and δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ (Fig. 6h) and the increase in the δ¹⁵N-${\mathrm{NO}}_{\mathrm{2}}^{-}$ was correlated with the abundance of denitrifying bacteria in the reservoirs (Fig. 6i). The gene used as a marker for denitrifying bacteria (i.e., nirS) encodes the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reductase that catalyses the reduction of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ during denitrification. Thus, it acts directly on the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ pool. Furthermore, the abundance of the nirS gene in the water column was correlated with the dissolved N₂O, as we also detected in a survey of twelve Mediterranean reservoirs (León-Palmero et al., 2023). These results suggest that denitrification was the main pathway of N₂O production, and it resulted in a characteristic isotopic imprint in the remaining ${\mathrm{NO}}_{\mathrm{2}}^{-}$ pool.

In addition, the cumulative Chl a concentration, which is a proxy for the vertical export of the autochthonous organic matter produced by primary producers in the whole water column, was significantly correlated with the abundance of the nirS gene and the dissolved N₂O concentration (Fig. 6b, e). This is also consistent with our previous study in twelve reservoirs (León-Palmero et al., 2023), and may indicate that denitrification is enhanced by particulate material derived from the phytoplankton community in the water column. Several studies in marine waters have described that denitrification was affected by the quantity and quality of organic matter (Babbin et al., 2014; Ward et al., 2008). Dalsgaard et al. (2012) found that the higher denitrification rates were all found at marine stations with high Chl a levels in the overlying water, suggesting a subducted and potentially decaying algal bloom. In general, this organic matter export represents a high-quality carbon source, but also sinking particles with a surface for microbial colonization, an environment where both oxic and anoxic/low oxygen microenvironments coexist, and they even increase the probability of contact between bacteria and nitrogen (Liu et al., 2013; Xia et al., 2017).

4.3 Implications for N₂O concentration and fluxes

The highest total N₂O production in Cubillas coincided with the highest N₂O concentration at the deepest depth in both months (Fig. 3). In the deeper reservoir, Iznájar, the highest production was measured at the oxycline, where there is a strong potential for N₂O fluxes, while the highest N₂O concentrations were detected in the hypolimnion (Fig. 4). In both reservoirs, the N₂O turnover time at the oxycline was the lowest in the profile. In Iznájar, the N₂O turnover time at the oxycline was as low as 13 d in July and 8 d in September (Table S3), suggesting that the N₂O produced at this location does not accumulate there. Instead, an important fraction of the N₂O produced at the top of the oxycline may be consumed or diffuse to the top layer. This diffusive flux, together with the N₂O produced in situ in the epilimnion by microbial activity and photochemidenitrification (Leon-Palmero et al., 2025), determines the large N₂O fluxes found previously in this reservoir, reaching up to 3.6 mg N N₂O m⁻² d⁻¹, and even exceeding the CO₂ equivalent warming potential from CO₂ and CH₄ emissions combined (León-Palmero et al., 2020b).

Table 2 Total DIN loss, and N₂O produced from July to September in Cubillas and Iznájar reservoirs. Details on the calculations and assumptions are provided in the Supplement.

4.4 Scaling up to the reservoir level: how much nitrogen did the reservoirs lose?

An important feature observed in the water column of these reservoirs over the summer was the substantial decrease in the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration, suggesting an active N filter for the high N loadings. Microbial activity in the water column and the sediments of reservoirs can reduce the excess of N through emissions of N₂, primarily produced during denitrification and anammox. In this study, N₂O emissions also constitute an important loss of fixed N. Total DIN loss calculations from July to September showed that Cubillas lost 468 kg N d⁻¹, while Iznájar lost 5337 kg N d⁻¹, representing a 45 % and 11 % decrease, respectively (Table 2). The DIN loss rates (2.4 and 0.7 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) were similar or even higher than those calculated in other lakes or in the Baltic Sea (Seitzinger, 1988). Normalized to reservoir surface area, the N loss was slightly higher in Cubillas. N₂O production was two orders of magnitude higher in Iznájar than in Cubillas in terms of kg N per day, but production rates were more similar when normalized to area. In the water column of Iznájar, the percentage of the N₂O production per DIN loss was higher than in Cubillas, at 1.9 % and 0.6 %, respectively. These percentages only refer to the biologically produced N₂O in the water column and may increase if the N₂O produced in the sediments, or the N₂O produced abiotically by photochemodenitrification, which was initially described in the surface waters of these reservoirs (Leon-Palmero et al., 2025), are also incorporated in the calculation. These estimates represent a major seasonal N loss event rather than annual rates. They are based on DIN concentration differences between July and September, without considering whether the reservoirs received N inputs from their watersheds during that period. Since summer is the dry period, and drawdown of the reservoirs exceeded any input via rain or runoff, N inputs from the watersheds were likely minimal during the study period. Further details on the calculations and assumptions are provided in the Supplement (Sect. S2).

Zhou et al. (2019) described a decrease of 97 % in the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration in the water column of Zhoucun reservoir during spring (2 months), and they related the N losses to aerobic denitrification occurring in the water column. Brezonik and Lee (1968) estimated that the hypolimnion of Lake Mendota lost 312 kg N d⁻¹. Beaulieu et al. (2011) found that <1 % of denitrified N was converted to N₂O in streams. Thus, these reservoirs act as important sinks for fixed N during the summer at the landscape scale, particularly within agricultural and urban watersheds, and sources of N₂O to the atmosphere. Denitrification significantly contributed to dissolved N loss and N₂O production in the water column. Although N₂O production per unit of DIN loss was less than 2 %, the absolute amount of N₂O produced in the water column and likely emitted into the atmosphere is substantial.

5 Conclusions

Our study shows that reservoir water columns actively process and remove fixed N while producing N₂O, with denitrification as the dominant pathway. This is supported by changes in DIN and N₂O concentrations, ¹⁵N isotope tracer experiments, presence of functional genes, and patterns in natural abundance of N and O isotopes in the DIN and N₂O pools. N₂O was produced from both NH${}_{\mathrm{4}}^{-}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$, with higher rates from the latter, especially in oxycline layers. Total N₂O production, and concentration were significantly correlated with nirS gene abundance. In addition, nirS abundance and N₂O concentration were correlated with the cumulative Chl a concentration, suggesting that organic matter fuels intense denitrification and N₂O production. The patterns in natural abundance isotopes further support the predominance of denitrification. δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ was positively correlated with δ¹⁵N-N₂O, and δ¹⁵N-NO${}_{\mathrm{2}}^{-}$ increased with N₂O concentration and nirS abundance. Elevated δ¹⁸O-N₂O and δ¹⁸O:δ¹⁵N slope near 2.5, along with the detection of nosZ genes suggest active N₂O consumption in several layers, such as the hypolimnion of Iznájar reservoir. Cubillas showed the highest N₂O production and concentration at depth, likely followed by surface release during summer drawdown. In Iznájar, N₂O accumulated substantially in the hypolimnion over the summer, with peak production at the oxycline, where there is a strong potential for N₂O fluxes. Both reservoirs acted as substantial N sinks for fixed N in the water column during the summer, losing 468 and 5337 kg N d⁻¹ (Table 2). Therefore, the role of reservoirs as N₂O emitters should be characterized in more detail in future studies, especially considering their the global expansion and growing importance in N₂O budgets over the past century (Li et al., 2024; Wang et al., 2023).