The field trial was conducted at the First Industrial Park Experimental Station (41°43′40”N, 111°36′6″E, and elevation: 1,417.2 m) of Xinyu Seed Industry Co., Ltd. in Wulanhua Town, Siziwang Banner, Inner Mongolia Autonomous Region. This region features a mid-temperate continental monsoon climate with annual precipitation of 313.8 mm and 3,084–3,286 h of sunshine. The soil is classified as Kastanozem (chestnut soil), with a sandy loam texture, deep profile development, and loose structure. The soil exhibits a pH of 7.8, and an organic matter content of 11.12 g/kg. These soil characteristics make it ideal for potato production. During the 2024 growing season, the average temperature was 14°C and wind speed ranged from 3–4 on the Beaufort scale.

The field experiment adopted a single-ridge, unmulched cultivation system with drip irrigation, using ridges spaced 90 cm apart and plants 16 cm apart within each 6-meter ridge. Each plot contained five ridges, yielding a planting density of 69,444 plants/ha. Four treatments replaced 0% (CK), 40% (T40), 60% (T60), and 80% (T80) of synthetic nitrogen with organic fertilizer (OF; 7.06% N, 1.15% P₂O₅, 1.14% K₂O), as detailed in Table 1. Synthetic fertilizers included urea (46% N), diammonium phosphate (DAP; 18% N, 46% P₂O₅), and potassium sulfate (SOP; 52% K₂O). To maintain constant nutrient inputs (250 kg N, 200 kg P₂O₅, 300 kg K₂O/ha), OF nutrient contributions were calculated first. Phosphorus deficits were initially supplemented with DAP, and remaining nitrogen and potassium with urea and SOP. For T80, DAP was applied only until total nitrogen reached 250 kg/ha, with further phosphorus from monopotassium phosphate (PDP; 52% P₂O₅, 34% K₂O), and remaining potassium from SOP. The potato cultivar V7 was planted mechanically and fertilized manually. The experiment followed a randomized complete block design with three replicates per treatment.

Rhizosphere soil samples were collected during three critical stages of potato development: seedling (June 24, 2024), tuber initiation (July 15, 2024), and tuber bulking (August 15, 2024). For each of the four treatments, three potato plants exhibiting similar growth vigor were randomly selected and excavated using a spade to approximately 20 cm depth and 5–10 cm distance from the stem to ensure complete root system recovery with adhering soil. After removing loosely attached soil aggregates and excising aboveground plant material, the intact root systems with closely associated rhizosphere soil were immediately preserved on dry ice, yielding 12 total samples across all treatments. In the laboratory, soil adhering to root surfaces was gently removed using a fine brush and collected onto sterile paper, with the three replicate samples from each treatment thoroughly homogenized and subdivided into three equal portions (5–6 g each) before being placed in separate sterile plastic bags with appropriate labeling (Morgan et al., 2018). These samples were then stored at −80°C for the extraction of soil DNA and metabolite analysis (Cindy et al., 2012).

The total DNA was extracted from approximately 0.5 g of rhizosphere soil sample using the E. Z. N. A.® Soil DNA Kit (Omega Biotek, D5625-01). DNA purity and concentration were evaluated using a NanoDrop 2000 spectrophotometer and 2% agarose gel electrophoresis. The V3-V4 regions of bacterial 16S rRNA genes were amplified via PCR using primers 341F (5’-CCTACGGGNGGCW GCAG-3′) and 806R (5’-GGACTACHVGGGTWTCTAA T-3′), while the fungal ITS1 region was amplified with primers ITS1F (5’-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS1R (5’-GCTGCG TTCTTCATCG ATGC-3′). Amplicons (400–450 bp for bacteria, 310 bp for fungi) were purified via agarose gel electrophoresis. Sequencing libraries were prepared using the Illumina TruSeq Nano DNA LT Kit, incorporating dual-index adapters for multiplexed analysis.

Paired-end sequencing (2 × 250 bp) was conducted using the Illumina NovaSeq 6,000 platform, generating raw 16S rRNA and ITS reads for bacterial and fungal communities, respectively. Demultiplexing and merging of the raw reads were performed using FLASH (v1.2.11) with a 10% mismatch tolerance, resulting in consensus sequences of 400–450 bp for bacteria and approximately 310 bp for fungi. Following initial processing, sequences underwent quality filtering with fastp (v0.19.6) and further merging with FLASH (v1.2.7). High-quality sequences were subsequently de-noised using the DADA2 plugin within the QIIME2 (v2021.4) pipeline, which resolves amplicon sequence variants (ASVs) at single-nucleotide resolution based on sample-specific error profiles. Taxonomic classification of bacterial ASVs was performed against the SILVA 16S rRNA database (v138) using the RDP classifier (v2.2), while fungal sequences were annotated using the UNITE v9.0 database.

Microbial diversity analyses were conducted at a uniform sequencing depth. Alpha diversity metrics, including observed ASVs, Chao1 richness, and Shannon index, were calculated using Mothur (v1.30.1), while beta diversity was assessed based on weighted UniFrac distances using Bray-Curtis-based PCoA ordination and PERMANOVA validation. The top 10 most abundant genera were profiled to elucidate treatment-induced shifts in community composition. Multivariate analyses, including ANOVA and PERMANOVA, were employed to evaluate the effects of experimental conditions on dominant phyla and genera. All sequencing and bioinformatic analyses were performed in triplicate to ensure reproducibility.

Potato harvesting was carried out on October 2, 2024, at physiological maturity. From each plot, three uniformly growing potato plants were randomly selected and excavated, yielding a total of nine plants across the three plots. To ensure data integrity, care was taken to minimize tuber damage during excavation, and adhering soil was gently removed from the tubers. For each plant, both the number and total weight of tubers were measured and recorded. The total potato yield was calculated, with total yield expressed as the average tuber weight per plant multiplied by the number of plants in the plot (kg/m²).

The polynomial regression analysis examining the relationship between potato yield (Y) and organic fertilizer application (X) was conducted using Python, employing the model formulation Y = β₀ + β₁ X + β₂ X² + ℇ. Model parameters (β₀, β₁, β₂) were estimated via ordinary least squares (OLS) using the “numpy.polyfit” function, with model fit assessed through the coefficient of determination (R²). The overall statistical significance of the regression model was evaluated using analysis of variance (ANOVA) with F-tests, where p < 0.05 served as the threshold for rejecting the null hypothesis. Correlation analyses utilized Pearson’s correlation for normally distributed variables (confirmed by Shapiro–Wilk tests with p > 0.05) and Spearman’s rank correlation for non-parametric data, with all correlation matrices subjected to Benjamini-Hochberg false discovery rate (FDR) correction (q < 0.1) to adjust for multiple comparisons. The relationships between microbial community and potato yield were further analyzed through SparCC cooccurrence networks, applying thresholds of∣ρ∣ > 0.6 and p < 0.01 for network edge inclusion.

High-throughput sequencing of bacterial 16S rDNA (V4 region) and fungal ITS1 regions across 12 rhizosphere soil samples generated robust datasets for microbial community analysis. Bacterial sequencing with 41F-806R primers produced 1,470,323 raw reads, yielding 1,216,234 high-quality non-chimeric ASVs after filtering and denoising, with an average of 101,352 ASVs per soil sample (Supplementary Table S1). Fungal sequencing with ITS1F-ITS1R primers generated 1,717,472 raw reads, resulting in 1,617,588 non-chimeric ASVs, averaging 134,799 ASVs per soil sample (Supplementary Table S2). These results provide reliable datasets for downstream diversity and community structure analyses.

Alpha diversity of potato rhizosphere microbial communities was assessed across seedling, tuber initiation, and tuber expansion stages under different organic fertilizer treatments (Table 2). Good’s Coverage values exceeded 0.995 for bacteria and 0.999 for fungi across all treatments and growth stages, confirming adequate sequencing depth for reliable diversity assessment.

Bacterial richness indices (Chao1 and ACE) varied temporally, peaking during tuber expansion, followed by the seedling stage, with lowest values at tuber initiation. Within each growth period, bacterial richness increased progressively with organic fertilizer application rates (T80 > T60 > T40 > CK). This positive relationship between bacterial richness and fertilizer application remained consistent across all developmental stages. Shannon diversity indices remained stable across treatments and growth stages, indicating that fertilization affected species richness but not overall community evenness.

Fungal communities exhibited distinct patterns compared with bacterial assemblages. During the seedling stage, T60 showed highest fungal richness, followed by T80, while T40 had lower diversity. At tuber initiation, T60 maintained peak fungal richness, followed by T40, while T80 showed the lowest values. By tuber expansion, T80 exhibited substantially increased fungal richness, reaching the highest observed values. T60 maintained consistently high fungal richness throughout all developmental stages.

Furthermore, Fungal Shannon indices were consistently lower than bacterial indices across all growth stages, approximating half the bacterial community values. CK and T40 maintained lower diversity indices, while T60 and T80 showed significantly higher values (T80 > T60 > T40 > CK).

Principal Coordinate Analysis (PCoA) based on ASV composition at the genus level revealed distinct compositional patterns in bacterial and fungal communities across organic fertilization treatments and developmental stages. Sample clusters separated according to fertilization regimes, demonstrating pronounced effects of organic fertilizer application on microbial community structure. For bacterial communities, PCoA1 and PCoA2 explained 49.55 and 15.29% of total variance, respectively (Figure 1A). Fungal community ordination showed PCoA1 and PCoA2 explaining 32.49 and 19.17% of variance, respectively (Figure 1B). These high variance explanations indicate that the ordination effectively captures microbial community differentiation patterns.

Microbial communities exhibited stage-specific temporal dynamics in response to organic fertilization. During the seedling stage, bacterial and fungal communities showed pronounced dispersion across fertilization treatments, indicating differential establishment of rhizosphere assemblages. Bacterial communities progressively converged through developmental stages, clustering cohesively during tuber expansion in the upper ordination space (Figure 1A). The control treatment (CK) maintained distinctly segregated positions throughout all stages, particularly in bacterial communities where CK samples consistently occupied the negative PCoA1 region. Fungal communities exhibited greater heterogeneity than bacterial assemblages, especially during tuberogenesis, suggesting taxon-specific sensitivities to organic inputs (Figure 1B).

The taxonomic composition of bacterial and fungal communities in the potato rhizosphere exhibited distinct temporal dynamics across different growth stages and Organic fertilizer treatments. Bacterial communities were dominated by four phyla: Proteobacteria, Acidobacteriota, Gemmatimonadota, and Bacteroidota. Proteobacteria dominated both CK and T80 during the seedling stage but declined progressively through tuber formation and expansion stages. Conversely, Planctomycetota increased markedly during tuber expansion (Figure 2A). At the genus level, Sphingomonas, Vicinamibacteraceae, Flavobacterium, and RB41 were the most abundant genera. The relative abundances of Sphingomonas, Pseudomonas, and Massilia decreased during potato development, while Vicinamibacteraceae increased, particularly under T60 and T80 during the tuber expansion. Flavobacterium and RB41 remained stable across developmental stages and treatments (Figure 2B).

Fungal communities consisted primarily of Ascomycota, Mortierellomycota, Basidiomycota, and Olpidimycota, collectively accounting for ~95% of total abundance. Ascomycota alone represented 75–90% of fungal communities. During potato development, Olpidimycota decreased sharply while Ascomycota and Basidiomycota increased gradually. Glomeromycota emerged prominently in T60 and T80 during tuber expansion. Ascomycota abundance followed the pattern CK > T40 > T60 > T80 at both seedling and expansion stages, while Mortierellomycota remained unchanged (Figure 2C). At the genus level, Plectosphaerella, Fusarium, Cephalotrichum, and Mortierella dominated. Increasing organic fertilizer application correlated with declining Cephalotrichum and increasing Fusarium and Plectosphaerella during the seedling stage. This trend continued through tuber formation, with Cephalotrichum declining and Fusarium increasing. During tuber expansion, both Plectosphaerella and Fusarium increased further, ranking T40 > T60 > CK > T80. The T80 consistently enriched Scutellinia throughout the growth period (Figure 2D).

LEfSe analysis revealed stage-specific shifts in potato rhizosphere microbial communities under varying organic fertilization regimes. During the seedling stage, organic treatments (T40, T60, and T80) showed reduced bacterial enrichment compared with the control (CK). The CK was dominated by f_Sphingomonadaceae, o_Sphingomonadales, and g_Massilia. Organic amendments preferentially selected for p_Chloroflexi and p_Gemmatimonadota phyla, likely due to fertilizer-induced increases in soil polysaccharides. These conditions favored cellulose-degrading genera such as g_Levilinea and g_Roseiflexus (Figure 3A). At tuber formation, T60 showed limited bacterial enrichment, with only g_Flavobacterium and g_Sphingomonas reaching significance (LDA > 3.5). In contrast, T80 promoted p_Actinobacteriota and f_Pseudomonadaceae lineages. These included nitrogen-cycling o_Micrococcaceae and antibiotic-producing g_Pseudomonas, which aligned with increased nitrogen demand and improved soil structure (Figure 3B). During tuber expansion, T60 enriched p_Acidobacteriota, o_Vicinamibacterales, and f_Sphingomonadaceae. Meanwhile, T80 selection pressure favored specialized p_Chloroflexi taxa adapted to high organic carbon environments (Figure 3C).

Fungal communities displayed treatment-dependent functional succession throughout potato development. Each organic amendment rate shaped distinct microbial assemblages at different growth stages. At the seedling stage, T40 enhanced plant growth-promoting fungi, including Mortierella and Mortierellomycetes (Figure 3D). During tuber formation, T60 favored ecological specialists such as Fusarium equiseti and Ascomycota. T80 led to enrichment of saprophytic and weakly pathogenic taxa, including Volutella and Nectriaceae (Figure 3E). By the expansion phase, T60 promoted phosphorus acquisition through Glomeromycota mycorrhizae and maintained Mortierella-driven decomposition of organic intermediates (Figure 3F). Collectively, these results demonstrate that T60 achieves an optimal balance of microbial community structure and function, fostering beneficial fungi that support nutrient cycling and potato growth across developmental stages.

Moderate organic fertilizer application significantly enhanced potato yield, exhibiting a quadratic response whereby yield increased with organic fertilizer substitution up to 60% and then declined at higher rates. Specifically, T60 produced the highest yield (64,886 kg/ha), which was markedly greater than the control (49,824 kg/ha), while T40 and T80 yielded 60,507 kg/ha and 58,901 kg/ha, respectively. Although these three treatments did not differ significantly from one another, all outperformed the control (Figure 4).

Regression analysis revealed a quadratic relationship between potato yield and organic fertilizer application using Python’s statsmodels package. Yield increased with fertilizer rate until reaching an optimal point, then declined at higher rates. The model predicted maximum yield (81,020 kg/ha) at 51.25% organic fertilizer substitution (Figure 5). Beyond this threshold, yields declined progressively with further increases in fertilizer application, following a quadratic pattern where greater deviations from the optimum resulted in steeper yield reductions.

Microbial diversity exhibited similar patterns. Bacterial richness (Chao1) peaked at approximately 60% fertilization. Bacterial diversity correlated strongly with yield, likely through enhanced enzyme activity that improves nutrient cycling. Fungal diversity (Shannon) also demonstrated a quadratic relationship with yield, with optimal yield occurring at Shannon index values of 6.60. This suggests a balance between beneficial fungal symbiosis and risks from competition or disease at higher fungal diversity levels. The optimal fertilizer application rate for both yield and microbial diversity ranged between 55–60%. At this level, bacterial nitrifiers provided most plant-available nitrogen while fungal hyphae delivered acquired phosphorus to roots.