Dynamic root microbiome sustains soybean productivity under unbalanced fertilization

The study investigates how long-term, unbalanced fertilization (omitting N, P, or K) shapes the quantitative assembly and temporal dynamics of soybean root-associated microbiomes and how these dynamics relate to plant performance. Root microbiomes are dynamic across plant development, largely driven by root exudates, and can influence nutrient acquisition. Conventional relative microbiome profiling lacks absolute abundance information, whereas QMP can reveal microbial load changes under environmental perturbations. Agriculture depends on N, P, and K fertilizers, but inputs can disrupt beneficial plant–microbe interactions. Plants often face N and P limitation; cereals can recruit specific taxa (e.g., Oxalobacteraceae under low N), while legumes rely on rhizobial N fixation that is suppressed by mineral N. P and K can limit legume production and influence N fixation. The authors hypothesize that QMP will reveal distinctive microbiome dynamics, that lacking N, P, or K will differentially affect root-associated bacterial dynamics, and that functional adaptation of the rhizosphere microbiome under nutrient deficiency benefits soybean growth.

Prior work shows plant development structures root microbiomes, with drought reducing diversity and disrupting temporal dynamics. Intensive fertilization can reduce microbial biomass and temporal succession rates and impair beneficial plant–microbe interactions. Maize recruits Oxalobacteraceae under low N to improve performance. In legumes, nodulation-derived N fixation is reduced by N fertilization, reflecting tradeoffs in carbon use versus fixed N acquisition; P and K can limit legume performance and support N fixation. Rhizosphere microbiomes enrich functions for mineral nutrient metabolism compared to bulk soil via root exudate-mediated priming, but causal roles of specific microbes in compensating nutrient deficiencies remained unclear, motivating QMP-based, quantitative, temporal investigations.

Field trial at Daowai District, Harbin, China (Mollisol; soybean–maize–wheat rotation since 1979) included four fertilization regimes: Control (NPK), -N, -P, and -K, each with three 36 m² plots in a randomized block design separated by concrete walls. For soybean, application rates for Control were 75 kg ha⁻¹ N, 150 kg ha⁻¹ P₂O₅, and 75 kg ha⁻¹ K₂O. Bulk soils were sampled pre-sowing (2020) for chemistry (AHN, available P and K, pH, DOC, SOM). Soybean variety Heinong 84 was sown in early May. Rhizosphere soil and root endosphere samples were collected at 1, 4, 7, 14, 28, 42, 60, and 72 days after germination (three plants per plot; nine biological samples per treatment per stage). Rhizosphere soil was detached via PBS sonication and centrifugation; roots and nodules were surface sterilized and processed for endosphere microbiome. Total DNA was extracted (FastDNA SPIN Kit). QMP was implemented by spiking synthetic DNA (SynSpike) during PCR to allow absolute quantification. 16S rRNA gene (V5–V7; primers 799F/1193R) and rhizobial rpoB (rpoB1479F/rpoB1831R) amplicons were sequenced (MGISEQ-2000, SE400). Bioinformatics used QIIME2 with DADA2 to generate ASVs; taxonomy via Greengenes (16S) and rpoB reference database. Non-bacterial ASVs and Bradyrhizobiaceae ASVs in endosphere were filtered for 16S analyses. Phylogenies were built with MAFFT and FastTree2. For QMP, absolute abundances were computed without rarefaction; α- and β-diversity metrics calculated (vegan). RMP analyses used rarefaction (1030 sequences per sample) for comparison. Metagenomics (D1, D42, D72 rhizosphere samples) were sequenced (DNBSEQ-T7, PE100). Reads were filtered (KneadData), assembled (MEGAHIT), mapped (Bowtie2), ORFs predicted (Prodigal), non-redundant gene catalog built (CD-HIT), and gene abundance quantified (Salmon TPM). Functional annotation used KEGG, COG (eggNOG-mapper) and CAZy (DIAMOND). Binning used MaxBin2 and MetaBAT2; MAG quality assessed (CheckM); dereplication (dRep) yielded 171 strain-level MAGs (140 species-level). Co-occurrence networks: Core rhizosphere ASVs (relative abundance >0.1% in ≥1 sample) were used. Spearman correlations on absolute abundances (|ρ|>0.8, FDR-adjusted P<0.05) constructed stage-specific and aggregated networks per treatment; properties computed (igraph); robustness assessed; modules detected (Gephi); keystone hubs identified via within-module (Zi>2.5) and among-module connectivity (Pi>0.62). Culturing and SynComs: 1011 rhizosphere strains were isolated on LB and R2B; isolates identified (16S 27F/1541R). Core ASVs were matched to isolates (≥95% identity). From module#2 (low-nitrogen-enriched, LNE cluster), seven genera were assembled as SynCom7 (Rhodococcus, Lysobacter, Terrabacter, Arthrobacter, Phyllobacterium, Bosea, Aeromicrobium). Module#1 provided SynCtrl (Brevundimonas, Sediminibacterium, Mycobacterium, Herbaspirillum, Sphingomonas). SynCom5 (Rhodococcus, Lysobacter, Terrabacter, Arthrobacter, Phyllobacterium) excluded Bosea and Aeromicrobium to avoid nodulation confounders and match SynCtrl diversity. Greenhouse validation: Surface-sterilized soybean seeds were grown in sterilized vermiculite with nutrient solution (MFP) with or without N (NH₄NO₃). Four treatments: water control, SynCtrl (5 strains), SynCom7 (7 strains), SynCom5 (5 strains). Strain suspensions adjusted to OD 0.1 per strain; 1 mL inoculated at 4 days after planting and again one week later. At 2 and 3 weeks post first inoculation, shoot height and total biomass were measured; nitrate-N in whole plants was quantified (colorimetric salicylic acid-H₂SO₄/NaOH method). PGP traits (IAA production, ACC deaminase activity, inorganic P solubilization, nitrogenase activity via acetylene reduction) were assayed in vitro. Statistics: Normality (Shapiro–Wilk) and homoscedasticity (Bartlett) tests guided use of ANOVA with Dunnett’s/LSD or Kruskal–Wallis with Dunn’s. Linear mixed models (lme4) for drivers of α-diversity and abundance. PERMANOVA (Adonis2) for β-diversity. Temporal-decay regressions for Bray–Curtis vs temporal distance and treatment dissimilarity vs stage. Paired Wilcoxon tests for taxa, core ASVs, and functional gene differences (stage-averaged). FDR corrections applied except paired Wilcoxon. Visualizations via ggplot2 and ComplexHeatmap.

• Soil chemistry under long-term unbalanced fertilization: relative to Control, available nutrients were reduced by 35% (alkaline-hydrolyzable N to 112 mg kg⁻¹) in -N, 95% (available P to 7 mg kg⁻¹) in -P, and 61% (available K to 146 mg kg⁻¹) in -K (P<0.01). In -N, soil pH increased by 1.1 units and DOC decreased by 5.6 mg kg⁻¹ (P<0.001); SOM was stable among treatments.
• Yield and nodulation: Soybean yield showed no significant difference between Control and -N across 2017 and 2020 and slightly increased by ~9%; -P consistently reduced yield by 25–29% (P<0.05). -N increased root nodule number and diameter; -P decreased both. Rhizobial abundance (rpoB) increased in -N and decreased in -P relative to Control.
• Community structure drivers: QMP-based β-diversity was primarily driven by compartment (R²=0.241, P<0.001), then developmental stage (R²=0.114, P<0.001), and fertilization (R²=0.040, P<0.001). -N rhizosphere and bulk communities separated from other treatments; endosphere was less responsive.
• Temporal dynamics: Rhizosphere temporal turnover slope (Bray–Curtis vs Δd) was lower in -P (0.0014 d⁻¹) than Control (0.0043 d⁻¹), indicating delayed succession; endosphere turnover was lower overall (0.0024–0.0031 d⁻¹). Dissimilarity to Control over time was stable in -N (slope 0.0005 d⁻¹) but increased sharply in -P and -K.
• QMP vs RMP: 38.1–48.3% of rhizosphere ASVs showed distinct temporal patterns between QMP and RMP, evidencing biases when using relative abundance alone.
• Absolute bacterial loads: Loads increased with development, peaking at day 72 at 2.3×10¹⁰ copies g⁻¹ (rhizosphere) and 4.7×10⁹ copies g⁻¹ (endosphere) from early-stage levels of 5.8×10⁹ and 1.5×10⁹, respectively. -P reduced rhizosphere bacterial loads by 54%, 61%, and 75% at days 42, 60, and 72 (P<0.05), with little consistent effect in the endosphere.
• Taxon-specific dynamics (absolute): Actinobacteria increased 2.1× (rhizosphere) and 18.8× (endosphere) despite decreasing relative abundance; Proteobacteria increased 5.4× (rhizosphere) and 119.2× (endosphere). Bacteroidetes underwent late-stage expansion, multiplying 48.8× (rhizosphere) and 460.7× (endosphere) by day 72 from day 1.
• Fertilization effects by phylum: -N reduced Acidobacteria in bulk and rhizosphere (P<0.01); -P broadly reduced rhizosphere phyla, notably Bacteroidetes; -K had minimal effects (Saccharibacteria reduced in rhizosphere, P<0.01).
• Functional metagenomics: -N enriched N mineralization genes (ureC, URE, aspQ) and depleted N reduction (narB, narG, narH), denitrification (nirK), inorganic P solubilization (gcd, ppa), and K transport (KdpA, KtrB, KefB). -P enriched genes for inorganic P solubilization and P starvation responses (confirmed across MAGs), while depleting N mineralization, N reduction, and K transport; nitrification and denitrification were enriched. -K specifically enriched K transport gene kup with few other changes.
• Network properties: Rhizosphere co-occurrence networks showed higher complexity (average degree 22.48 in -N, 25.12 in -K vs 7.48 in Control) and greater robustness in all unbalanced treatments (P<0.001). Only -P network contained a module hub, implying group-level rather than single-taxon control in others.
• Ecological modules and SynCom validation: Three modules were identified; module#2 (Actinobacteria-rich) constituted a low-nitrogen-enriched (LNE) cluster increased in -N, while module#3 decreased. In -P, all modules’ loads were reduced; in -K, module loads resembled Control. From LNE, SynCom7 and SynCom5 promoted soybean growth in sterile substrate with or without N: at 2 weeks, SynCom5 increased shoot height by 15% (+N) and 6% (-N), and dry weight by 38% (+N) and 18% (-N) vs SynCtrl. At 3 weeks, plant nitrate-N increased by 39.5–81.5% (+N) and 36.6–51.4% (-N). SynComs exhibited higher IAA and ACC deaminase activities; two strains solubilized inorganic P; no N-fixation detected in strains, suggesting potential complementarity with rhizobia in situ.

QMP uncovered that soybean root-associated bacterial communities expand in absolute load and undergo strong succession across development, with marked late-stage increases in Bacteroidetes and Proteobacteria while Actinobacteria also accumulate despite declining relative abundance. Rhizosphere communities were more responsive to fertilization than endosphere, reflecting their interface role. P deprivation delayed rhizosphere microbiome maturation, reduced α-diversity changes, lowered temporal turnover, and suppressed bacterial loads, aligning with reduced nodulation and yield; yet functional enrichment for P acquisition (inorganic P solubilization and P starvation genes) suggests microbiome adaptation to scarcity that may mitigate impacts. In contrast, long-term N omission produced an early divergence in rhizosphere composition but maintained yield, consistent with increased nodulation and functional enrichment for N mineralization, indicating alternative N provisioning pathways (biological N fixation plus organic N mineralization). Network analyses indicated increased complexity and robustness under -N and -K, consistent with intensified metabolic cooperation under nutrient stress, possibly enhancing resilience. The LNE ecological cluster (module#2) overlapped with ASVs increased under -N and, when assembled into SynComs, directly promoted soybean growth and improved plant nitrate status irrespective of external N, via multiple PGP traits (IAA, ACC deaminase and some P solubilization). These quantitative and functional insights connect fertilization regimes, microbiome dynamics, and plant performance, supporting the concept that specific ecological clusters within the root microbiome can sustain legume productivity under reduced N inputs.

This work establishes a quantitative, temporal framework (QMP) for profiling root-associated microbiomes and demonstrates that soybean rhizosphere bacteria exhibit strong succession and load increases across development. Long-term unbalanced fertilization differentially reshapes these dynamics: P deprivation delays microbiome maturation and reduces loads, while N omission leads to distinct yet productive microbiomes with enhanced N mineralization signatures and nodulation. Network analyses identified an Actinobacteria-rich, low-nitrogen-enriched ecological cluster whose SynComs consistently promoted soybean growth and nitrate status without added N, highlighting a tractable microbial resource for sustainable agriculture. Future research should extend beyond bacteria to include fungi and protists and dissect tripartite interactions (e.g., rhizobia–mycorrhiza–plant) and field-scale SynCom applications across soils and climates to translate ecological clusters into robust bioinoculants.

The study focuses on bacterial components of the root-associated microbiome; other organisms (e.g., filamentous fungi and protists) likely contribute to network structure and plant-beneficial functions through top-down effects and tripartite mutualisms, which require further investigation. In SynCom experiments, no nitrogen fixation activity was detected among the strains and nodulation was not observed, indicating that interactions with rhizobia and broader community members were not evaluated in the synthetic systems.