A six-site network experiment was established in October 2021 in the Shajiang Black Soil (Vertisol) region of the North China Plain (Figure 1). The experimental sites were located in Funan (FN), Yongxing (YX), Kantong (KT), Guoyang (GY), Mengcheng (MC), and Wuhe (WH), all within the central distribution area of Shajiang Black Soil. To reduce experimental heterogeneity, all sites were selected on level terrain. The region has a warm temperate monsoon climate, with mean annual temperatures ranging from 14.9 to 16.5 °C. Mean annual precipitation varied across sites from 823.9 to 985.2 mm. All sites were rain-fed, with no irrigation applied. Although all sites were classified as Vertisols (Staff, 2010), their initial soil physicochemical properties and textures differed substantially. Detailed site-specific characteristics are provided in Table 1.

The experiment was conducted using a randomized complete block design with three replicates and four treatments: no fertilizer (Control), pure chemical fertilizer (NPK), and chemical fertilizer partially substituted with 15% (15%M) or 30% (30%M) pig manure. The control treatment received no fertilizer inputs-neither chemical nor organic-after the experiment began, although all sites had previously followed conventional local fertilization practices. For the NPK treatment, application rates were calculated based on initial soil fertility, crop nutrient demand, and nutrient balance using baseline data (Table 1). To ensure consistency, an average rate derived from all six sites was applied uniformly.

In the 15%M and 30%M treatments, 15% and 30% of the chemical nitrogen (N) supplied by pig manure, following the principle of nitrogen equivalence. All fertilizer inputs were standardized and reported as pure N, P₂O₅, K₂O, and pig manure amounts (Table 2). Each year, the composted pig manure applied to wheat was identical to that applied to maize and was sourced from the same farm to ensure consistency in manure quality throughout the three-year experiment. The organic carbon and total nitrogen contents of the composted pig manure were 19.9% and 1.4% in 2021, 23.9% and 2.1% in 2022, and 21.1% and 1.7% in 2023, respectively. Manure application rates were based on dry weight, with a consistent moisture content of approximately 12% each year.

Plot sizes were at least 160 m² (20 m × 8 m), with larger dimensions where required for machinery access. After maize harvest, soils were rotary-tilled before wheat sowing in October, whereas maize was directly sown following wheat harvest in June. Maize was planted at a density of 0.25 m × 0.6 m (intra- and inter-row spacing, respectively). All organic fertilizer, phosphorus, potassium, and 60% of nitrogen were applied as basal dressing, while the remaining 40% of nitrogen was top-dressed at the jointing stage of wheat and the large-flower-bud stage of maize in each season.

During crop maturity (2021 to 2023), undisturbed soil cores (5.0 cm diameter and 5.0 cm height) were collected from each plot at 7.5–12.5 cm depth and oven-dried at 105 °C to determine bulk density (BD), representing the 0–20 cm soil layer. Composite soil samples (~1 kg) were collected from the 0–20 cm layer using a five-point sampling method. Soil properties: pH was measured in a 1:5 soil-to-water suspension. Soil organic matter (SOM) was determined using the H₂SO₄–K₂Cr₂O₇ oxidation method. Total nitrogen (TN) and available nitrogen (A.N) were measured using a Kjeldahl nitrogen analyzer with standard and semi-micro digestion methods. Total phosphorus (TP) was determined by HClO₄; digestion and molybdate colorimetry, while available phosphorus (A.P) was extracted with 0.5 M NaHCO₃ and measured by UV spectrophotometer. Total potassium (TK) and available potassium (A.K) were analyzed using flame photometry. Soil texture was determined by laser diffraction particle size analyzer. All analytical procedures followed (Lu, 1999). Elemental stoichiometric: Elemental ratios (C:N, C:P, C:K, N:P, N:K, and P:K) were calculated as molar ratios (Liu et al., 2017). Maize yield was determined from the central 10 m × 1.2 m area (two rows) of each plot, standardized to 14% grain moisture (Jiang et al., 2025). Wheat yield was assessed from a 1 m × 1 m quadrat per plot and adjusted to 13% grain moisture.

Soil quality was assessed using the soil quality index (SQI), which was derived from a minimum data set (MDS) method (Rinot et al., 2019). All indicators were first standardized, after which the MDS was identified using principal component analysis (PCA) combined with correlation analysis of nine soil variables. Five principal components (PCs) with eigenvalues greater than 1 were retained (Table 3), collectively accounting for 92.55% of the total variance. Within each PC, variables with absolute loadings within 10% of the maximum and ≥0.5 were selected. Where multiple indicators were retained within a single PC, those with lower weights and correlation coefficients exceeding 0.6 were excluded to avoid redundancy (Li et al., 2021a). Consequently, SOM, pH, BD, TP, and a.K were selected for inclusion in the MDS (Table 3).

Following MDS selection, indicator weights were determined based on the PCA using Equation 1, and all indicators were then normalized to 0 and 1 using non-linear scoring functions shown in Equation 2 and Equation 3:

Where w_i is the weight of the i-th indicator, S_i is its normalized score, and m is the number of retained principal components. γ_k and a_kj represent the eigenvalue and loading of the k-th principal component, respectively. Equation 2 was used when higher values indicated better soil quality; Equation 3 when they indicated poorer quality. The final SQI was computed as the weighted sum of all normalized scores with Equation 4.

Crop yield sustainability index (SYI) was used to evaluate the stability and sustainability of grain yield within the wheat-maize system, which was calculated as Equation 5:

Where Ŷ is the mean yield of wheat or maize, σ is the standard deviation of yield across seasons, and Y_max is the maximum observed yield.

Data were processed in Excel 2013 and visualized using OriginPro 9.0. ANOVA test was performed in SPSS 19.0. Linear regression assessed relationships among variables. To further quantify the relative importance of soil physical properties, nutrients, soil quality, and ecological stoichiometry for effect on wheat and maize yields, a random-forest approach was used by using the “ randomForest “ packages (Number of trees in each forest = 1000, Number of observations in trees’ terminal note = 2, and Number of mtry features = 4) in R (v3.6.3) software. The Gini coefficient computed from the random-forest structure was calculated to indicate the relative importance (%) of each explanatory variable (Islam et al., 2022).

Structural equation modeling (SEM) was subsequently performed in R using the lavaan package to evaluate the effects of multiple factors on wheat and maize yields and their stability. Six latent variables—namely bulk density (BD), pH, soil organic matter (SOM), nutrient content (total N, P, and K), nutrient availability (available N, P, and K), and nutrient ratios (C:N, C:P, C:K, N:P, N:K, and P:K)—were specified as drivers of soil quality (SQI), which subsequently influenced wheat and maize yields, and ultimately their yield stability (SYI). Prior to modeling, all data were tested for normality and standardized to a mean of 0 and a standard deviation of 1 (Jiang et al., 2024).

Fertilization significantly affected soil bulk density (BD), with consistent reductions under partial organic substitution relative to both control and NPK treatments. Across sites, BD were significantly declined by 1.5% to 6.3% under 15%M and 1.3% to 17.1% under 30%M, with the strongest decreases at GY and MC (P < 0.05; Figures 2a, b). Soil pH remained stable across treatments (5.1 to 7.4), within the optimal range for wheat-maize systems (Figures 2c, d). SOC increased under all fertilized treatments compared with the control, with greater gains under partial organic substitution (4.6%–22.5%) than NPK (1.5%–13.7%). Relative to NPK, with consistent patterns across wheat and maize seasons Compared to NPK, SOC significantly rose further by 2.8%–5.9% in wheat and 0.77%–8.9% in maize, particularly at KT and GY (P < 0.05; Figures 2e, f). SOC accumulation was generally higher under the 30%M than 15%M, albeit with greater variability (4.6%–17.5% for 15%M; 3.7%–22.5% for 30%M) (P < 0.05; Figures 2e, f).

Fertilizer addition increased soil total and available nutrients compared with control, though effects differed by nutrient type, substitution rate, and site (Figures 3, 4). Furthermore, nutrient concentrations were often higher under the 15%M treatment than under 30%M (Figures 3, 4). Total N and P rose by 6.1%–23.7% in wheat and 3.6%–14.9% in maize, with the highest values under NPK, followed by 15%M and then 30%M (P < 0.05);. Total K showed no significant treatment effect. Available nutrients responded more strongly than total pools. Available N rose by 7.7%–42.5% in wheat and 0.9%–50.2% in maize, while available P and K increased by 7.1%–45.9% and 0.4%–33.2%, respectively (P < 0.05). Across sites, 15%M generally outperformed 30%M for available nutrient supply, particularly at YX, KT, and GY, whereas 30%M occasionally fell below NPK, suggesting potential nutrient constraints at higher substitution levels (Figure 4).

Elemental ratios were markedly influenced by fertilization. Under the control, C:N ratios were lowest (8.98–9.81) at MC, GY, and KT, but increased under fertilized treatments to 10 to 12, whereas sites with already high C:N (FN, YX, WH) showed slight declines (Table 4). Compared with the control, C:N ratios increased under NPK, 15%M, and 30%M in KT, GY, and MC, where baseline ratios were low, but decreased in FN, YX, and WH, where initial ratios were relatively high (Table 4). C:P and N:P ratios indicated excessive P availability (17.2–29.8 and 1.66–3.00 respectively). Both NPK and organic substitution increased these ratios, improving stoichiometric balance at most sites. Similarly, C:K, N:K, and P:K ratios rose under fertilization, with stronger effects under organic substitution-particularly 30%M. However, declines under 30%M at YX suggested localized nutrient imbalance.

Soil quality index (SQI) increased markedly under fertilization, with the greatest improvements under partial substitution (Figure 5; Table 5). Relative to the control, SQI significantly rose by 230%–249% under 15%M and 278%–416% under 30%M, exceeding the gains from NPK (136%–198%) (P < 0.05). Enhancements were greater in maize than wheat, following the order: 30%M > 15%M > NPK. Crop yields were extremely low under the control (2.4–2.5 Mg ha⁻¹ for wheat, 2.8–3.2 Mg ha⁻¹ for maize). All fertilized treatments significantly improved yields, with partial substitution frequently outperforming NPK. Wheat and maize yields increased by up to 251% and 194%, respectively, with site-specific differences: 30%M was generally superior, except at KT and MC where 15%M yielded better performance. Yield stability, as measured by the sustainability yield index (SYI), improved under all fertilization regimes. On average, SYI rose by 13.9% under 30%M, 9.4% under NPK, and 8.9% under 15%M (P < 0.05).

Correlation analysis revealed strong negative effects of BD on SQI (r = –0.60 in wheat, –0.43 in maize). SOC correlated positively with TN, available N and P, and C:P and N:P ratios (r = 0.46–0.89). Soil pH was strongly associated with K availability and K-related stoichiometry, while total nutrient contents correlated closely with their available forms (Figure 6). Random forest models identified SQI as the dominant predictor of yield (26.4% for wheat, 30.5% for maize), followed by available N (17%–22%). Additional contributors included C:K ratio, total K, soil pH, and N:K ratio for wheat, and total N, P, K, and available K for maize (Figure 7a). SEM analysis confirmed that partial organic substitution improved elemental ratios, which enhanced SQI and in turn promoted yield and stability. Key pathways involved nutrient availability, total nutrient content, BD reduction, and SOC accumulation, with nutrient-related effects stronger than those of SOC (Figure 7b).

Across all six Vertisol sites, partial substitution of chemical fertilizer with organic manure (15%M and 30%M) significantly reduced bulk density, with the strongest effects at GY and MC, where initial compaction and clay content were greatest. These reductions were closely linked to SOC increases, as organic inputs promoted aggregate formation, enhanced pore connectivity, and reduced compaction (Fang et al., 2021). In initially, as fresh organic fertilizers decompose, microorganisms (e.g., fungi and actinomycetes) secrete polysaccharides and mucilage that bind fine soil particles into microaggregates (Wen et al., 2021). As decomposition progresses, stable humic substances (e.g., humic and fulvic acids) form and, together with clay minerals and cations (e.g., Ca²⁺), promote the formation of water-stable macroaggregates (Wen et al., 2021). This hierarchical pore structure enhances soil porosity and reduces bulk density. Similar long-term evidence shows that sustained organic amendments can decrease bulk density by 2.29–13.74% and increase SOC by 15.7–132.7% in Vertisols (Guo et al., 2019), reinforcing the structural benefits of organic inputs in fine textured soils.

Nutrient dynamics also reflected the interactive roles of chemical and organic fertilization. While total N, P, and K contents under partial substitution remained comparable to NPK, available nutrient pools were consistently higher (Figure 4). This finding is particularly significant, as more than 4 million hectares of Vertisols in the Huang–Huai–Hai Plain have long experienced severe soil nutrient imbalances, characterized by nutrient depletion, low nutrient availability, slow accumulation of available nutrients, and deficiencies in phosphorus and nitrogen, as well as micronutrients such as zinc, boron, and molybdenum (Xiong et al., 2020; Li et al., 2011). The implementation of an integrated organic–inorganic fertilization strategy represents an effective approach to concurrently ameliorate macronutrient (P and N) limitations and their availability in this region. This enhancement can be attributed to several processes: (i) competition between organic carbon and nutrient ions (NH₄⁺, PO₄³⁻, K⁺) for mineral sorption sites, which reduces fixation (Jilling et al., 2021; Xu and Tsang, 2024); (ii) organic acid-mediated mobilization of native nutrients; and (iii) stimulation of microbial activity, enzymatic processes, and nutrient turnover (Wu et al., 2006; Beermann et al., 2015). These mechanisms collectively increase nutrient bioavailability, particularly at YX, KT, and GY where organic substitution showed the greatest effects. Thus, partial substitution enhanced not only soil physical conditions but also nutrient accessibility beyond that achieved by chemical fertilization alone (Singh et al., 2024).

Stoichiometric ratios provide integrative indicators of soil nutrient balance and plant–microbe interactions (Liu et al., 2017). In this study, C:N ratios were below the optimal 10–12:1 in MC, GY, and KT under unfertilized conditions (Bing et al., 2015), but increased to optimal ranges under fertilization, driven primarily by SOC accumulation. Vertisols, inherently low in nutrients, suffer from a severe shortage of exogenous organic inputs in the absence of organic fertilizer application (Xiong et al., 2020). Continuous carbon removal through grain and straw harvest exacerbates this deficit, driving the soil carbon cycle into a persistent imbalance (Wang et al., 2021). Under these conditions, soil microorganisms accelerate the decomposition of the limited soil organic matter (SOM) to obtain energy (Huang et al., 2018). However, the lack of fresh carbon inputs restricts replenishment of the soil organic carbon (SOC) pool, leading to slow accumulation or continual decline (Huang et al., 2018). Meanwhile, microbial decomposition of SOM releases mineral nitrogen, resulting in rapid carbon depletion coupled with relative nitrogen accumulation and maintaining a low soil C:N ratio (typically below 10:1; Table 3). In contrast, manure application supplies substantial amounts of organic carbon and nitrogen, stimulating microbial activity—particularly bacterial proliferation—and promoting SOC accumulation (Wu et al., 2006). The increase in SOC arises both from the carbon contained in the manure and from microbial residues formed during decomposition (Zhang et al., 2024), thereby enhancing SOC stability and restoring the soil C:N ratio to an optimal range.

Meanwhile, C:P and N:P ratios revealed a regional excess of phosphorus, consistent with surveys across the North China Plain (Xiong et al., 2021). Both chemical and organic fertilization increased these ratios by elevating soil C and N while P remained relatively stable, thereby alleviating P-induced stoichiometric imbalance. This finding provides a cost-effective strategy for medium- to low-yield lime concretion Vertisols through organic fertilizer substitution aimed at “increasing carbon and supplementing nitrogen.” By enhancing soil C and N pools, this approach corrects phosphorus-driven imbalances without resorting to costly phosphorus-reduction measures, optimizing C:P and N:P ratios while ensuring adequate nutrient availability for crops (Table 4). Moreover, it mitigates the risk of phosphorus runoff and leaching, simultaneously improving soil fertility and protecting the environment.

Ratios involving K (C:K, N:K, P:K) were generally lower than crop uptake requirements, suggesting limitations of C, N, and P relative to high background K (Li et al., 2022). Partial organic substitution, especially 30%M, improved these ratios more effectively than NPK, though imbalances persisted at some sites (e.g., declines at YX). Given that K is less incorporated into microbial biomass yet critical for crop uptake (Wang et al., 2024), maintaining balanced stoichiometry through combined organic-inorganic fertilization is essential for long-term sustainability in Vertisols (Chen et al., 2020).

Crop yield and soil quality responded strongly to fertilization, with partial organic substitution outperforming chemical fertilizer alone. Improvements in the Soil Quality Index (SQI) under 15%M and 30%M corresponded to significant yield gains, particularly in maize, reflecting the close link between soil quality and productivity in this system (Figure 5, Table 5). This is consistent with (Engels, 1993) who reported that maize, as a C4 crop, undergoes a pronounced “explosive growth phase” from jointing to filling, during which nitrogen demand accounts for over 70% of total growth. The gradual nitrogen release from organic fertilizer mineralization aligns with this peak demand. In contrast, wheat, as an overwintering crop, experiences markedly reduced mineralization rates during its reviving stage due to low temperatures, creating a mismatch between nutrient supply and demand; consequently, wheat relies more heavily on chemical fertilizers to meet its nitrogen requirements during this critical period (Page et al., 1978).

Random forest analysis identified SQI as the dominant predictor of yield (26–31%), followed by available N, emphasizing the importance of soil quality. Structural equation modeling indicated that yield gains were mediated through improvements in stoichiometry, nutrient availability, bulk density, and SOC, with nutrient pathways exerting the strongest influence (Figures 6, 7). These results align with previous studies showing that integrated nutrient management enhances both soil function and crop performance more effectively than sole chemical fertilization (Chen et al., 2020; Singh et al., 2024).

Importantly, the observed site-specific variation—for example, stronger benefits for bulk density at GY and MC, or superior yield performance of 15%M at KT and MC—highlights the need to tailor substitution rates to local soil constraints. Overall, a moderate organic substitution rate (~30%) appears effective in improving soil structure, nutrient availability, and crop yield in Vertisols. Nevertheless, site-specific conditions such as compaction and organic matter content should guide local adjustment of substitution rates. Future practices should adopt tailored organic–inorganic strategies to sustain soil fertility and crop productivity.