From 2012 to 2021, we conducted various fertilizer treatments on a maize field located at Gansu Agricultural University in Gansu Province, NW China (35°28′N, 104°44′E). Within the 2021 cropping season, soil samples was collected for microbial data analysis and gathered data related to maize crops. The study site is characterized by a semi-arid climate on the Loess Plateau, with 140 frost-free days annually and an elevation of 1971 meters above sea level. The region experiences an average annual precipitation of 400 mm, a mean annual temperature of 10.8°C, an evaporation rate of 1,531 mm, and an average annual radiation of 5,930 MJ m⁻². The soil in this area originates from aeolian deposits, containing at least 50% sand and classified as Calcaric Cambisol according to FAO (1990) guidelines. Additional physiochemical properties of the soil before the start of the experiment can be found in Supplementary Table S1.

The experiment utilized a randomized complete block design with four treatments and three replicates per treatment. The nitrogen fertilizer fertilizers rates were 100 (N1), 200 (N2), and 300 (N3) kg N ha⁻¹, while no nitrogen was applied for control (N0). Each treatment was replicated three times. All plots received an application of 150 kg P₂O₅ ha⁻¹ at the time of application of first nitrogen dose. Due to the high potassium (K) content (220 mg kg⁻¹) in the soil of the region, no additional K fertilizer was applied, as it was deemed sufficient for crop growth (Yue et al., 2022). Nitrogen were supplied as urea (16% N) and triple superphosphate (16% P₂O₅) respectively. The experiment comprised 16 plots, each measuring 3 m × 14.2 m in area. Nitrogen fertilizer was applied in two stages: one-third of the total nitrogen fertilizer was evenly spread on the soil surface before maize sowing, while the remaining two-thirds were applied to the soil at the six-leaf stage of the maize. Prior to sowing, nitrogen (N), and phosphorus (P) fertilizers were hand-applied via broadcasting and incorporated by shallow cultivation, followed by harrowing. However, during the six-leaf stage of maize growth, N application was carried out using a handheld injection device positioned alongside each row of maize plants due to the presence of plastic mulch.

In late April, maize seeds of the Pioneer 335 cultivar were planted at a density of 52,500 plants per hectare. Maize Pioneer 335 cultivar are disease and pest resistant. It is also estimated to mature within a time frame of 5 to 6 months with an estimated grain yield of 34,382 kg ha⁻¹. The field was prepared with wide (0.7 m) and narrow (0.4 m) ridges covered with transparent plastic film. The transparent plastic mulch (polyethylene film) was 1.01 mils thick, 4 feet wide, in rolls 2,000 feet long. It was available in a clear (transparent) color. Holes were created in the film over the furrows to facilitate the collection of precipitation. The seeds were sown within the furrows (Supplementary Figure S1). After placing the film over the soil, it was perforated using a handheld device (Lamptey et al., 2019). Manual weeding (hand) was performed to control weeds throughout the growing season. The maize grain was harvested in late September.

During the flowering stage of maize, soil samples were collected at a depth of 20 cm from the field in the 2021 cropping season using a 5 cm diameter auger. Randomly selected bulk soil samples were collected from various locations within the field, while rhizosphere soil samples were obtained by sampling the soil adhered to the crown root. A total of 12 soil samples (three replicates for each of the four treatments) were collected. To create a uniform sample, 10 soil cores were combined from each plot. The samples underwent processing to remove stones and surface debris and were then sieved using a 2 mm sieve mesh. To prevent cross-contamination, the auger and sieve mesh were cleaned with ethanol and clean tissue paper before collecting new samples. The soil samples were promptly placed on dry ice for transportation to the laboratory. Half of each sample was air-dried for subsequent chemical analysis measurements, while the other half was stored at −80°C for microbial community analysis.

To analyze the soil pH, a mixture of deionized water and soil was prepared in a ratio of 1:2.5 (mass to volume). The pH of the resulting extract was measured using a pH meter (Mettler Toledo FE20, Shanghai, China) (Amadou et al., 2020). Total nitrogen (TN) content was determined using the Kjeldahl method, while the organic carbon content (SOC) was measured using the Walkley-Black wet oxidation method (Bao, 2000). The concentration of dissolved organic nitrogen (DON) was quantified using the multi-N/C 2100 s Analyzer (Analytik Jena, Germany). Ammonium nitrogen (NH₄⁺-N) and nitrate nitrogen (NO₃⁻-N) concentrations were determined using a spectrophotometer at a wavelength of 550 nm and 204 nm, respectively, (UV-1800, Mapada instruments, Shanghai, China) after extracting the soil with 2 M KCl (Bremner, 1996). The available phosphorus (AP) concentration was determined via the molybdenum-blue method after extracting the soil with sodium bicarbonate (Olsen et al., 1954). Soil water content (SWC) was assessed by subjecting the soil to oven-drying at 105°C for 24 h (Fudjoe et al., 2022).

The soil’s potential denitrification activity (PDA) was assessed during the maize cropping season of in 2021 using a modified version of the acetylene inhibition method (Philippot et al., 2011). In this approach, approximately 4 grams each of soil (equivalent to the dry weight) from a total of 12 soil samples (four treatments with three replicates) were incubated in an incubator with a solution consisting of KNO₃ (50 μg NO₃⁻-N g⁻¹ dry soil), glucose (0.5 mg C g⁻¹ dry soil), and sodium glutamate (0.5 mg C g⁻¹ dry soil) in a 150 mL sterile flask. The mixture was gently mixed and then incubated at 28°C in an incubator. The atmosphere within each sterile flask was evacuated to establish anaerobic conditions and inhibit N₂O-reductase activity, and purged with a 90:10 He-C₂H₂ gas combination. Gas samples were collected at the start of the incubation after 2 h and analyzed for N₂O concentrations using a gas chromatograph (Agilent, 7890A, United States) equipped with an electron capture detector. The PDA value was determined as ng N-N₂O produced per hour per gram of dry soil.

The aboveground dry biomass was assessed by subjecting it to oven drying until a constant weight was achieved during the maize cropping seasons of 2020 and 2021. The grain yield was measured after harvesting and air-drying all maize cobs obtained from the plot. The nitrogen use efficiency (NUE) was calculated by subtracting the nitrogen uptake in the treatment without nitrogen fertilizer from the nitrogen uptake in the treatment with nitrogen fertilizer and then dividing this difference by the nitrogen application rate (Liu et al., 2018b).

Analyses of functional genes (nirS and nosZ) were done with clean chimera tags from the 12 samples (four treatments × three replications). To obtain the total genomic DNA from the soil samples, 0.5 g of the dry weight equivalent from a total of 12 soil samples each was processed using the Power Soil^® DNA Isolation Kit (MoBio, Carlsbad, CA, United States). After extraction, the Wizard DNA Clean-Up System (Axygen Bio, United States) was employed to purify the DNA. The purity level of the DNA was at a ratio of 1.8. The DNA samples were then stored at −80°C until analysis. The copy numbers of the nirS and nosZ genes were determined using quantitative polymerase chain reaction (qPCR) with specific primer sets listed in Supplementary Table S2. For qPCR, the 20 μL reaction mixture consisted of 7.2 μL of aseptic water, 0.4 μL of each primer (10 mM), 10 μL of GoTaq^® qPCR Master Mix (Promega, United States), and 2 μL of template DNA. Each qPCR reaction was carried out in triplicate. Following incubation of the qPCR products at 72°C for 5 min, electrophoresis was conducted on a 2 percent agarose gel with ethidium bromide staining to facilitate detection. The known copy numbers of the target gene were used to create standard curves through a 10-fold serial dilution of plasmids containing the gene. The purified amplicons were combined in equal proportions and processed using an Illumina Miseq^® PE300 platform (Illumina, San Diego, United States). The amplification efficiencies and r² values exceeded 90 and 0.99%, respectively.

DNA sequencing was utilized to examine the abundance and composition of nirS and nosZ genes. The forward primers were modified with a unique 7 bp barcode sequence, and the concentration of the PCR products was measured using a TBS-380 fluorometer. Subsequently, the PCR products were diluted and subjected to paired-end sequencing on an Illumina MiSeq sequencer (Shanghai Personal Biotechnology, Co., Ltd., Shanghai, China). Additional details on the primer pairs, reaction mixtures, and thermal cycling conditions for amplifying all six genes can be found in Supplementary Table S2. Following amplification, the PCR products of all six genes were isolated from agarose gels and purified using a universal DNA Purification Kit. Quality screening of the raw sequences was performed using Quantitative Insights Into Microbial Ecology (QIIME) to ascertain any low-quality sequences (Caporaso et al., 2010). To identify chimeric assembled sequences, the Usearch tool was employed, while the FrameBot tool from the Ribosomal Database Project (RDP) was used to screen for chimeric sequences (Edgar, 2013). The FunGene Pipeline, as described by Edgar and Flyvbjerg (2015), was utilized to exclude low-quality sequences. Operational Taxonomic Units (OTUs) were defined using the CD-HIT approach within MOTHUR, with a 3% difference threshold applied to nucleotide sequences (Schloss et al., 2009). The nirS and nosZ sequences obtained spanned 31,446–56,507 bp and 34,219–59,361 bp of the valid codes, respectively. MiSeq^® sequencing of nirS obtained chimera-free reads, and average length of valid tags ranging from 223.89–224.91 bp while nosZ had 392.63–519.26 bp. The sequences of the nirS and nosZ genes were deposited in the NCBI Sequence Read Archive (SRA) database and can be accessed using the specific accession numbers SRR18481671 and SRR18481638, respectively.

The data from the bulk and rhizosphere soil were subjected to one-way analysis of variance (ANOVA) using SPSS (version 22, IBM Corporation, Chicago, United States, 2013). Duncan’s multiple range tests (DMRT) were applied to determine significant differences (p ≤ 0.05) among the treatment means for the copy number of nirS and nosZ composition, maize grain and biomass. Alpha diversity indices, including the Shannon index, Simpson index, and Chao1 richness for nirS and nosZ genes, were determined using R software (version 3.5.3). The effects of soil physicochemical properties and nirS and nosZ Composition were evaluated using redundancy analysis (RDA) with the “vegan” package in R software. A correlation test examined the association between soil chemical properties, soil water content (SWC), maize grain and biomass. Additionally, the significant taxa in the nirS and nosZ Composition in the bulk and rhizosphere soil were identified using a co-occurrence network.

The results of four fertilization treatment samples with three replications were pooled together. The Operational Taxonomic Units (OTUs) present in each treatment replicate were chosen for network analysis. Together, Pearson correlation, Bray–Curtis, and Kullback–Leibler dissimilarities were utilized. A true co-occurrence network was defined as a statistically significant association between species when the correlation coefficient (r) exceeded 0.8 or was below −0.8, with a value of p of 0.01. Permutation and bootstrap distributions were generated with 1,000 iterations to assess the reliability of the network edges. The network was visualized using the Fruchterman–Reingold algorithm implemented in Gephi (version 0.9.2). Various topological properties of the network, including the number of nodes and edges, average clustering coefficient, average degree, average path length, closeness centrality, network centrality, and modularity, were calculated. OTUs with higher degrees and closeness centrality were identified as prospective keystone taxa. Modules are defined as the structure of networks which measures the strength of division of microbial communities (also called groups or clusters) (Berry and Widder, 2014).

Additional analysis was done using random forest modeling to ascertain the important predictors of nirS and nosZ composition and maize yield. These predictors included soil variables, and the unexpected forest package was employed (Liaw and Wiener, 2002). The significant forecasters obtained from the random forest analysis were then used to analyze the direct and indirect influences of physiochemical soil properties on biomass, network modules, nirS and nosZ communities, NUE, and maize productivity using AMOS 21.0 in SPSS (SPSS, Inc., Chicago, IL). Before modeling, the data distribution was tested for normality. A structural equation model (SEM) was employed, and the model’s fitness was assessed using the chi-square test (χ², p > 0.05), root mean square error of approximation (RMSEA), and goodness-of-fit index (GFI) (Sahoo, 2019).

The analysis of variance revealed significant differences in most soil indices (TN, NO₃⁻-N, AP, SOC, DON, and SWC) among different application rates of N fertilization in both the rhizosphere and bulk soil, except for NH₄⁺-N. Considering the alterations in pH due to physicochemical factors, no fertilizer (N0) treatment had a higher pH than the fertilizer treatments (N1, N2, and N3). Furthermore, biotic factors such as root exudates alter the rhizosphere (8.66) pH as compared to bulk soil (8.56). The soil pH significantly varied across N fertilization rates, ranging from 8.32 to 8.66 as compared to the initial pH (8.37) of the soil before experimentation (Table 1; Supplementary Table S1). The levels of SOC, NO₃⁻-N, AP, and DON tended to be higher in the rhizosphere compared to the bulk soil and increased with N2 and N3 treatments. Across the bulk and rhizosphere soil, the N3, N2, and N1 treatments led to an increase (13.7–27.5%; 33.5–36.1%), (12.6–21.4%; 22.7–28.9%) and (12.2–13.1%; 13.7–16.8%) compared to N0 in SOC and NO₃⁻-N, respectively. However, N3 and N2 significantly increased TN, while the N2 treatment had higher AP levels than the N0 treatment in bulk and rhizosphere soil. The SWC was significantly higher in the bulk soil than in the rhizosphere soil across the N fertilization treatments (Table 1).

The application of N fertilization treatments resulted in a significant increase (p < 0.05) in maize productivity compared to the no fertilizer (N0) treatment during the 2020 and 2021 cropping seasons (Supplementary Table S3). Grain yield in both seasons showed a substantial increase under N3 and N2 treatments, with (64.8, 63.2%) and (61.7, 60.8%), respectively, compared to the N0 treatment. The aboveground biomass during the same cropping seasons exhibited a similar trend, with the highest increase observed under the N3 treatment, followed by N2 and N1 treatments, with increases of (62.8, 64.9%), (58.9, 61.9%), and (48.1, 51.2%) compared to the N0 treatment, respectively (Table 2). However, there was no significant difference in grain and biomass yield between the N2 and N3 treatments (p > 0.05; Supplementary Table S3).

The nitrogen use efficiency (NUE) exhibited a similar pattern as maize productivity. The N3 and N2 treatments resulted in a fold increase in NUE of (2.16, 1.95%) and (1.38, 1.35%) compared to the N1 treatment (Supplementary Table S3). Significant differences (p < 0.05) were observed among the treatments for grain yield, biomass, and NUE across the 2 years. The variation in grain yield, biomass, and NUE primarily stemmed from the year effect and the interaction between the year and treatment (Supplementary Table S3).

The PDA index in both bulk and rhizosphere soil showed a significant increase (p < 0.05) under N3 (35.7, 42.2%), N2 (45.5, 37.1%), and N1 (24.6, 19.3%) fertilization treatments compared to the N0 treatment, respectively, in the 2021 cropping seasons (Figure 1).

The highest peaks of N₂O flux emissions were observed in July, while the lowest levels were recorded in October and September across all N fertilization rates during the 2020 and 2021 cropping seasons (Figure 2A). Moreover, the release of N₂O flux emissions was significantly higher (p < 0.05) under nitrogen fertilizer application (N1, N2, and N3) treatments compared to the no N fertilizer (N0) treatment during the growing seasons (Figure 2A). During the 2020 and 2021 cropping seasons, the highest N₂O emission flux was observed under N3 and N2 treatments, measuring (126.4 and 115.7 mg m⁻² h⁻¹) and (94.2 and 86.1 mg m⁻² h⁻¹), respectively. However, the lowest N₂O flux was recorded under N0 treatment, at 48.4 and 17.3 mg m⁻² h⁻¹, respectively (Figures 2A,B). Additionally, the cumulative N₂O emissions in the 2020 cropping season were 32.4% higher with N3 treatment compared to N0, and in the 2021 cropping season, they significantly increased under N3, N2, and N1 treatments by 33.3, 37.3, and 22.8%, respectively, relative to N0 treatment. The emission rates ranked in the following order: N3 > N2 > N1 > N0 (p < 0.05; Figures 3A,B).

The rhizosphere soils exhibited higher abundance and diversity in the nirS and nosZ composition compared to the bulk soil under N fertilization rates (Figures 4A–D, 5A–D, and Table 2). Specifically, the abundance of nirS-harboring denitrifiers significantly increased by 20.4, 27.6, and 28.6% under N1, N2, and N3 treatments, respectively, compared to N0 treatment in the bulk soil (Figure 4A). In the rhizosphere soil, N3 and N2 treatments showed significant increases in the abundance of nirS and nosZ-harboring denitrifiers compared to N0 treatment (Figures 4C,D). Moreover, the abundance of nosZ-harboring denitrifiers under N2 treatment (20.5%) was significantly higher than N3 (15.4%) and N1 (11.7%) treatments in the bulk soil under N fertilization rates (Figure 4B). According to the diversity Chao1 and Shannon indices, there were substantial differences in the nirS- and nosZ-denitrifiers across N fertilization rates, with N3 and N2 treatments predominating in both bulk and rhizosphere soil (p < 0.05; Table 2). Compared to the bulk soil, the rhizosphere soil indicated the nirS denitrifier OTU diversity was considerably higher in the N0, N1, N2, and N3 treatments (p < 0.05; Table 2). The nosZ-denitrifier OTUs in the bulk soil were 30.4% (N3), 15.6% (N2), and 24.1% (N1) more diverse than in the N0 treatment, indicating a significant difference (p < 0.05; Table 2). However, no significant distinction existed between the N fertilization treatments in the rhizosphere soil (p > 0.05; Table 2).

Across the range of N fertilization rates, the genera Cupriavidus (22.5%), Bradyrhizobium (18.0%), Rhodanobacter (12.7%), Azospira (9.4%), Herbaspirillum (7.7%) and Zoogloea (7.1%) were the most prevalent nirS-harboring denitrifier community in the bulk soil (Figure 5A). When compared to N0 and N2 treatments, the relative abundance of genus Cupriavidus and Zoogloea significantly (p < 0.05) increased under N3 and N1 treatments, whereas genus Bradyrhizobium and Rhodanobacter showed the reverse trend in the bulk soil (Figure 5A). Whereas in the rhizosphere soil, genera Luteimonas (19.3%), Rubrivivax (16.5%), Anaerolinea (12.6%), Arenimonas (10.8%), Acidovorax (7.9%), Tepidiphillus (7.3%), and Rhodanobacter (5.3%) were the major nirS-harboring denitrifiers. Compared to N0 and N3 treatments, the relative abundance of genera Luteimonas was significantly (p < 0.05) higher under N2 and N1 treatments. In the rhizosphere soil, genera Anaerolinea and Rubrivivax established dominance under the N3 treatment (Figure 5C).

The nosZ-harboring denitrifiers in the bulk soil were mainly dominated by the genera Azospirillum (23.2%), Mesorhizobium (17.2%), Burkholderia (16.0%), and Herbaspirillum (14.0%) across N fertilization rates (Figure 5B). Figure 5B illustrates that genera Azospirillum and Mesorhizobium were significantly (p < 0.05) more abundant under N3 and N2 treatments than other N fertilizer rates, while genera Burkholderia showed the opposite pattern in the bulk soil. In the rhizosphere soil, the most prevalent nosZ-harboring denitrifiers were genera Rhizobium (14.2%), Microvirga (12.9%), Bradyrhizobium (10.2%), Burkholderia (8.8%), Mesorhizobium (8.3%), and Cupriavidus (6.3%; Figure 5D). Genera Rhizobium and Burkholderia were significantly (p < 0.05) more abundant in the N2 and N0 treatments compared to the N3 and N1 treatments. On the other hand, genera Microvirga and Bradyrhizobium exhibited higher abundance in soils treated with N3 compared to other N fertilizer rates (Figure 5D).

Redundancy analysis (RDA) results demonstrated that NO₃⁻-N (19.1%), SOC (16.6%), pH (15.7%), and N₂O emissions (13.2%) significantly (p < 0.05) influenced the composition of the nirS denitrifier composition in the bulk soil (Supplementary Table S4). The genus Cupriavidus exhibited a strong positive relationship with SOC and N₂O emissions while showing a negative correlation with NO₃⁻-N and SWC. Furthermore, NO₃⁻-N was positively connected with the genera Rhodanobacter. The SWC was positively associated with the genus Bradyrhizobium, and SOC exhibited a negative relationship (Supplementary Figure S2A). The nirS-harboring denitrifier composition structure in the rhizosphere soil showed that NO₃⁻-N (17.8%), SWC emissions (17.2%), pH (12.4%), DON (11.8%), and N₂O emissions (10.1%) significantly increased (Supplementary Table S5, p < 0.05). Genera Luteimonas and Rubrivivax had a positive relationship with pH and DON and opposite relationships with SWC, NO₃⁻-N and N₂O emissions, respectively. However, the genera Arenimonas and Acidovorax demonstrated the reverse pattern (Supplementary Figure S2C).

The findings of RDA in the nosZ-harboring denitrifier composition showed that SOC (20.5%), TN (18.2%), NO₃⁻-N (15.6%), and pH (13.5%) in the bulk soil impacted significantly (p < 0.05; Supplementary Table S4). Genera Azospirillum had a positive relationship with NO₃⁻-N and an opposite relationship with pH, whereas Mesorhizobium had the reverse pattern. Genera Burkholderia was also positively associated with SOC (Supplementary Figure S2B). RDA study in the rhizosphere soil revealed that TN, soil pH, NO₃⁻-N, SOC, and SWC accounted for 10.5, 10.3, 14.2, 13.4, and 12.8% of the variation in the nosZ denitrifying community, respectively (Supplementary Table S5; p < 0.05). Genera Rhizobium and Cupriavidus were identified to have a positive relationship with TN and pH. Genera Microvirga, on the other hand, exhibited a positive relationship with SOC (Supplementary Figure S2D).

Distinct modules were formed in co-occurrence networks based on treatments to analyze the soil denitrifying community. These networks were used to investigate the relationships between modules and functional groups of taxa in the bulk and rhizosphere soil (Figures 6A–D, 7A–D). Only the OTUs with a relative abundance above 0.01% in at least one replicate were included in the network. Supplementary Tables S6, S7 provide information about the module nodes and edges based on the topological characteristics of the nirS- and nosZ-harboring community in the bulk and rhizosphere soil. The nirS- and nosZ-harboring community networks in both soil types exhibited a higher number of positive associations (415 and 936 edges; 539 and 1,160 edges) compared to negative associations (8 and 34 edges; 77 and 67 edges), respectively (Supplementary Tables S5, S6). In the bulk soil, Module I exhibited a strong interconnectedness among OTUs within the nirS-harboring network, surpassing the relationships observed in other modules (Figure 6B). In the case of the nosZ-harboring network in the bulk soil, all four modules demonstrated a high level of OTU relationships, with Modules II and IV particularly prominent (Figure 7B). Additionally, in the rhizosphere soil, the nirS- and nosZ-harboring network communities displayed high OTU relationships within Modules I and II, while Modules III, IV, and V exhibited comparatively weaker relationships (Figures 6D, 7D).

The nirS-harboring denitrifier network in the bulk soil revealed three potential keystone taxa: genera Cupriavidus, Rhodanobacter, and Bradyrhizobium, indicating their critical roles within the network (Figure 6A). Similarly, the nosZ-harboring denitrifier network in the bulk soil indicated five keystone taxa: genera Azospirillum, Mesorhizobium, Burkholderia, Ensifer, and Pseudomonas emphasizing their significance within the network (Figure 7A). Furthermore, in the rhizosphere soil, the nirS-harboring denitrifier network was primarily influenced by the genera Luteimonas, Rubrivivax, and Anaerolinea, while the nosZ-harboring network community identified the genera Microvirga, Rhizobium, Burkholderia, and Mesorhizobium based on their network centrality and closeness centrality through the OTU modules (Figures 6C, 7C).

Using random forest modeling, our study aimed to ascertain the association between N₂O emissions, maize yield, and potential predictors related to denitrifier communities containing nirS and nosZ genes. We examined biotic variables (such as composition, abundance, and diversity) and abiotic factors (including soil physiochemical properties). The results of the random forest modeling revealed that in the bulk soil, pH (5.2–6.1%, p < 0.05), SOC (8.0–8.9%, p < 0.01), NO₃⁻-N (6.9–7.5%, p < 0.05), and DON (10.5–11.5%, p < 0.01) played a crucial role as indicators of abiotic variables for N₂O emissions and maize yield (Supplementary Figures S3A,B). However, in the rhizosphere soil, the significant abiotic drivers affecting N₂O emissions and maize productivity were pH (7.7–10.3%, p < 0.01), SOC (6.3–10.1%, p < 0.01), SWC (8.7–9.3%, p < 0.05), and NO₃⁻-N (8.1–9.8%, p < 0.01) (Supplementary Figures S3C,D).

Additionally, the primary (biotic) factors influencing N₂O emissions and maize yield in the bulk soil, specifically within the denitrifier communities of nirS and nosZ genes, were found to be abundance (8.7–9.6%, p < 0.05 and 6.7–7.7%, p < 0.01), composition (3.5–6.9%, p < 0.05 and 1.2–2.1%, p > 0.01) and Module I (6.1–6.4%, p < 0.05 and 3.4–4.3%, p < 0.01), respectively (Supplementary Figures S3A,B). Conversely, in the rhizosphere soil, the significant (biotic) factors affecting N₂O emissions and maize productivity within the denitrifier communities harboring nirS and nosZ genes were identified as abundance (10.6–11.1%, p < 0.01 and 7.6–8.5%, p < 0.05), Diversity (5.5–7.4%, p < 0.05 and 2.3–5.0%, p > 0.05), Module I (4.6–6.0%, p < 0.05 and 2.5–5.2%, p > 0.05), and Module II (6.7–9.7%, p < 0.05; and 5.7–7.4%, p < 0.05), respectively (Supplementary Figures S3C,D).

We established a structural equation model to establish further connections between the potential predictors, including composition, abundance, diversity, and network Modules of nirS- and nosZ-harboring denitrifiers and the abiotic drivers represented by soil properties. The aim was to examine their impact on N₂O emissions, potential denitrification activity (PDA), maize productivity, and nitrogen use efficiency (NUE). Overall, in the bulk soil, the soil physiochemical properties (such as SOC, pH, NO₃⁻-N, and DON) exhibited significant positive effects on the denitrifier community harboring the nosZ gene through abundance (r = 0.46, p < 0.05). Similarly, the nirS-harboring denitrifier community was influenced by abundance, composition, and Module I (r = 0.81, p < 0.01). Additionally, the soil properties positively affected maize productivity (r = 0.72, p < 0.01), and PDA (r = 0.49, p < 0.05; Figure 8A). The abundance of nosZ in the denitrifier community had a significant positive impact on PDA (r = −0.39, p < 0.05) and maize productivity (r = 0.65, p < 0.01). The nirS-harboring denitrifier community, on the contrary, had a significant positive effect on PDA (r = 0.63, p < 0.05) but a significant opposite influence on maize productivity (r = −0.58, p < 0.05) via abundance, composition, and Module I. Furthermore, a statistically significant positive relationship existed between PDA and N₂O emissions (r = 0.71, p < 0.01). Maize yield also showed a positive association with nitrogen use efficiency (NUE) (r = 0.42, p < 0.05; Figure 8A).

The rhizosphere soil showed physiochemical characteristics of the soil (i.e. SOC, pH, NO₃⁻-N and SWC) had significant positive impacts on the composition of the nirS-harboring denitrifier community through abundance, diversity, Module I and Module II (r = −0.29, p < 0.05; Figure 8B). Additionally, the nosZ-harboring denitrifier community was influenced by abundance and Module V (r = 0.49, p < 0.05), while maize yield (r = 0.71, p < 0.05), and PDA (r = 0.38, p < 0.05) showed significant positive relationships with the same physiochemical characteristics (Figure 8B). Furthermore, the combined effects of abundance, diversity, Module I, and Module II had a substantial positive impact on the PDA (r = 0.59, p < 0.05) and Yield (r = 0.45, p < 0.05) in the nirS-denitrifier community (Figure 8B). The abundance and Module V revealed a significant positive impact on maize productivity (r = 0.67, p < 0.01), while PDA had a significant negative influence (r = −0.52, p < 0.05) within the nosZ-harboring denitrifier community (Figure 8B). NUE and maize productivity are positively associated (r = 0.81, p < 0.01). As shown in Figure 8B, PDA had a statistically significant positive connection (r = 0.65, p < 0.05) with N₂O emissions. The results from both the bulk and rhizosphere soil indicate that, in contrast to the nosZ-harboring denitrifier community, the nirS-harboring denitrifier community exhibited a strong positive effect with PDA activity and made a significant contribution to N₂O emissions.

To maintain the sustainability of the microbial community composition and nitrogen-cycling processes, it is crucial to examine how nitrogen fertilization and soil physiochemical dynamics affect denitrification activities in the bulk and rhizosphere soil (Yang et al., 2018; Pang et al., 2019; Schmidt et al., 2019). Significant variations were detected in the direction of the rhizosphere soil compared to the bulk soil related to the abundance and diversity of nirS and nosZ gene copy species. This finding aligns with prior research highlighting how roots and soil pH contribute to competitive filtering, increasing biodiversity, and species richness of functional genes linked to nitrogen cycling (Han et al., 2020; He et al., 2020). The levels of nitrogen application impacted the contents of SOC, TN, and DON, which contributed to the prevalence of denitrifiers carrying the nirS gene, making them more notable than the nosZ gene. This suggests their high soil nutrient cycling and impact on the functional aspects of denitrification processes (Cui et al., 2016; Ouyang et al., 2018; Huang et al., 2020b). Substantially SOC and NO₃⁻-N in optimal N fertilization emerged as a significant factor contributing to the elevation of denitrifier gene copies (Tao et al., 2018; Fudjoe et al., 2022). Various nitrogen fertilizer rates and environmental factors exerted distinct effects on the diversity of denitrifying communities associated with nirS and nosZ, underscoring their specialized enhancement linked to specific ecological niches that lead to heightened species richness (Yoshida et al., 2009; Jang et al., 2018; Shi et al., 2019).

Our findings indicate that alterations in the structure of denitrification communities were predominantly influenced by varying rates of nitrogen fertilization and soil properties, contributing to enhanced soil fertility. The carbon and nitrogen contents and their decomposition in the soil played a crucial role in shaping the structure of the denitrifying community (Huang et al., 2019; Chen et al., 2019a). Distinct differences were noted in the composition of denitrifying community genera between the bulk and rhizosphere soils. In our study, the bacteria genera Cupriavidus, Rhodanobacter, Bradyrhizobium, and Zoogloea were prominent in the composition of the nirS-harboring community in the bulk soil under the N3 and N1 fertilizer treatments. The increase in SOC and DON under different nitrogen fertilizer rate treatments underscored the significant contribution of the nirS denitrifying community to N₂O emissions (Sun et al., 2017; Tao et al., 2018). Cupriavidus possesses metabolic pathways that utilize SOC and NO₃⁻-N as alternative electron acceptors (Fudjoe et al., 2022). In alkaline soil, the genus Rhodanobacter facilitates nitrogen losses during denitrification (Ma et al., 2018). However, in the rhizosphere soil under N2 and N1 fertilizer treatments, the dominant bacteria in the nirS-harboring denitrifier community composition were genera Luteimonas, Anaerolinea, and Rubrivivax. Genus Luteimonas is recognized for its capacity to produce N₂O and significant role in nitrogen biogeochemical cycles, rendering it promising for potential applications in bioremediation (Finkmann et al., 2000; Wen et al., 2017). Genus Anaerolinea, which belongs to the phylum Chloroflexi, is an oligotrophic bacterium found in alkaline soils that can decompose various carbon substrates under anaerobic conditions (Matsuura et al., 2015).

Furthermore, the nosZ denitrifying community stimulates the conversion of N₂O to N₂, reducing N₂O emissions (Jang et al., 2018; Shi et al., 2019). The application of different rates of N fertilizers effect on the contents of DON, TN, SWC, and NO₃⁻-N significantly increased the abundance of nosZ denitrifying communities in the bulk and rhizosphere soil (Ouyang et al., 2018; Han et al., 2020; Yang et al., 2021a). In the bulk soil, the main genera within the nosZ denitrifying communities under N2 and N3 treatments were Mesorhizobium, Azospirillum, and Herbaspirillum. Specifically, genera Mesorhizobium and Azospirillum were pivotal for nitrogen fixation, playing a vital role in the production of nutrient biomolecules such as organic compounds, enzymes, adenosine phosphate groups, and polynucleotides (Huang et al., 2020a; Fudjoe et al., 2022). The usage of nitrogen fertilizers resulted in the decomposition of available nutrients such as SOC and NO₃⁻-N in the rhizosphere soil resulting in the notable prevalence of the genus Rhizobium, Burkholderia, and Microvirga among denitrifiers containing nosZ. Genus Rhizobium and Microvirga are recognized for their distinct interactions involving nitrogen fixation, antibiotic production, and phytoalexins synthesis in maize root microbiota (Sun et al., 2017; Tao et al., 2018). These findings reinforce the notion that denitrifying microorganisms residing in the rhizosphere soil play significant roles in attaining microbiome exudates and regulating N₂O emission, consequently impacting factors related to the environment (Yang et al., 2018; Schmidt et al., 2019).

Nitrogen fertilization and environmental factors are vital in increasing maize production and improving soil fertility in the semiarid Loess Plateau region. Nonetheless, this practice significantly influences the release of the emission of N₂O from the soil, which is closely linked to plant growth and soil quality elements (Ward et al., 2017; Shi et al., 2019; Yue et al., 2022). During our research, N₂O emissions were higher in soils treated with nitrogen fertilizer rates (N2 and N3), with the highest levels occurring in July. Moreover, plots without fertilization (N0) had comparatively lower N₂O emissions than those with nitrogen fertilization. The study reveals that climatic conditions are crucial in determining N₂O emission rates. The growth rate of crops, soil microbial activity, and nitrogen cycling in the soil is substantially impacted by factors such as soil temperature and moisture conditions (Peng et al., 2018; Fudjoe et al., 2022). Plastic mulch increased soil moisture and temperature (ranging from 21 to 31°C), leading to variations in drying and wetting cycles within the same period (June and July). These variations enhanced the soil water-holding capacity and temperature, improving soil aeration that affects N₂O emissions by modifying the populations of denitrifiers (Ouyang et al., 2018; Zhang et al., 2021). Earlier studies have indicated that significant soil water stress might hinder the positive effects of temperature on nutrient accessibility to microbes (Han et al., 2020). Previous research also indicated that microorganisms regulate the allocation of SOC between cell growth and stress tolerance, thereby influencing their participation in nutrient cycling in N₂O emissions (Cui et al., 2016; Jang et al., 2018).

Furthermore, employing various nitrogen fertilizer rates increased potential denitrification activity (PDA) in the bulk and rhizosphere soil compared to the absence of fertilizer treatment (N0). This variation can be attributed to the denitrification mechanisms that contribute to the reduction of nitrate (NO₃⁻-N) and changes in the physiological activity of individual cells within the soil’s denitrifying communities (Yoshida et al., 2009; Yang et al., 2018). The application of N1, N2, and N3 treatments creates a favorable environment for denitrifiers primarily driven by increased quantities of organic C and N compounds, accelerating the decomposition that elevates soil nutrient turnover. Consequently, this can augment the activities of soil microbiota, thereby impacting the production of PDA and N₂O emissions (Shi et al., 2019; Ullah et al., 2020). Our studies revealed a connection between elevated levels of nutrients such as SOC, DON, and NO₃⁻-N and the emissions of N₂O (Cui et al., 2016; He et al., 2020). The application of the N1 treatment had a relatively minor impact on soil N₂O emissions compared to N3 and N2 fertilizer treatments. However, the decomposition of nitrogen fertilizer via hydrolysis led to the release of heavily oxidized nitrate substrates (NO₃⁻-N), which subsequently caused a rise in N₂O emissions (Ouyang et al., 2018; Fudjoe et al., 2021).

Nitrogen fertilizer stimulates crop growth, increasing root exudates and activating soil organic carbon (Chen et al., 2019b). Our study observed that the N200 and N300 treatments increased yields, NUE, and cumulative N₂O emissions. However, there was no significant difference between these two treatments (N2 and N3), indicating that applying nitrogen fertilizer beyond 200 kg N ha⁻¹ did not result in further yield increase. Our findings suggest that applying nitrogen at a 200 kg N ha⁻¹ rate is advantageous in ensuring crop productivity while minimizing environmental impacts, as it achieves a high yield per unit of N₂O emissions. Earlier research has shown that the application of 200 kg N ha⁻¹ on a maize field in the semiarid Loess Plateau leads to a substantial increase in yield and nitrogen use efficiency (NUE) and also effectively limits the emission of greenhouse gases (Yue et al., 2022; Wang et al., 2023).

The significance of soil denitrifying communities in regulating the ecological nitrogen cycle and soil properties is essential (Ullah et al., 2020; Fudjoe et al., 2022). The utilization of N fertilizers influences traits of the plant–soil system, encompassing maize yield, nitrogen use efficiency (NUE), N₂O emissions, and soil properties. These elements are crucial in establishing the connection between biodiversity and ecosystem function within semiarid regions (Chen et al., 2016; Liu et al., 2018a; Han et al., 2020). The presence of denitrifiers harboring the nirS gene positively correlated with soil nutrients (soil pH, SOC, NO₃⁻-N and DON), PDA, and N₂O emissions in the bulk and rhizosphere soil, indicating the differential effects. By examining the connection between maize productivity, NUE, and soil factors (abiotic and biotic), we gained insight into the primary mechanisms behind denitrification and microbial-driven nutrient cycling in both the bulk and rhizosphere soil (Yin et al., 2014; Sun et al., 2017).

In contrast, the denitrifier containing the nosZ gene exhibited an opposite pattern. The abundance of nosZ gene denitrifier in the bulk and rhizosphere soil showed a positive relationship with maize productivity and nitrogen use efficiency (NUE) and with variables such as NO₃⁻-N, SOC, pH, SWC, and module V. The differences observed in the nosZ-denitrifying community could be attributed to the diffusive transport of organic substrates such as soil organic carbon and nitrogen nutrients involved actively in the cycling of nitrogen communities (Yin et al., 2014; He et al., 2020; Huang et al., 2020b). The prevalence of DON functions as a nitrate reducer leads to a significant increase in N₂O emissions, which aligns with previous research (Yang et al., 2018; Schmidt et al., 2019). The root exudates released by plants contain a variety of soil nutrients (NO₃⁻-N, DON, and SOC) that enrich niche distinction and functional taxonomic species of denitrifiers in the rhizosphere soil, affecting maize productivity (Chen et al., 2016; Kuang et al., 2018).

Additionally, we discovered that the network of nirS and nosZ denitrifying compositions consisted of distinct microbial modules comprising closely related species, which we identified through co-occurrence network analysis (Williams et al., 2014; Herren and McMahon, 2018). When examining the structural variances and functional processes associated with different nitrogen fertilization rates, we observed that the edges and nodes in the nosZ denitrifying network exhibited denser connections compared to the nirS-denitrifying networks in the bulk and rhizosphere soil. This discrepancy could be attributed to environmental variations and niche differentiations within the denitrifying networks (Han et al., 2020; Yang et al., 2021b). The presence of a higher proportion of positive relationships compared to negative edges in the bulk and rhizosphere soil networks indicates a mutual association that reduces competition among species within the modules of nirS- and nosZ-containing denitrifiers (Schmidt et al., 2019; He et al., 2020). The nirS-carrying denitrifier network exhibits a high level of centralization, facilitating efficient information transfer among discrete components of the soil microbiomes and highlighting their effective performance within the rhizosphere modules (Chen et al., 2019a; Han et al., 2020). Additionally, the nirS- and nosZ-containing denitrifiers had higher rhizosphere modularity compared to the bulk soil due to a more likely resilient links within the functional alignments and extensive network connectivity (Barberán et al., 2012; Berry and Widder, 2014; Banerjee et al., 2016). The presence of keystone species within the modules significantly impacts shaping the soil microbiome’s structure, function, and ecological steadiness (Mamet et al., 2019; Chen et al., 2019b).

Furthermore, the dynamics of nitrogen (N) fertilizer application stimulate crop growth, leading to increased root exudates from the microbiome network and improved soil physicochemical properties (Yue et al., 2022). This could contribute to the superior adaptability of biodiversity in the rhizosphere soil network compared to the bulk soil involved in nitrogen cycling (Amadou et al., 2020; Yang et al., 2021a). This distinction might be attributed to elevated soil pH, SOC, NO₃⁻-N, maize root exudation, biomass, and litter, all contributing to enhanced soil mineralization (Schmidt et al., 2019). Therefore, our findings imply that the augmentation in N₂O emissions could be attributed to the specific influence of soil nutrients and microbial root respiration, alongside potential competitive interactions induced by substrate from keystone taxa. These factors collectively shape the diversity of denitrification-related soil microbiomes (Ward et al., 2017; Wen et al., 2017). However, caution is necessary when making conclusions regarding the potential contribution of keystone species on denitrification communities to the overall networking system; additional stable isotope tracking is required to disentangle the underlying mechanisms further.