Distinctive Microbial Processes and Controlling Factors of Nitrous Oxide Emission in an Agricultural River Network: Perspective in Riparian Zone Type and Season

TL;DR

This study explores how different riverbank types and seasons affect nitrous oxide emissions in agricultural rivers, focusing on microbial processes.

Contribution

The study identifies riparian zone type and season as key factors influencing microbial processes and N2O emissions in agricultural river networks.

Findings

N2O emissions were significantly higher in spring compared to autumn.

Natural riparian zones emitted more N2O than artificial ones in spring.

Wind-based models overestimated N2O emissions due to inflated IPCC factors.

Abstract

The emission of nitrogen oxides (N2O) in rivers is an important source of potent greenhouse gases. However, the mechanism at the interface between rivers and riverbanks remains unclear. This study quantified N2O emissions from natural and artificial riparian zones across seasons and explored the microbial mechanisms affecting N2O production and consumption in an intensive agricultural river network in China. Significant seasonal variability in N2O emission rates was observed (p < 0.05), with mean values of 0.56 ± 0.09 mmol·m−2·h−1 in autumn and 1.13 ± 0.32 mmol·m−2·h−1 in spring. In spring, emissions from natural riparian zones (1.38 ± 0.28 mmol·m−2·h−1) were significantly higher than those from artificial riparian zones (0.89 ± 0.05 mmol·m−2·h−1). All wind-based models significantly overestimated N2O emissions (p < 0.05) due to inflated IPCC emission factors (EF5r), exceeding measured values by 1.76–3.09 times. Dissolved organic carbon and nitrite nitrogen were identified as key environmental drivers of N2O emissions. Nitrogen fixation and ammonification accounted for 82.3% of N2O production. Network analysis revealed a dominant microbial niche containing nitrifiers, sulfate-reducing bacteria, and carbohydrate-degrading taxa. Partial least squares path modeling indicated that riparian zone type altered DOC and NO2− availability, regulated nifH and ureC gene abundances, and enhanced N2O production. These findings underscore the importance of riparian-zone-specific microbial regulation of riverine N2O emissions and demonstrate the necessity of refining EF5r estimates for agricultural river networks.

Figures

N2O emission rate (a), dissolved concentration (b), saturation values (c) and emission factor EF5r−e (d) of the overlying water with the natural and artificial riparian zones in different seasons. Differences among groups were tested using one-way analysis of variance (ANOVA), with significance set at p < 0.05. Error bars represent the standard deviation (SD) of independent biological replicates (n = 2–3). In panel (d), the dashed line indicates the IPCC default EF5r value (0.26%).

Estimated N2O emission rate based on the IPCC and field-measured N2O methods. Gas transfer velocities (k, cm h−1) were calculated using seven widely accepted wind-based models (LM1986, RH2006, W1992, F2007, CC1998, N2000, and RC2001). *, and * mean significance levels at p < 0.05, p < 0.01, and p < 0.001, respectively.

(a) The relative abundances of the microbial community at the phylum level. (b) Euclidean-distances-based PCA according to season. (c) The relative abundances of nitrifying bacteria at each sampling site. (d) Euclidian-distances-based PCA according to riparian zone type and habitat.

Co-occurrence network analysis of species based on the strong (Spearman’s ρ > 0.6) and significantly positive (p < 0.05) correlation. The nodes represent phylum (a), modules (b) and N-related bacteria at the genus level (c). The color of each connection between the two nodes represents a positive or negative Spearman correlation (blue for positive correlation, red for negative correlation).

(a) Schematic diagram of the nitrogen cycle; (b–h) Relative abundances of nitrogen functional genes in the collected samples. Differences among groups were tested using one-way analysis of variance (ANOVA), with significance set at p < 0.05. Error bars represent the standard deviation (SD) of independent biological replicates (n = 2–3).

(a) Correlation between environmental factors and N-cycle-related genes and genus. (b) PLS-PM showing the direct and indirect effects of different factors on N2O emission rate. *, , and * indicate p < 0.05, p < 0.01, and p < 0.001, respectively.

Full text

1. Introduction

Nitrous oxide (N_2_O) is a major contributor to stratospheric ozone depletion. It also has a much higher global warming potential than carbon dioxide at approximately 265 times that of CO_2_ [1]. Atmospheric N_2_O concentrations have increased markedly over the past century, reaching 329 ppb in 2019, which is 22% higher than the pre-industrial level of 270 ppb [2]. Inland waters are essential sources of atmospheric N_2_O emissions [3]. It has been estimated that global riverine N_2_O emissions range from 0.68 to 0.9 Tg·y^−1^, accounting for approximately 10–17% of total anthropogenic N_2_O emissions [4]. The Intergovernmental Panel on Climate Change (IPCC) estimated that microbial conversion of agriculturally derived nitrogen to nitrous oxide was the largest source of anthropogenic nitrous oxide in the atmosphere.

The river network N_2_O emission factor (EF_5r_) was one of the three factors defined by the IPCC (2006) [5]. The concept of EF_5r_ is hypothesized based on consideration of the nitrification or denitrification ratio in the aquatic environment and subsequent production of N_2_O. After several evaluations and revisions, the EF_5r_ value was adjusted to the now commonly used value of 0.26% in 2019 [6]. When estimating EF_5r_, the IPCC assumes that aquatic N_2_O was released to the atmosphere by diffusion, with its spatial and temporal variability depended by dissolved N_2_O concentration, atmosphere equilibrium, temperature changes and gas exchange between water and air [4,7]. However, the derivation at least partially neglected the influence of various potentially decisive conditions such as the surrounding environment [3,4], water body type [8] and different forms of nitrogen [9]. Actually, significant differences between the calculated EF_5r−e_ and the IPCC default EF_5r_ have been initially identified in different environments. For example, some studies have reported that the IPCC approach may underestimate N_2_O emissions in river systems affected by human activities [4,10]. On the contrary, it has also been suggested that N_2_O production may be more likely to be overestimated [9,11,12]. Together, this information supported the critical hypothesis that EF_5r_ should incorporate more variables for adjustment when refining national GHG inventories and quantifying global N_2_O budgets based on the different surrounding environments [7].

Microbial communities are highly diverse, comprising bacteria, archaea, fungi, and protozoa. These microorganisms play pivotal roles in the microbial N cycling [13]. Microbial N transformations represent the primary biological pathways responsible for N_2_O production, including nitrification, incomplete denitrification, nitrifier denitrification, and dissimilatory nitrate reduction to ammonium (DNRA) [10,14,15,16]. Among them, incomplete denitrification is likely the globally dominant process to produce N_2_O [16]. The nitrite reductase genes nirS and nirK mediate the reduction of NO_2_^−^ to NO during denitrification, directly preceding N_2_O formation, whereas nosZ encodes nitrous oxide reductase, responsible for N_2_O reduction to N_2_ [17]. Heterotrophic denitrification is strongly regulated by carbon and nitrogen availability, with organic carbon (OC) serving as the primary electron donor and oxidized nitrogen species (NO_3_^−^ and NO_2_^−^) acting as terminal electron acceptors. Variations in OC quantity and quality, nitrogen substrates, and dissolved oxygen (DO) conditions jointly influence microbial respiration, redox dynamics, and electron allocation within denitrifying communities, thereby shaping the balance between N_2_O production and consumption [7,10,18]. For example, denitrification is advantageous when high nitrogen and low OC concentrations are associated with low oxygen [16]. In the case of insufficient organic carbon supply, nitrous oxide reductase (Nos) faces intense electron competition from the upstream electron pool, leading to N_2_O accumulation [19].

Riparian zones, as key interfaces between terrestrial and aquatic systems, exert strong control over nitrogen inputs, organic carbon availability, and redox gradients. Artificial riparian zones, particularly impermeable cement banks, can significantly alter the distribution of nirS- and nirK-harboring denitrifying communities, thereby potentially affecting N_2_O production [20,21]. Considering the rapid development of global agricultural intensification and increasing human activities, more studies are needed to quantify the contribution of artificial riparian zones to global river N_2_O emissions. Therefore, in this study, 16S rRNA gene amplicons sequencing and N-related functional gene analysis were conducted to assess the N_2_O emission rate and the vital microbial processes and factors controlling N_2_O emissions in the rivers with artificial and natural riparian zones in the Hangjiahu Plain, an intensive agricultural plain with a river network. This study aimed to (1) understand the influence of environmental factors on N_2_O emission from agricultural river networks; (2) assess the applicability of the IPCC methodology (using default EF_5r_ values to calculate N_2_O emission rates); and (3) explore the potential biological mechanisms of river N_2_O release in different riparian zone types by co-occurrence network and functional gene analysis. This study improves our understanding of the unique microbial processes underlying N_2_O emissions from different riparian zone types and illustrates the need to modify the EF_5r_ estimated for N_2_O in agricultural rivers.

2. Materials and Methods

2.1. Sample Collection and Physicochemical Index Measurement

The study was conducted in Jiashan county (120°44′22″ E–121°1′45″ E, 30°45′36″ N–31°1′12″ N) in the plain river network of the Yangtze Delta (Figure S1). Jiashan lies in the East Asian monsoon region and experiences four distinct seasons. The area is characterized by a dense network of agricultural rivers and intensive farmland drainage, making it highly representative of typical agricultural river systems in eastern China. Artificial riparian zones are widespread in the region and constitute a substantial proportion of riverbanks, providing an ideal setting to compare natural and artificial riparian influences on riverine processes. The river water area is 54 km^2^, accounting for about 11% of the Jiashan County area. The annual average temperature ranges from 14 to 22 °C, with a mean annual precipitation of approximately 1200 mm (Zhejiang Meteorological Bureau, http://zj.cma.gov.cn/, accessed on 16 January 2026). The studied rivers receive continuous runoff and drainage inputs from surrounding farmlands, offering favorable conditions for examining the transport and transformation of agriculturally derived nitrogen. Sampling was conducted during two representative seasons, from 6–9 October 2021 (autumn) and 22–24 April 2022 (spring). Spring and autumn were selected to represent two contrasting but hydrologically stable periods in agricultural river systems. Spring typically corresponds to enhanced microbial activity and increased nitrogen inputs following agricultural fertilization, whereas autumn reflects more stable nutrient conditions and reduced biological activity. For each site and sampling period, water and sediment samples were grouped into six categories, with two to three independent biological replicates per category. Sediment samples collected from natural riparian zones in spring and autumn were denoted as Spr-SN and Aut-SN, respectively. Overlying water samples from natural riparian zones were designated as Spr-WN and Aut-WN, while overlying water samples from artificial riparian zones were designated as Spr-WA and Aut-WA. The artificial riparian zones were constructed with impermeable cement revetments, had been stabilized for more than 10 years, and exhibited limited lateral hydrological connectivity with surrounding soils. Consequently, sediment accumulation at the river bottom was absent. Representative landscapes of natural and artificial riparian zones are shown in Figure S2.

Temperature, pH, redox potential (ORP), DO concentration and conductivity (HQ3d, HACH, Loveland, CO, USA), and wind speed were measured in the field. The surface water samples were collected at each sampling site using a portable water sampler (0.5 m below the air-water interface). During each sampling process, all on-site parameters were measured three times, and the average value was calculated. Surface water samples of 1 L were collected at each site using a portable water sampler, and sediment samples were collected from the top 0–10 cm layer using a grab sampler. The overlying water was pre-filtered through a 200 μm screen to remove large organisms and debris. Three surface sediment samples (<10 cm in depth) were collected with a grab sampler and homogenized to eliminate potential errors resulting from the different sampling depths. Both the overlying water and the sediment were divided into two parts. For high-throughput sequencing, one part of the water sample was filtered through a 0.22 μm pore size filter (Millipore Company, Billerica, MA, USA). The filtrated membrane and sediment were stored at −80 °C, and DNA was extracted for high-throughput sequencing and real-time qPCR analysis. In the other part, water samples and sediments were stored at 4 °C. Total organic carbon (TOC), dissolved organic carbon (DOC) and total nitrogen (TN), ammonium nitrogen (NH_4_^+^-N), nitrite nitrogen (NO_2_^−^-N) and nitrate nitrogen (NO_3_^−^-N), total phosphorus (TP), total dissolved nitrogen (TDN), total dissolved phosphorus (TDP) and phosphate (PO_4_^3−^) and suspended solids (SS) was measured using standard methods [22].

2.2. Estimation of N2O Emission Rate and EF5r

Dissolved N_2_O concentrations were determined using the headspace equilibrium technique. Detailed procedures for measuring dissolved N_2_O concentrations are provided in Supporting Method S1. N_2_O emission rates were estimated using diffusion-based approaches and the IPCC methodology. Specifically, emission rates were calculated using the thin boundary layer model based on the concentration gradient between surface water and the atmosphere, together with gas transfer velocities (k) derived from seven widely accepted wind-based models. The detailed calculation procedures and formulations for gas transfer velocity are described in Supporting Method S1 and summarized in Table S1. The equilibrium concentration of N_2_O in water (C_eq_) was estimated using the global mean atmospheric N_2_O concentration of 329 ppb. For the IPCC-based approach, the dissolved N_2_O concentration in water (C_w_) was calculated by multiplying the field-measured NO_3_^−^-N concentration by the default IPCC emission factor EF_5r_ (0.26%) [23]. This method assumes equilibrium with the global average atmospheric N_2_O concentration and applies a standardized gas exchange between water and air formulation under steady-state conditions.

The dissolved N_2_O saturation was calculated as follows:

$$[eqn]$$

where C_w_ is the measured N_2_O dissolved concentration, and C_eq_ is the corresponding concentration of N_2_O that is in equilibrium with the N_2_O in the ambient atmosphere.

The N_2_O emission rate was calculated as follows:

$$[eqn]$$

where k (m·h^−1^) represents the gas transfer velocity. Seven widely used wind-based models were applied to estimate k, following equations proposed in previous studies (Table S1).

In addition, EF_5r_ calculated (EF_5r−e_) from agricultural rivers in this study was calculated according to established IPCC methodology using the following equation:

$$[eqn]$$

where N_2_O-N and NO_3_^−^-N are the dissolved concentrations in overlying water (mg-N L^−1^).

2.3. DNA Extraction and qPCR

According to the manufacturer’s instructions, environmental DNA was extracted from the sample using an E.Z.N.A.^®^ soil DNA kit (Omega Bio-Tek, Norcross, GA, USA). The DNA concentration was measured using a Nanodrop 2000 spectrophotometer (Thermo Fisher, Waltham, MA, USA), with A260/A280 ratios between 1.8 and 2.0 and DNA yields ranging from 10 to 150 ng·g^−1^. Amplification was conducted in a 100 nL reaction system on a Wafergen SmartChip Real-time PCR system (Wafergen, Fremont, CA, USA). Standard curves were generated using serial dilutions of plasmid DNA, yielding amplification efficiencies between 90% to 105% and correlation coefficients (R^2^) greater than 0.99. The N metabolism gene abundances were quantified by real-time qPCR using a CFX96 Real-Time PCR Detection System (BioRad Laboratories, Hercules, CA, USA). Detailed information regarding the primer sets and PCR conditions is shown in Table S2. The specific details are shown in Supporting Method S2.

2.4. Statistical Analysis

The relative abundances of nitrogen-cycling functional genes were calculated by normalizing their copy numbers to the total 16S rRNA gene copy numbers in each sample, thereby accounting for differences in microbial biomass. The statistical analyses, including principal component analysis (PCA), Student’s t-test, Wilcoxon’s rank sum test, Kruskal–Wallis’ tests and α-diversity, were performed using the following R packages: stats, vegan and MASS. The co-occurrence network was explored by calculating all pairwise Spearman’s rank coefficients (ρ) among 0.97-OTUs with a relative abundance greater than 0.01% across all samples. All significant and robust correlations were visualized and explored in networks using the Gephi 0.9.2 software. The stepwise multiple regression model and the variance inflation factor were conducted using the SPSS 24.0 software (IBM Corp., Chicago, IL, USA). Partial least squares path modeling (PLS-PM) was established through the smartPLS (version 3.0) software to analyze the influence of environmental factors and functional genes on N_2_O emissions in the agricultural river network.

3. Results

3.1. N2O Emission Rate and Physicochemical Properties of the Overlying Water

N_2_O-related values of the overlying water with the different riparian zones are summarized in Figure 1. There were noticeable seasonal differences in the measured N_2_O emission rate (p < 0.05), averaging 0.56 ± 0.09 mmol·m^−2^·h^−1^ in autumn and 1.13 ± 0.32 mmol·m^−2^·h^−1^ in spring. Moreover, there were significant differences between different riparian zones in spring, averaging 1.38 ± 0.28 mmol·m^−2^·h^−1^ in the natural riparian and 0.89 ± 0.05 mmol·m^−2^·h^−1^ in the artificial riparian zones. Correspondingly, the N_2_O dissolved concentration and saturation values in the natural riparian zones (99.24 ± 13.14 μmol·L^−1^, 465.15 ± 57.13%) were slightly higher (p > 0.05) than those in the artificial riparian zones (91.88 ± 16.03 μmol·L^−1^, 430.96 ± 65.79%). The N_2_O dissolved concentrations were higher than the saturation level in all survey sites, ranging from 339.08 to 575.83% (448.05 ± 65.16%). N_2_O supersaturation was commonly observed in the agricultural river network of the Yangtze River Delta, suggesting a tendency for net N_2_O efflux to the atmosphere under the observed conditions.

N2O emission rate (a), dissolved concentration (b), saturation values (c) and emission factor EF5r−e (d) of the overlying water with the natural and artificial riparian zones in different seasons. Differences among groups were tested using one-way analysis of variance (ANOVA), with significance set at p < 0.05. Error bars represent the standard deviation (SD) of independent biological replicates (n = 2–3). In panel (d), the dashed line indicates the IPCC default EF5r value (0.26%).

The physicochemical characteristics of the overlying water were significant factors affecting the production and release of N_2_O (Figure S3). According to Pearson’s correlation analysis, the N_2_O emission rate was significantly positively correlated with wind speed (r = 0.83, p < 0.001), pH (r = 0.66, p < 0.05), and TP (r = 0.68, p < 0.05). In contrast, the N_2_O emission rate was significantly negatively correlated with DOC (r = −0.84, p < 0.001), TOC (r = −0.76, p < 0.01) and the C/N ratio (r = −0.77, p < 0.01). Furthermore, the N_2_O dissolved concentration (r = −0.83, p < 0.001) and saturation values (r = −0.64, p < 0.05) were strongly negatively correlated with SS content. Multiple stepwise regression models incorporating environmental factors were established to represent the N_2_O emission rate (Table 1). The results demonstrated that DOC, NO_2_^−^ and DO could predict the N_2_O emission rate jointly (R^2^ = 0.932, p < 0.001). Correspondingly, the most critical factor impacting the N_2_O emission rate was DOC (t1 = −9.528), which had a stronger correlation with the N_2_O emission rate than NO_2_ (t2 = 4.832) and DO (t3 = 2.929). DOC, NO_2_^−^, and DO differed significantly between natural and artificial riparian zones (Figure S4). The concentrations of DOC and NO_2_^−^ in artificial riparian zones (28.21 ± 5.27 mg·L^−1^ and 0.16 ± 0.03 mg-N·L^−1^) were significantly higher than those in natural riparian zones (18.65 ± 6.94 mg·L^−1^ and 0.05 ± 0.01 mg-N·L^−1^). However, DO in the artificial riparian zones (6.81 ± 1.18 mg·L^−1^) was significantly lower than that in the natural riparian zones (9.40 ± 1.19 mg·L^−1^). In addition, NO_3_^−^ accounts for 33.4–44.9% of TN in the artificial riparian zones, while NO_3_^−^ accounts for 60.6–80.3% of TN in the natural riparian zones.

3.2. Riverine N2O Emission Factor and N2O Emission Rates

The field-measured N_2_O emission rate was significantly lower than that measured by the corresponding IPCC methods (1.76~3.10-fold, p <0.001) (Figure 2). The field-measured results were 2.63 ± 0.43-fold lower than IPCC results in natural riparian zones and were 2.11 ± 0.33-fold lower in artificial riparian zones. In this study, F2007 was the latest and most widely accepted model for predicting the N_2_O emission rate. Among the other six N_2_O emission rate prediction models tested, only W1992 (0.62 ± 0.13 μmol·m^−2^·h^−1^) showed no significant difference with F2007 (0.75 ± 0.36 μmol·m^−2^·h^−1^) (p > 0.05). LM1986 (0.28 ± 0.05 μmol·m^−2^·h^−1^) and RH2006 (0.46 ± 0.09 μmol·m^−2^·h^−1^) were significantly lower than F2007. However, N2000 (1.01 ± 0.18 μmol·m^−2^·h^−1^) and CC1998 (3.41 ± 0.57 μmol·m^−2^·h^−1^) and RC2001 (4.17 ± 0.69 m^−2^·h^−1^) were significantly higher than F2007.

*Estimated N2O emission rate based on the IPCC and field-measured N2O methods. Gas transfer velocities (k, cm h−1) were calculated using seven widely accepted wind-based models (LM1986, RH2006, W1992, F2007, CC1998, N2000, and RC2001). , ** and *** mean significance levels at p < 0.05, p < 0.01, and p < 0.001, respectively.

EF_5r−e_, e.g., the ratio of dissolved N_2_O/NO_3_^−^, ranged from 0.101 to 0.161% (0.130 ± 0.019) and about half of the most recent value from the IPCC (0.26%; 2019) (Figure 1d). There was no significant difference in EF_5r−e_ in autumn (p > 0.05), but there was a significant difference between natural and artificial riparian zones in spring. The average EF_5r−e_ was 0.105 ± 0.007% in the natural riparian zones and 0.154 ± 0.007% in the artificial riparian zones. Furthermore, according to the stepwise multiple regression model, the 48.9% variation in the EF_5r−e_ value in each sampling site can be mainly explained by the variation in NO_3_^−^/TN (Table S3).

3.3. Microbial Community Diversity, Structure and Interaction Relationship

High-throughput sequencing was used to analyze the microbial diversity and community composition in the river network of the Yangtze River Delta Plain to unravel the biological factors controlling N_2_O emissions. There was a strong relationship between α-diversity and habitat (Table 2). Specifically, α-diversity (Richness, Evenness, Diversity) in the sediment was significantly higher than that in the overlying water (p < 0.05). In the overlying water, the Simpson and Shannon diversity values in autumn (0.98–0.98, 5.05–5.89) were significantly higher than in spring (0.96–0.97, 4.02–4.44). Notably, the planktonic community diversity of natural riparian zones was more stable than that of artificial riparian zones. Especially in spring, the diversity of microbial communities was only one-third of that in the autumn along the natural riparian zones.

The relative abundances of the microbial community at the phylum level are shown in Figure 3a. The top five phyla accounted for 58.6–67.0% of the total population in the sediment: Proteobacteria, Actinobacteria, Cyanobacteria, Bacteroidota and Chloroflexi. However, the top five phyla covered 85.0–96.8% of the total population in the overlying water. The phyla Proteobacteria, Cyanobacteria and Bacteroidota had obvious seasonal changes in the overlying water. The relative abundances of Proteobacteria, Cyanobacteria and Bacteroidota were 33.2 ± 3.4%, 28.0 ± 1.5% and 8.6 ± 0.1% in spring and 47.1 ± 7.6%, 4.6 ± 1.0% and 16.1 ± 0.9% in autumn, respectively. The abundances of Actinobacteriota and Chloroflexi were noticeably affected by the riparian zone type in the overlying water (p < 0.01). The relative abundances of Actinobacteriota and Chloroflexi were 18.2 ± 2.5% and 0.38 ± 0.2% in the natural riparian zones, while they were 23.6 ± 8.4% and 2.1 ± 1.1% in the artificial riparian zones. Regarding the dissimilarities in microbial communities, PCoA showed that the first two axes covered 59.6% of the explanation (Figure 3b,d). The succession range of planktonic bacteria was wider than that of benthic bacteria. Moreover, the most crucial influence factor of the microbial community composition was the season, the riparian zone type and the habitat type.

The high-throughput sequencing results indicated that nitrifying bacteria included AOB (ammonia-oxidizing bacteria) and NOB (nitrite-oxidizing bacteria) in this study (Figure 3c). The nitrifying bacteria were found to belong to the phyla Proteobacteria, Nitrospirota and Verrucomicrobiota. The nitrifying bacteria levels in the sediment (2.2 ± 0.6%) were much higher than those in the overlying water (0.25 ± 0.15%). In the overlying water, the amount of nitrifying bacteria in autumn (0.42 ± 0.33%) was much higher than that in spring (0.04 ± 0.01%). The abundance of nitrifying bacteria in artificial riparian zones (0.27 ± 0.17%) was much higher than that in natural riparian zones (0.03 ± 0.01%). Among nitrifying bacteria, the relative abundance of Ellin6067, Nitrospira and Candidatus_Nitrotoga were the top three in the sediment, which were 1.3 ± 0.1%, 0.4 ± 0.3% and 0.1 ± 0.9%, respectively. In the overlying water, Ellin6067 was also the most abundant (0.06 ± 0.03%), following oc32 (0.04 ± 0.03%) and Nitrospira (0.02 ± 0.02%). AOB accounted for 73.2 ± 10.5%, 70.4 ± 8.7% and 86.3 ± 13.3% of the nitrifying bacteria in SN, WN and WA, respectively, including the genus Ellin6067, mle1-7, MND1, oc32, Nitrosomonas, GOUTA6 and Nitrosospira. Ellin6067 was the predominant nitrifying bacteria in SN, WN and WA (57.9 ± 6.2%, 45.5 ± 25.2% and 26.9 ± 31.3%). NOB accounted for 25.9 ± 11.1%, 29.6 ± 8.7% and 13.3 ± 13.4% of the nitrifying, including the genus Nitrospira, Nitrotoga, and CI75cm.2.12.

(a) The relative abundances of the microbial community at the phylum level. (b) Euclidean-distances-based PCA according to season. (c) The relative abundances of nitrifying bacteria at each sampling site. (d) Euclidian-distances-based PCA according to riparian zone type and habitat.

The co-occurrence (i.e., positive and negative) associations between species (i.e., 0.97-OTUs) were explored between different riparian zone types and seasons. The average clustering coefficient (0.732) and modularity (1.155) in the co-occurrence network were higher than those in the Erdos–Rényi random networks (0.0100.021, 0.3060.460), which indicated that the microbial community network in the agricultural river network has a modular structure and “small-world” properties, e.g., high interconnectivity. The resulting co-occurrence network consisted of 154 nodes (i.e., 0.97-OTUs) and 3976 edges (i.e., correlation). The positive network consisted of 132 nodes (i.e., 0.97-OTUs) and 2721 edges (i.e., correlation), while the negative network included 22 nodes and 1254 edges, which indicated that the interaction relationship was dominated by symbiosis in the agricultural river network. Overall, the networks of OTUs belonged to Proteobacteria (39.61%), Actinobacteriota (22.08%), Bacteroidota (10.39%) and Cyanobacteria (9.09%) (Figure 4a). Network partitioning divided the network into three modules, which can be considered as three niches. Twenty of the negative species were found in Module I. A negative species was included in both Module II and Module III. There were many riparian zones with different species in Module II. For example, Limnohabitans were significantly enriched species in natural riparian zones, followed by Rhodoluna, and Rhodobacter (p < 0.05). They were 3.43-, 4.92-, and 7.59-fold that of the former in artificial riparian zones, respectively (Figure S5). Mycobacterium were the most abundant species of artificial riparian zones (p < 0.05), followed by Alsobacter and Methylocystis, which were 15.8-, 14.0-, and 5.3-fold that of the former, respectively. In Module III, 90% of Cyanobacteria belonged to the phylum. The abundance of Cyanobacteria in spring (28.0 ± 1.5%) was significantly higher than in autumn (4.6 ± 1.0%).

Co-occurrence network analysis of species based on the strong (Spearman’s ρ > 0.6) and significantly positive (p < 0.05) correlation. The nodes represent phylum (a), modules (b) and N-related bacteria at the genus level (c). The color of each connection between the two nodes represents a positive or negative Spearman correlation (blue for positive correlation, red for negative correlation).

3.4. Nitrogen Metabolism Functional Genes

To assess nitrogen metabolism, marker genes encoding key enzymes involved in nitrogen transformation processes were analyzed (Figure 5a). The organic nitrogen mineralized marker gene (gdhA) was significantly higher in spring sediments (4.12 ± 0.89%) than in autumn sediments (1.78 ± 0.66%), but there was no significant difference in overlying water (p > 0.05). ureU, a marker gene for ammonification, was significantly higher in spring sediments (6.24 ± 1.21%) than in autumn sediments (1.34 ± 0.51%). A nitrogen fixation marker (nifH) was significantly affected by seasons but had no significant relationship with riparian zone type. The nifH relative abundance was around 0.13 ± 0.05% in spring, while the relative abundance was around 0.09 ± 0.03% in autumn. In addition, the abundance of comammox amoB in artificial riparian zones (0.25 ± 0.15%) in spring was significantly higher than in other riparian zone types (0.12 ± 0.09%). The abundance of nirS in sediment was 1–2 orders of magnitude higher than that in water (7.9-fold), with no significant difference in water (p > 0.05). The abundance of nirS was 49.3–392.9-fold that of nirK in sediments and 31.52~38.00-fold that in overlying water, indicating that nirS denitrification was dominant in both natural and artificial riparian zones. The nosZ gene, responsible for the conversion of nitrous oxide to nitrogen, was significantly higher in spring sediments (2.16 ± 0.68%) than in autumn sediments (0.35 ± 0.11%), which was not significantly different in overlying water. The ratio of the N_2_O-producing gene (nirK, nirS) to the N_2_O-consuming gene (nosZ) was utilized to indicate the N_2_O-producing potential in denitrification. In this study, although there was no significant difference between nirK + nirS/nosZ in different riparian zones, different types of denitrification have significantly different effects on the N_2_O emission rate. Although nirS denitrification played an essential role in N_2_O production potential, the actual net N_2_O production importance of nirK/nosZ should not be ignored (r = 0.67, p < 0.05) (Figure S2b). Overall, according to Pearson’s correlation and multiple stepwise regression, the N_2_O emission rate was most correlated with ureC (r = 0.73, p < 0.05) and nifH (r = 0.79, p < 0.01) in overlying water, which together explained 82.3% of the variation in the N_2_O emission rate (Table 3).

(a) Schematic diagram of the nitrogen cycle; (b–h) Relative abundances of nitrogen functional genes in the collected samples. Differences among groups were tested using one-way analysis of variance (ANOVA), with significance set at p < 0.05. Error bars represent the standard deviation (SD) of independent biological replicates (n = 2–3).

Mantel tests were used to assess correlations between environmental variables, nitrogen cycling functional genes, and the genus matrix (Figure 6a). The results showed that genes were significantly correlated with the N_2_O emission rate, TN, TOC, DOC, NO_3_^−^-N, NO_2_^−^-N, C/N and NO_3_^−^/TN. In contrast, the bacterial genus was significantly correlated with ORP, pH, temperature (T), TDP, NO_3_^−^-N, SS and NO_3_^−^/TIN. Most significantly, the N_2_O emission rate was only significantly correlated with the functional genes. Among the genes, ureC was negatively correlated with TN, TOC, DOC, NO_2_^−^-N and C/N and positively correlated with NO_3_^−^/TN. The gene nifH was negatively correlated with pH, temperature (T) and SS and positively correlated with DOC. The PLS-PM model was constructed to reveal the influence of riparian zone type on N_2_O release, likely facilitating the changes in physicochemical properties and microbial community caused by the DOC concentration and then affecting N_2_O release (Figure 6b). ureC and nifH together account for the 82.3% variance in the N_2_O emission rate. The PLS-PM results showed that DOC contributed 80.3% of the variance in ureC, and DOC and NO_2_^−^ contributed 70.2% of the variance in nifH. The statistical results highlighted the influence of ammonification and nitrogen fixation on N_2_O based on different riparian zone types.

*(a) Correlation between environmental factors and N-cycle-related genes and genus. (b) PLS-PM showing the direct and indirect effects of different factors on N2O emission rate. *, *, and *** indicate p < 0.05, p < 0.01, and p < 0.001, respectively.

4. Discussion

4.1. Evaluation of IPCC-Based N2O Emission Estimates for the Yangtze River Delta River Network

Increased anthropogenic N_2_O emissions have become a global issue. Both the N_2_O dissolved concentration and N_2_O emission rate (48.78 ± 7.45 nmol·L^−1^, 0.75 ± 0.31 μmol·m^−2^·h^−1^) in this study were lower than the estimated global mean values for rivers and streams (51.9 μmol·L^−1^, 1.72 μmol·m^−2^·h^−1^) [24] but comparable to other agricultural systems (Table 4). In particular, our results were consistent with observations from moderately eutrophic rivers such as the Yangtze River mainstem [25] and the Wensum River [26] but substantially lower than values reported for highly nutrient-enriched watersheds, including the Taihu [10] and Chaohu basins [6]. Rivers and watersheds exhibiting higher N_2_O emission rates and EF_5r_ values (e.g., Taihu and Chaohu) are typically characterized by intense eutrophication, elevated nitrogen inputs, and strong organic matter enrichment, which promote high N_2_O supersaturation and emission (Table 4). In contrast, river systems with lower N_2_O fluxes, such as the Yangtze River mainstem and rivers on the Tibetan Plateau, generally experience lower nutrient concentrations, higher dissolved oxygen levels, or reduced residence times, which constrain incomplete denitrification and N_2_O accumulation. Methodological differences further influence estimated N_2_O emissions. While diffusion-based models parameterized with field-measured N_2_O concentrations (e.g., W1992 and F2007) produced emission estimates consistent with observed fluxes, the IPCC approach relies on a fixed EF_5r_ value combined with a global average atmospheric N_2_O concentration. The IPCC methodology inherently neglected local hydrological, biogeochemical, and riparian zones controls, which can introduce substantial bias when applied to heterogeneous river systems. In our study, N_2_O emission rates estimated using the IPCC methodology were higher than those derived from field-based measurements (Figure 2). The consistently lower EF_5r−e_ values observed in both natural and artificial riparian zones relative to the IPCC default EF_5r_ value further indicate that riverine N_2_O emission factors were highly context-dependent (Figure 1d). Moreover, mechanistic modeling studies of entire river networks have demonstrated that achieving the magnitude of the IPCC default EF_5r_ was not kinetically feasible in most river systems [27]. These findings suggested that future refinements of EF_5r_ should explicitly incorporate geographical context and riparian zone characteristics to better capture spatial heterogeneity in N_2_O production and exchange, thereby reducing uncertainty in interwatershed N_2_O emission assessments.

Riparian zone type and season were key factors regulating N_2_O emissions in human-impacted freshwater systems. N_2_O emission rates were significantly higher in spring than in autumn, and emissions from natural riparian zones exceeded those from artificial riparian zones during spring. This seasonal pattern likely reflected the combined effects of increasing temperature, enhanced microbial activity and greater availability of reactive nitrogen in spring. Warmer conditions stimulated microbial growth and enzyme activity, accelerating nitrogen transformation processes such as nitrification and denitrification, while increased agricultural runoff during this period further enhanced substrate supply for N_2_O production. In artificial riparian zones, N_2_O production was largely confined to the overlying water, where higher oxygen concentrations favor nitrification but limit denitrification [28]. Moreover, the natural riparian zones provided more river–riparian exchange than artificial riparian zones. Impervious/artificial riparian zones significantly improved the nitrification rate but reduced the denitrification rate [21]. The enhanced river–riparian exchange and favorable denitrification conditions explained the higher fluvial N_2_O emission rates observed in natural riparian zones, particularly during spring.

4.2. Primary Process of N2O Production in the Agricultural River Network

Microbial α-diversity patterns provide important context for understanding spatial and seasonal variability in N_2_O emissions. Higher α-diversity in sediments compared to overlying water indicated greater functional redundancy among sediment-associated communities, which may stabilize nitrogen transformation processes [29]. In contrast, the lower and more seasonally variable diversity observed in the water column, particularly during spring, reflects a more disturbance-sensitive community (Table 2). The pronounced seasonal decline in diversity observed in artificial riparian zones may further reduce functional stability and favor incomplete nitrogen transformations, thereby promoting transient N_2_O accumulation.

Functional gene analysis supported denitrification as the dominant microbial pathway underlying N_2_O production. qPCR results showed that nirS-type denitrification contributed most to N_2_O formation across sites (Figure 5). This was consistent with previous observations in freshwater systems, where the content of the nirS gene in the sediment is usually higher than that of the nirK gene (nirK/nirS = 0.003–0.020) [30,31]. In this study, nirS abundance was highest in spring, particularly in sediments of natural riparian zones, indicating that nirS-type denitrifiers dominated nitrite reduction under organic-rich and more reduced conditions. In contrast, nirK exhibited relatively higher abundances in artificial riparian zones during spring, suggesting that nirK-type denitrifiers were favored under more oxic or redox-fluctuating environments.

Seasonal variation further influenced the balance between nirK- and nirS-type denitrification. The nirK/nirS ratio was significantly higher in spring (0.34–0.65) than in autumn (0.02–0.19), likely reflecting shifts in redox conditions, substrate availability, and microbial community structure associated with rising temperatures and enhanced hydrological connectivity. nirK-type denitrifiers are often associated with aerobic or facultative aerobic microorganisms and may operate more efficiently under oxic or redox-variable conditions due to higher oxidative stress tolerance [32]. Although nirS remained numerically dominant across most samples, the enrichment of nirK in spring highlighted a complementary denitrification pathway that may enhance incomplete denitrification and promote N_2_O accumulation, particularly in riparian environments subject to seasonal hydrological and biogeochemical disturbances. The nirK/nirS ratio varied from 0.05 to 0.65 in natural riparian zones and from 0.02 to 0.35 in artificial riparian zones, indicating distinct denitrification regimes between riparian types. Despite the dominance of nirS-type denitrifiers, nirK-containing taxa played a non-negligible functional role in nitrite reduction, as supported by the significant positive relationship between nirK/nosZ and N_2_O emission rates (Figure S2). Although our study focused primarily on sediment-associated microbial communities, these non-bacterial pathways may represent an additional N_2_O source in eutrophic river systems and warrant further investigation. Recent studies have demonstrated that photosynthetic microorganisms, including microalgae such as Chlamydomonas, can produce N_2_O from NO_2_^−^ via nitrate-reductase-dependent pathways, particularly under oxic or fluctuating redox conditions [33]. These non-bacterial pathways may therefore represent an additional, seasonally relevant source of N_2_O in eutrophic river systems and warrant further investigation.

4.3. Importance of N-Related Bacteria in the Co-Occurrence Patterns

The “small-world” and highly modular network structure indicated that nitrogen cycling is driven not by isolated taxa but by coordinated microbial assemblages operating within distinct ecological niches (Figure 4). Module I represented a cooperation-dominated functional regime, as it contained several well-documented denitrifying genera, including Flavobacterium, Hydrogenophaga, Arenimonas, and Rhodobacter [34,35]. The dense positive associations within this module suggested enhanced microbial cooperation, communication, and coordinated activity in riparian zone environments. Such tightly integrated networks were likely to enhance nitrogen transformation rates and promote N_2_O accumulation, particularly when electron donor availability and redox conditions favor incomplete denitrification. The coordinated regulation at the community level may amplify N_2_O production beyond what would be expected from individual taxa alone [35]. Consistent with these observations, PICRUSt1 functional predictions indicated greater capacities for flagellar assembly, cell motility, and secretion in microbial communities from natural riparian zones compared to artificial riparian zones (Figure S6).

In contrast, Module II reflected a competition-driven and functionally partitioned regime, characterized by a predominance of negative correlations and pronounced niche differentiation. This module included nitrifiers, sulfate-reducing bacteria, and fermentative taxa that indirectly influence nitrogen cycling by regulating electron flow and substrate availability. For example, seven nitrifying bacteria were identified within Module II (Figure 4c), along with members of the phyla Firmicutes and Desulfobacterota. Sulfate-reducing genera such as Geothermobacter and Desulfatiglans may cooperate syntrophically with denitrifiers by using sulfate as an electron shuttle to facilitate organic matter degradation and nitrogen removal [36,37]. However, the dominance of competitive interactions within this module suggested weaker metabolic coupling, leading to more fragmented nitrogen transformation pathways and potentially limiting sustained N_2_O production despite the presence of functionally relevant taxa [38].

Module III represents a phototrophy-driven regime dominated by Cyanobacteria and strongly shaped by seasonal dynamics. The enrichment of Cyanobacteria in spring indicates a close coupling between primary production and downstream nitrogen cycling. Through photosynthetic exudation and biomass turnover, Cyanobacteria supply labile organic carbon, which can indirectly stimulate heterotrophic denitrification in sediments and overlying water [39]. In addition, diel oxygen fluctuations associated with photosynthesis and respiration generate transient redox gradients that favor incomplete denitrification and episodic N_2_O release. Beyond these indirect effects, growing evidence suggests that phototrophic microorganisms can directly produce N_2_O from nitrite under oxic or redox-fluctuating conditions, highlighting Module III as a potential seasonal source of N_2_O.

4.4. Mechanism of N2O Production Affected by the Riparian Zone Types

The riparian zone types affected soil exchange at the water–sediment interface, which may affect microbial community composition, metabolism and consumption, thus determining N_2_O emissions. In the first step of the N_2_O emissions, nitrogen from N_2_ is fixed to NH_4_^+^ (Figure 4a). nifH, a gene involved in nitrogen fixation, played an essential role in nitrogen input. The PLS-PM model showed that the relative abundance of the nifH gene ranged from 0.09 to 4.89%, 70.2% of which can be explained by DOC and NO_2_^−^ (Figure 6b). As a prominent electron donor, organic carbon is a major growth-limiting factor for heterotrophic or mixotrophic bacteria [40,41]. Higher N_2_O concentrations in spring than in autumn can therefore be attributed to lower river discharge, higher nitrogen availability, and a greater proportion of labile, protein-like DOM, which together promote microbial activity and N_2_O production. In contrast, increased discharge and more refractory DOM inputs in autumn dilute nitrogen substrates and suppress N_2_O accumulation [42]. Furthermore, NO_2_^−^ is an essential substrate of anammox and denitrification. The importance of NO_2_^−^ on the nifH gene may demonstrate that the nitrogen loss process facilitates the nitrogen fixation process [43].

In the second step of nitrogen transformation, organic nitrogen is decomposed by microorganisms through ammonification, releasing ammonia, while assimilation converts ammonium back into organic nitrogen. The functional genes associated with these processes were ureC (ammonification) and gdhA (assimilation). The relative abundance of ureC was significantly higher than that of gdhA, exceeding it by 9.41–64.45 times, and it was strongly regulated by DOC, which explained 80.3% of its variation (p < 0.001). Together with nitrogen fixation, ammonification represents an upstream process that controls the supply of reduced nitrogen substrates (e.g., NH_4_^+^ and NO_2_^−^), thereby exerting a strong influence on downstream N_2_O production potential. Although these pathways did not directly generate N_2_O, their variability affects N_2_O emissions by regulating substrate availability and redox conditions. Consistent with this mechanism, nitrogen fixation and ammonification emerged as the most important functional pathways influencing N_2_O release in this study. The abundances of nifH and ureC together explained 82.3% of the variability in N_2_O emissions (Table 3 and Figure 6b).

The third step of nitrogen transformation is nitrification, during which ammonia is oxidized to nitrate (NH_4_^+^ → NO_2_^−^ → NO_3_^−^) (Figure 4a). Despite their measurable abundance, nitrifying bacteria were not the primary drivers of N_2_O emissions in this study. However, except for the nitrification gene amoB, the gene amoA, hao, and nxrA concentrations of all samples were present at very low abundances or were undetectable (<0.4%) in this study. This suggested that nitrification played a limited direct role in controlling N_2_O production.

Finally, denitrification (NO_3_^−^ → NO_2_^−^ → NO → N_2_O → N_2_) converts nitrates into nitrous gas and nitrogen. Key denitrification genes, including nirS, nirK, and nosZ, were consistently detected in both natural and artificial riparian zones. Among them, nirS-type denitrification was dominant, and nirS abundance was positively correlated with TOC and DOC (p < 0.05), reflecting the importance of organic carbon as an electron donor for denitrifying bacteria. The consistently high nirS abundance across sites and seasons indicated that denitrification represents a baseline and persistent N_2_O-producing capacity in the system. However, because nirS abundance showed relatively limited spatial and temporal variability, it explained less variation in N_2_O emissions than upstream functional genes such as nifH and ureC (Figure 6b).

Denitrification is more likely to promote N_2_O accumulation under DOC-limited conditions [44]. High DOC availability favors complete denitrification, resulting in N_2_ rather than N_2_O as the final product, and it can also enhance NO_3_^−^ immobilization, thereby reducing the potential for incomplete denitrification. In contrast, low DOC conditions promote incomplete denitrification and N_2_O accumulation. Consistent with this mechanism, previous syntheses have shown that DOC. NO_3_^−^ ratios strongly regulate aquatic N_2_O dynamics, with lower ratios favoring higher N_2_O accumulation [45]. Beyond bacterial denitrification, emerging evidence suggests that photosynthetic microorganisms may also contribute to N_2_O production from nitrite. Microalgae and cyanobacteria have been shown to reduce NO_3_^−^ to N_2_O via nitrate reductase-dependent pathways under oxic or fluctuating redox conditions [33]. Given the seasonal enrichment of phototrophic taxa observed in this study, particularly in spring, phototrophic N_2_O production from nitrite may represent a complementary and previously underappreciated pathway contributing to N_2_O emissions in agricultural river networks. In contrast, genes associated with anammox and dissimilatory nitrate reduction to ammonium (DNRA), such as hzsA, nirB, and nrfA, were not detected, indicating that these pathways played a negligible role in the studied system.

By integrating microbial functional genes with environmental variables, the PLS-PM analysis revealed that regulating external carbon and nitrogen inputs, especially DOC availability and nitrogen loading, may help reduce N_2_O emissions from agricultural river networks. These results suggest that regulating external carbon and nitrogen inputs, particularly DOC availability and nitrogen loading, may help reduce N_2_O emissions from agricultural river networks. Maintaining or restoring functional riparian zones could further stabilize microbial nitrogen transformations and limit conditions that favor incomplete denitrification and NO_2_^−^ accumulation. The strong variability in EF_5r_ across riparian zone types and seasons highlights the need to refine emission factors by incorporating riparian structure, nutrient status, and microbial functional characteristics. Such refinements would improve the accuracy of regional and global greenhouse gas inventories. This study also has limitations. The limited spatial sampling density within each riparian zone type may constrain the resolution of fine-scale heterogeneity in microbial processes and N_2_O emissions. The future studies with denser spatial sampling would help improve statistical power and capture local variability. In addition, N_2_O fluxes were primarily estimated using diffusion-based approaches, and combining these with direct water–air and sediment–water flux measurements would strengthen flux validation and reduce uncertainty. Although absolute N_2_O emission magnitudes may vary among regions due to differences in climate, hydrology, and land use, the mechanisms identified here provide a transferable framework that can be tested and refined across broader spatial and temporal scales.

5. Conclusions

This study provides a comprehensive assessment of microbially mediated N_2_O production in agricultural river networks and demonstrates that N_2_O emissions were strongly regulated by coupled carbon and nitrogen dynamics. Notably, when compared with seven commonly used wind-based gas transfer models, N_2_O emission rates estimated using the IPCC default EF_5r_ value (0.26%) were significantly overestimated in the studied agricultural river network. In contrast, field-based estimates yielded much lower EF_5r−e_ values, averaging 0.105 ± 0.007% in natural riparian zones and 0.154 ± 0.007% in artificial riparian zones, underscoring the limitations of applying a uniform emission factor to heterogeneous river systems. The PLS-PM analysis further revealed that DOC and NO_2_^−^ jointly explained 70.2% of the variability in nifH abundance and 80.3% of the variability in ureC abundance and together accounted for 82.3% of the variation in N_2_O emission rates. Low DOC availability was associated with enhanced incomplete denitrification and increased N_2_O accumulation, whereas higher DOC levels favored more complete denitrification to N_2_. Meanwhile, the pronounced imbalance between ammonia-oxidizing bacteria (AOB) and nitrite-oxidizing bacteria (NOB), with AOB abundances exceeding those of NOB by approximately 2.74-fold, promoted NO_2_^−^ accumulation, thereby creating a critical bottleneck in nitrogen transformation. DOC and NO_2_^−^ contents regulate the abundance of enzyme-encoding genes (ureC and nifH) of denitrifiers, which further affect river N_2_O emissions, which all depend on the riparian zone type. In addition, nitrifying bacteria (genus Ellin6067, Nitrospira and C. Nitrotoga) SRB (genus Geothermobacter and Desulfatiglans) and carbohydrate-degrading bacteria (genus Clostridium) became the dominant niche to promote the dynamic succession of the microbial community. The modification of EF_5r_ is urgently needed in the future, as well as a comprehensive study of microbial N_2_O production processes, their interactions, and, most importantly, their impact on global warming potential.