Limited Effects of Fertilization Practices on Microbial Community Structure and Greenhouse Gas-Related Functional Potential in Paddy and Field Soils

Introduction

Agricultural lands are major anthropogenic sources of greenhouse gases (GHGs), including methane (CH₄) and nitrous oxide (N₂O), and exert a significant influence on global climate change. The production and consumption of these GHGs are primarily regulated by microbe-mediated biogeochemical processes such as methanogenesis, methane oxidation, nitrification, and denitrification. The structural and functional characteristics of soil microbial communities act as key regulators of agricultural GHG emissions. Therefore, elucidating the impacts of agricultural management practices on soil microbial communities and their functional potential is essential for establishing sustainable agricultural strategies to mitigate climate change. Among agricultural land-use types, paddy and field soils provide distinctly contrasting ecological environments. Paddy soils are characterized by limited oxygen availability and dominant reductive conditions due to flooding management, whereas field soils are relatively aerobic and frequently exposed to mechanical disturbances like tillage. These hydrological and physicochemical differences significantly influence the composition and assembly processes of microbial communities. Numerous previous studies have reported that flooding status is the dominant factor driving microbial community differentiation, surpassing the influence of crop species [6].

Soil microbial communities can also differentiate along environmental gradients within the vertical soil profile. Generally, variations in oxygen, organic carbon, and nutrients across depths induce microbial niche differentiation. However, in agricultural lands, solute movement via irrigation, vertical infiltration, and physical mixing by tillage can mitigate this vertical heterogeneity. Indeed, an increasing number of reports suggest that differences in microbial communities according to depth are limited in intensively managed paddy and field soils [12]. Fertilization practices are also critical management factors that can influence soil microbial communities. Fertilization alters nutrient availability and stoichiometric constraints in the soil, potentially inducing increases or decreases in specific microbial taxa or functional groups. However, findings regarding the impact of fertilization on overall microbial diversity or community structure remain inconsistent across studies. Even in long-term fertilization experiments, its effects are often reported to be limited compared to overarching environmental factors such as soil moisture status or land-use type [4]. This suggests that the effects of fertilization may manifest selectively at the level of specific taxa or functional niches rather than the community as a whole. Recently, there has been increasing interest in the "decoupling" between shifts in microbial community composition and functional potential. However, it remains insufficiently understood how this phenomenon manifests across different agricultural land-use types. In particular, whether taxonomic shifts observed due to fertilization extend to actual changes in the GHG-related functional gene pool remains a significant research gap. Accordingly, this study comprehensively evaluated the effects of land-use type, soil depth, and fertilization methods on soil microbial community structure and GHG-related functional potential in paddy and field soils. Through 16S rRNA gene-based community analysis and Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2 (PICRUSt2)-based functional prediction, we aimed to verify (i) the patterns of microbial community differentiation between paddy and field soils, (ii) community stability within the vertical soil profile, (iii) changes in microbial diversity and community structure according to fertilization methods, and (iv) whether taxonomic-level changes translate into shifts in functional potential at the GHG metabolic pathway level.

Results and Discussion

Community Differentiation by Land-Use Type and Stability within Vertical Profiles

In this study, the microbial community structures of paddy and field soils were compared using β-diversity analysis, which revealed statistically significant differences between the two land-use types (Fig. 1). This indicates that the distinct soil environmental conditions shaped by land use are major determinants of microbial community composition. Previous research has reported that paddy soils, influenced by flooding management, experience more restricted oxygen availability compared to field soils, leading to distinct redox gradients and unique chemical properties (e.g., pH fluctuations) that create microbial habitats differentiated from those of aerobic upland fields [6,11]. Interestingly, within each land-use type, the analysis of community structure across soil depths showed no significant vertical differentiation in either paddy or field soils (p>0.05; Supplementary Fig. S1 and Table S1). This absence of vertical stratification suggests that the management practices employed in each agricultural system induce microbial homogenization throughout the soil profile. In paddy soils, the limited effect of depth can be explained by hydrological connectivity under irrigation management. Vertical percolation and preferential flow pathways facilitate the active downward transport of dissolved organic matter (DOC) and soluble nutrients from the surface to deeper layers [22, 9]. This mitigation of resource imbalances across depths is thought to inhibit the vertical segregation of microbial communities [12]. Similarly, the microbial community structure in field soils did not differ significantly across depth intervals (p>0.05), likely due to mechanical disturbances such as tillage and plowing, which are common in non-flooded croplands. Periodic soil inversion and mixing disrupt the development of vertical environmental gradients and physically blend the microbial populations, maintaining community composition in the subsoil similar to that of the surface layer [2]. Taken together, these findings suggest that irrigation management in paddies and mechanical tillage in fields serve as major drivers of microbial community stability across soil profiles through hydrological and physical pathways, respectively.

Differences in Microbial Diversity by Fertilization Practices

Analysis of the impact of fertilization practices on microbial diversity in both paddy and field soils revealed that α-diversity (Shannon and Chao1 indices) did not exhibit statistically significant differences among treatments (Fig. 2). In both paddy and field soils, the Kruskal–Wallis test confirmed no significant influence of fertilization practices on any diversity indices (Table 1), suggesting that the structural diversity within the microbial communities remained robustly stable under the various fertilization regimes tested in this study. β-diversity analysis (principal coordinates analysis (PCoA) based on Bray–Curtis distance) further supported this stability, as no distinct separation of community structures was observed according to fertilization in either paddy or field soils (Fig. 3). This lack of systemic shifts was statistically reinforced by ANOSIM results (Paddy: R=0.05, p>0.05; Field: R=0.00, p>0.05) and pairwise Permutational Multivariate Analysis of Variance (PERMANOVA), which consistently yielded no significant differences among treatments (p>0.05; Supplementary Table S2). These findings align with previous research indicating that long-term agricultural management can have limited effects on the overall structural assembly of soil bacterial communities [4]. In summary, the impact of fertilization practices was not manifested as systemic changes in community-wide diversity or structure. Consequently, Section 3 utilizes Linear Discriminant Analysis Effect Size (LEfSe) analysis to identify selective responses at the specific taxonomic level, which may provide deeper insights into the subtle modulations within the functional microbial niches.

Differentially Abundant Taxa Identified via LEfSe Analysis relative to Standard Fertilization

Analysis of the impact of fertilization practices using LEfSe relative to standard fertilization (std) revealed significant differences in the relative abundance of specific microbial taxa depending on the treatment (Fig. 4). In paddy soils (Fig. 4A), among the differentially detected taxa, Methanocella, which is associated with methanogenesis, significantly increased in the unfertilized treatment (none) [14]. Conversely, the methanotroph Methylocystis showed a decreasing trend in the reduced deep fertilization treatment (deep2) [1] (Fig. 4). In field soils (Fig. 4B), Nitrospira, a key taxon involved in nitrification, was significantly enriched in the deep fertilization treatment (deep1) compared to the standard practice [7]. In contrast, Rhodanobacteraceae, which is associated with the nitrogen cycle, significantly decreased in the same treatment (deep1) [10] (Fig. 4). Interestingly, the aforementioned taxa that exhibited fluctuations in response to fertilization practices are known to play pivotal roles in soil methane and nitrogen cycling —specifically, in GHG metabolism. These shifts at the individual taxonomic level suggest that fertilization practices can exert selective influences on GHG-related microbial niches. To investigate how these microscopic changes relate to the overall functional potential of the soil, the next section compares the distribution of functional genes associated with GHG metabolic pathways to evaluate consistency with the observed taxonomic shifts.

Robustness of GHG-related Functional Potential to Fertilization Practices

To evaluate the impact of fertilization practices on the functional characteristics of the soil microbiome, we first quantified the total potential—a representative index for the cumulative capacity of metabolic pathways—for four core pathways involved in CH₄ and N₂O metabolism: nitrification, denitrification, methanogenesis, and methane oxidation. In both paddy and field soils, no statistically significant differences in total potential were detected among fertilization practices (Kruskal–Wallis test, p>0.05; Fig. 5, Supplementary Table S3). Furthermore, pairwise comparisons based on Dunn’s post-hoc test with Benjamini–Hochberg correction showed that all treatments were grouped under the same letter (‘a’), statistically confirming that fertilization practices did not significantly alter the aggregate indicators at the pathway level (Fig. 5, Supplementary Table S4). Subsequently, to determine whether this macroscopic stability masked any underlying shifts at the gene level, the relative abundances of GHG-related functional genes predicted via PICRUSt2 were log-scaled and compared (Supplementary Fig. S2). Consistent with the module-level observations, the overall distribution patterns of individual functional genes remained highly similar across fertilization treatments in both paddy and field soils, despite the differential abundance of specific taxa observed in the preceding LEfSe analysis (Supplementary Fig. S2). It should be noted that PICRUSt2 provides predictions of functional gene potential rather than direct measurements of in situ metabolic activity or GHG flux. Therefore, empirical validation through gas flux measurements or metatranscriptomic approaches would further strengthen the conclusions drawn here. These findings suggest that land-use type is a primary driver of soil microbial community structure and associated GHG-related functional potential, with fertilization practices exerting comparatively limited and taxon-selective effects. From a practical standpoint, this implies that broad-scale GHG mitigation strategies in agricultural soils may be more effectively achieved through land-use management decisions—such as conversion between paddy and upland systems—rather than through optimization of fertilization regimes alone. Nevertheless, the taxon-level shifts observed in response to fertilization (e.g., Methanocella, Nitrospira) indicate that targeted management of specific niches may still offer opportunities for fine-scale modulation of GHG-cycling microorganisms. In summary, while fertilization practices can induce selective shifts at specific taxonomic levels, these changes do not consistently extend to shifts in the community-wide gene pool or module-level functional indicators. This taxonomic-functional decoupling [3], where taxonomic variation does not translate into functional differentiation, can be interpreted alongside the high level of functional redundancy within the soil microbial community, wherein multiple phylogenetically distinct taxa share overlapping metabolic capabilities, allowing other community members to compensate when specific taxa such as Methanocella or Nitrospira fluctuate in relative abundance [4]. Moreover, it suggests that higher-level environmental constraints—such as flooding conditions in paddies or the stable physico-chemical environment of field soils—may exert a stronger influence on the maintenance of microbial functional configurations than fertilization management [8]. However, since these results are based on 16S rRNA-based predictions (PICRUSt2), future studies incorporating shotgun metagenomic sequencing and transcriptomics, alongside actual environmental data (e.g., pH, C/O ratios, Organic matter), are necessary to precisely validate the correlation between predicted functional potentials and empirical GHG emission rates.

Materials and Methods

Soil sample collection

Soil samples were collected from May to August at an experimental field located in Deokjin-gu, Jeonju-si, Jeonbuk, Korea (approx. 35.835° N, 127.050° E). The experimental design consisted of six fertilization treatments: standard fertilization (std), no fertilization (none), deep fertilization (deep1), deep fertilization with reduced rates (deep2), customary fertilization (cus), and organic fertilization (org). At each treatment plot, soil samples were collected at three distinct depths: topsoil (top, 0–20 cm), middle soil (middle, 20–40 cm), and subsoil (sub, 40–60 cm). Three independent experimental plots (n=3) were sampled per treatment to ensure biological replicates. The collected samples, intended specifically for DNA extraction, were cleared of plant residues and stones and immediately transported to the laboratory in an ice box (maintained below 4℃). Subsequently, the soil was aliquoted into sterile tubes and stored at −80℃ until further processing for DNA extraction.

Total DNA extraction and 16S rRNA gene amplicon sequencing

Total DNA was extracted from the soil samples using the Exgene^TM Soil DNA Mini Kit (GeneAll Biotechnology, Seoul, Korea). The concentration of the extracted DNA was quantified using a Qubit 2.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) with the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, MA, USA). The purified DNA was stored at −80℃ until subsequent analysis. The 16S rRNA gene libraries were prepared according to Illumina’s 16S Metagenomic Sequencing Library Preparation Protocol. The V3–V4 hypervariable regions of the 16S rRNA gene were amplified using the primers 341F (5′-CCTACG GGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′) with TaKaRa EX Taq Hot Start Version (Takara Bio, Shiga, Japan). The resulting PCR products were purified using HiAccuBeads (AccuGene, Korea), followed by the attachment of Illumina index adapters (barcodes). The indexed amplicons underwent an additional purification step before final pooling. Sequencing was conducted on the Illumina MiSeq platform (2×300 bp paired-end reads) at Macrogen Inc. (Seoul, Korea).

Data processing

To ensure an equal sequencing effort across all samples, the raw sequence dataset was initially subsampled to 10,000 reads per sample. Subsequently, the reads were processed using Mothur software (v1.48.2) [5]. After merging paired-end reads, low-quality reads were removed. The resulting reads were aligned to the SILVA reference database (v138.2) [16], and chimeric sequences were identified and removed with VSEARCH (v2.17.1) [17]. Operational taxonomic units (OTUs) were clustered at 97% similarity using the OptiClust algorithm [14] and classified against the Ribosomal Database Project (RDP) 16S rRNA training set (v19) [5].

Functional prediction

Functional profiles were inferred from 16S rRNA gene sequences using PICRUSt2 (v2.6.2) [2]. Representative sequences and the sample-by-OTU abundance table generated in mothur were used as inputs. The PICRUSt2 workflow (picrust2_pipeline.py) was applied to place sequences into a reference phylogeny, predict gene family abundances at the level of KEGG Orthologs (KOs), and summarize predicted functions at the pathway/module level. Stratified outputs were additionally generated to decompose predicted functions into OTU-level contributions. For GHG–related functional analyses, KO sets associated with methanogenesis, methane oxidation, nitrification, and denitrification were curated based on KEGG module definitions (accessed on December 2025).

Statistical analysis

Normality was evaluated using the Shapiro–Wilk test in R (v4.5.2) for each comparison, and non-parametric tests were applied when normality assumptions were not met. α-diversity indices (Shannon and Chao1) were calculated in Mothur (v1.48.2). Differences in α-diversity among treatments were assessed using the Kruskal–Wallis test, and p-values were adjusted for multiple testing using the Benjamini–Hochberg false discovery rate (FDR) procedure using the FSA package (v0.10.0) in R [14]. Because depth-specific biological replication (n=3) was incomplete for some fertilization treatments, α-diversity comparisons for the fertilization × depth interaction were not performed; the final sample sizes included in each group are summarized in Table S1. β-diversity was computed using Bray–Curtis dissimilarity and visualized by PCoA, then Group differences in community composition were tested by Analysis of Similarities (ANOSIM) and PERMANOVA (perm. 999) in R using the vegan package (v2.7-2) [15]. Pairwise PERMANOVA was additionally conducted, and p-values were adjusted using the Benjamini–Hochberg FDR method. Taxonomic differential abundance was assessed using the LEfSe() function implemented in Mothur (v1.48.2). The Mothur-generated OTU abundance table and the corresponding design file containing group assignments were used as inputs, and LEfSe was executed using the default settings. Comparisons of PICRUSt2-predicted functional profiles (KO- and module/pathway-level) among fertilization treatments were performed using predicted relative abundances. For visualization, KO-level profiles were transformed using log(x+1) and displayed as a heatmap with hierarchical clustering. Module/pathway-level total potential was computed separately for paddy and field soils as a weighted index derived from PICRUSt2 outputs by multiplying taxon relative abundance by genome function counts. GHG-related KOs/modules were compared across land-use and fertilization treatments using the Kruskal–Wallis test. Post-hoc comparisons were conducted using Dunn’s test, and group labeling was summarized using compact letter display (CLD) using rcompanion package (v2.5.1) [12] in R. Statistical significance was set at p<0.05.

Data Availability: All data are available in the main text or in the Supplementary Information.

Author Contributions: YJ collected the soil sample, NJ collected the data, performed the bioinformatic analysis, conducted statistical analysis and wrote the first manuscript, JH and JO reviewed the manuscript, TU and DL designed the research and finalized the manuscript.

Notes: The authors declare no conflict of interest

Acknowledgments: This work was carried out with the support of “Research Program for Agriculture Science and Technology Development (Project No. RS-2025-02633155)” Rural Development Administration, Republic of Korea.

Additional Information:

Supplementary information The online version contains supplementary material available at https://doi.org/10.5338/KJEA.2026.45.08

Correspondence and requests for materials should be addressed to Tatsuya Unno.

Peer review information Agricultural and Environmental Sciences thanks the anonymous reviewers for their contribution to the peer review of this work.

Reprints and permissions information is available at http://www.korseaj.org

Tables & Figures

Fig. 1. Principal coordinate analysis (PCoA) of soil microbial community structure in paddy and field soils based on Bray–Curtis dissimilarity.

Fig. 2. Microbial alpha diversity indices across different fertilization practices in paddy and field soils: (A) Shannon diversity in paddy soil, (B) Chao1 richness in paddy soil, (C) Shannon diversity in field soil, (D) Chao1 richness in field soil.

Table 1. Kruskal–Wallis test results for alpha diversity indices (Chao1 and Shannon) across fertilization practices in field and paddy soils

Fig. 3. Principal coordinate analysis (PCoA) of soil microbial community structure based on Bray–Curtis dissimilarity across different fertilization practices: (A) paddy soil, (B) field soil.

Fig. 4. LEfSe analysis of differentially abundant microbial genera across various fertilization practices: (A) paddy soil, (B) field soil. Genus-level groups designated as 'Gp' (e.g., Gp5, Gp6) represent members of the phylum Acidobacteria. Only taxa with an LDA score > 2.0 are shown.

Fig. 5. Comparison of the total functional potential for key greenhouse gas (GHG) metabolic pathways across fertilization practices in field soil (A) and paddy soil (B).

References

1. 1. Bowman, JP., & Sly,LI. ((1993)). Revised taxonomy of the Methanotrophs: Description of Methylobacter gen. nov., emendation of Methylococcus, validation of Methylosinus and Methylocystis species, and a proposal that the family Methylococcaceae includes only the group I Methanotrophs..
2. 2. Caicedo, HH., Hashimoto, DA., Caicedo, JC., Pentland, A., & Pisano,GP. ((2020)). PICRUSt2 for prediction of metagenome functions.. Nature Biotechnology 38. 669 -673.
3. 3. Chen, H., Ma, K., Huang, Y., Yao, Z., & Chu,C. ((2021)). Stable soil microbial functional structure responding to biodiversity loss based on metagenomic evidences.. Frontiers in Microbiology 12. 716764.
4. 4. Chen, H., Ma, K., Lu, C., Fu, Q., Qiu, Y., Zhao, J., Huang, Y., Yang, Y., & Schadt,CW. ((2022)). Functional redundancy in soil microbial community based on metagenomics across the globe.. Frontiers in Microbiology 13. 878978.
5. 5. Cole, JR., Wang, Q., Fish, JA., Chai, B., McGarrell, DM., Sun, Y., Brown, CT., Porras-Alfaro, A., Kuske, CR., & Tiedje,JM. ((2014)). Ribosomal database project: Data and tools for high throughput rRNA analysis.. Nucleic Acids Research 42. D633 -D642.
6. 6. Dai, X., Yuan, Y., & Wang,H. ((2016)). Changes of anaerobic to aerobic conditions but not of crop type induced bulk soil microbial community variation in the initial conversion of paddy soils to drained soils.. CATENA 147. 578 -585.
7. 7. Daims, H., Lebedeva, EV., Pjevac, P., Han, P., Herbold, C., Albertsen, M., Jehmlich, N., Palatinszky, M., Vierheilig, J., & null,null. ((2015)). Complete nitrification by Nitrospira bacteria.. Nature 528. 504 -509.
8. 8. Edwards, KR., Bárta, J., Mastný, J., & Picek,T. ((2023)). Multiple environmental factors, but not nutrient addition, directly affect wet grassland soil microbial community structure: A mesocosm study.. FEMS Microbiology Ecology 99. fiad070.
9. 9. Fuhrmann, I., Maarastawi, S., Neumann, J., Amelung, W., Frindte, K., Knief, C., Lehndorff, E., Wassmann, R., & Siemens,J. ((2019)). Preferential flow pathways in paddy rice soils as hot spots for nutrient cycling.. Geoderma 337. 594 -606.
10. 10. Kostka, JE., Green, SJ., Rishishwar, L., Prakash, O., Katz, LS., Mariño-Ramírez, L., Jordan, IK., Munk, C., Ivanova, N., & null,null. ((2012)). Genome sequences for six Rhodanobacter strains, isolated from soils and the terrestrial subsurface, with variable denitrification capabilities.. Journal of Bacteriology 194. 4461 -4462.
11. 11. Lee, SA., Kim, JM., Kim, Y., Joa, J-H., Kang, S-S., Ahn, J-H., Kim, M., Song, J., & Weon,H-Y. ((2020)). Different types of agricultural land use drive distinct soil bacterial communities.. Scientific Reports 10. 17418.
12. 12. Mangiafico,S. ((2016)). Rcompanion: Functions to support extension education program evaluation. 2.5.1.
13. 13. Moulton‐Brown, CE., Feng, T., Kumar, SS., Xu, L., Dytham, C., Helgason, T., Cooper, JM., & Moir,JWB. ((2022)). Long‐term fertilization and tillage regimes have limited effects on structuring bacterial and denitrifier communities in a sandy loam UK soil.. Environmental Microbiology 24. 298 -308.
14. 14. Ogle, DH., Doll, JC., Wheeler, AP., & Dinno,A. ((2015)). FSA: Simple Fisheries Stock Assessment Methods. 0.10.0.
15. 15. Oksanen, J., Simpson, GL., Blanchet, FG., Kindt, R., Legendre, P., Minchin, PR., O’Hara, RB., Solymos, P., Stevens, MHH., & null,null. ((2001)). Vegan: Community Ecology Package. 2.7-2.
16. 16. Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., Peplies, J., & Glöckner,FO. ((2012)). The SILVA ribosomal RNA gene database project: improved data processing and web-based tools.. Nucleic Acids Research 41. D590 -D596.
17. 17. Rognes, T., Flouri, T., Nichols, B., Quince, C., & Mahé,F. ((2016)). VSEARCH: A versatile open source tool for metagenomics.. PeerJ 4. e2584.
18. 18. Sakai, S., Imachi, H., Hanada, S., Ohashi, A., Harada, H., & Kamagata,Y. ((2008)). Methanocella paludicola gen. nov., sp. nov., a methane-producing archaeon, the first isolate of the lineage “Rice Cluster I”, and proposal of the new archaeal order Methanocellales ord. nov.. International Journal of Systematic and Evolutionary Microbiology 58. 929 -936.
19. 19. Schloss, PD., Westcott, SL., Ryabin, T., Hall, JR., Hartmann, M., Hollister, EB., Lesniewski, RA., Oakley, BB., Parks, DH., & null,null. ((2009)). Introducing mothur: Open-source, platform-independent, community-supported software for describing and comparing microbial communities.. Applied and Environmental Microbiology 75. 7537 -7541.
20. 20. Shi, Y., Gahagan, AC., Morrison, MJ., Gregorich, E., Lapen, DR., & Chen,W. ((2024)). Stratified effects of tillage and crop rotations on soil microbes in carbon and nitrogen cycles at different soil depths in long-term corn, soybean, and wheat cultivation.. Microorganisms 12. 1635.
21. 21. Westcott, SL., & Schloss,PD. ((2017)). OptiClust, an improved method for assigning amplicon-based sequence data to operational taxonomic units.. mSphere 2. e00073 -17.
22. 22. Xu, J., Yang, S., Peng, S., Wei, Q., & Gao,X. ((2013)). Solubility and leaching risks of organic carbon in paddy soils as affected by irrigation managements.. The Scientific World Journal 2013. 546750.
23. 23. Yuan, C., Zhang, L., Wang, J., Teng, W., Hu, H., Shen, J., & He,J. ((2020)). Limited effects of depth (0–80 cm) on communities of archaea, bacteria and fungi in paddy soil profiles.. European Journal of Soil Science 71. 955 -966.