The study was carried out at Xiyi Township, Longyang District, Baoshan City, Yunnan Province (99°30′E, 24°93′N), at an altitude of 1,602 m, which was a typical mountainous tobacco growing area with a subtropical monsoon climate, an average annual temperature of 16.20 °C and an average annual precipitation of 827.0 mm. The type of test soil was red soil with clay texture (17.44% sand, 26.90% silt and 55.66% clay). The basic physicochemical properties of the test soil were shown in Table 1.

In this study, three tillage treatments with different tillage depths were implemented during the fallow season prior to tobacco planting from 2019 to 2021: rotary tillage (RT) with a tillage depth of 15 cm, deep tillage with middle depth (DTM) with a depth of 25 cm, and deep plowing (DP) with a depth of 35 cm (Figure 1). The experiment followed a completely randomized design with three replicates per treatment, totaling 9 plots. Each plot measured 144 m² (20 m × 7.2 m). Each year, treatments were applied once before planting using machinery specific to the designated depth: a hand-held rotary tiller (6.3 kW, model 1WG6.3) for RT, a crawler-type automatic tiller (14.7 kW, model DH1115) for DTM, and a deep tiller (73.5 kW, model LR4A3L-23) for DP. Following these primary tillage operations, all plots were uniformly tilled one final time with the hand-held rotary tiller to ensure a consistent seedbed condition across all treatments. Agronomic practices were consistent for all plots. Prior to transplanting, base fertilizers were applied in strips, including 450 kg·hm⁻² of compound fertilizer (N 7%, P₂O₅ 16%, K₂O 26%), 300 kg·hm⁻² of nitrogen-potassium fertilizer (N 16%, K₂O 30%), and 1,200 kg·hm⁻² of commercial organic fertilizer (N 2.77%, P₂O₅ 1.48%, K₂O 1.15%). Tobacco seedlings were transplanted in early May each year, using a row spacing of 1.20 m and a plant spacing of 0.50 m. Each plot contained 240 plants, resulting in a planting density of approximately 16,667 plants hm⁻².

Soil samples were collected in August 2021 (during the tobacco harvest period) using a 5 cm diameter auger at two soil depths: 0–20 cm and 20–40 cm. Soil samples were taken from 10 sites arranged in an S-shape within each plot and composited into a single replicate sample after removing roots and debris. The soil samples were divided into three parts: one portion was stored at 4 °C for analysis of microbial biomass carbon (MBC), microbial biomass nitrogen (MBN), and enzyme activities; another portion (0–20 cm only) was stored at −80 °C for microbial DNA extraction and sequencing; and the remaining portion was air-dried and sieved through a 0.25 mm mesh for determination of physicochemical properties.

Soil organic carbon (SOC) was determined using the potassium dichromate oxidation method (Bao, 2000). Total nitrogen (TN) was measured by the Kjeldahl digestion method (Bao, 2000). Ammonium nitrogen (NH₄⁺-N) and nitrate nitrogen (NO₃⁻-N) were extracted with 1 mol/L KCl and measured using a continuous flow analyzer (AA3, SEAL Analytical, Germany). Available phosphorus (AP) was determined by the molybdenum-antimony anti-colorimetric method, and available potassium (AK) was analyzed using a flame photometer (Bao, 2000). Soil pH was measured in a 1:2.5 soil-to-water suspension using a pH meter. Bulk density (BD) and soil porosity (SP) were determined by the core method (Xue et al., 2018). Soil sucrase and cellulase activities were determined by the 3,5-dinitrosalicylic acid colorimetric method (Guan et al., 1986). Urease activity was measured by the indophenol blue colorimetric method (Guan et al., 1986). MBC and MBN were determined using the chloroform fumigation-extraction method. Fresh soil samples (10 g) were fumigated with ethanol-free CHCl₃ for 24 h at 25 °C in a vacuum desiccator. Both fumigated and non-fumigated samples were then extracted with 0.5 mol L⁻¹ K₂SO₄ (soil-to-solution ratio of 1:4) by shaking for 30 min. The extracts were filtered and analyzed using a TOC/TN analyzer (Multi N/C 3100, Analytik Jena, Germany). MBC and MBN were calculated based on the difference in extractable C and N between fumigated and non-fumigated samples, using conversion factors of k_EC = 0.45 and k_EN = 0.54, respectively.

Soil genomic DNA was extracted using the FastDNA® Spin Kit for Soil (MP Biomedicals, United States). The amoA genes were amplified using specific primers: Arch-amoAF (5-STAATGGTCTGGCTTAGACG-3) and Arch-amoAR (5-GCGGCCATCCATCTGTATGT-3) for AOA (Francis et al., 2005), and amoA-1F (5-GGGGTTTCTACTGGTGGT-3) and amoA-2R (5-CCCCTCKGSAAAGCCTTCTTC-3) for AOB (Rotthauwe et al., 1997). Detailed procedures for the PCR amplification process are provided in Supplementary methods. The purified amplicons were used to construct paired-end (PE) 2 × 250 bp libraries following the standard protocol for the Illumina MiSeq platform (Illumina, San Diego, United States). Sequencing was conducted on the Illumina MiSeq PE250 platform by Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). The bioinformatics analysis was performed using the Majorbio Cloud Platform.¹ The raw sequencing reads were demultiplexed and quality-filtered using fastp (version 0.19.6) and merged using FLASH (version 1.2.7). Operational Taxonomic Units (OTUs) were clustered at a 97% similarity cutoff using UPARSE (version 7.1), and chimeric sequences were identified and removed. The taxonomy of each representative OTU sequence was analyzed using the RDP Classifier (Wang et al., 2007) against the FunGene functional gene database (Fish et al., 2013) with a confidence threshold of 0.7. To eliminate the effect of sequencing depth on diversity measures, the dataset was rarefied to the minimum number of sequences per sample prior to downstream analysis.

Tobacco yield was determined by harvesting, flue-curing, and weighing the leaves from 50 plants in the center of each plot. After curing, leaves were graded into superior, medium, and inferior categories to calculate the proportion of each grade. The production value was then calculated by multiplying the weight of each grade by the corresponding market price for that year (Dong et al., 2025; Hu W. et al., 2021).

We used two-way ANOVA to assess the effects of tillage, soil depth, and their interaction on soil physicochemical properties, enzyme activities, and microbial biomass. One-way ANOVA was used to evaluate the effect of tillage on microbial alpha diversity indices. Significant differences among treatments were identified using Tukey’s honestly significant difference (HSD) test at a significance level of p < 0.05. Alpha diversity indices (Shannon and Sobs) of AOA and AOB were calculated, and their taxonomic composition was visualized by relative abundance. We used Mantel tests to correlate community structure with environmental variables and Pearson correlation analysis to explore relationships between microbial diversity and soil properties. All of these analyses were conducted using the “vegan” package (Oksanen et al., 2020). To infer microbial interactions, we constructed co-occurrence networks for AOA and AOB using the “igraph” package (Csárdi and Nepusz, 2006). The networks were built from Spearman correlation matrices, retaining only strong (|r| > 0.6) and significant (p < 0.05) correlations. Key topological parameters (e.g., average degree, modularity) were calculated to characterize network complexity. Finally, we performed a Partial Least Squares Path Modeling (PLS-PM) analysis to investigate the complex causal relationships between soil properties, microbial features, and tobacco production benefits. The model included six latent variables: (1) Soil physicochemical property (pH, SP, BD); (2) Soil nutrients content (TN, AK, AP, SOC, NH₄⁺-N); (3) Microbial community (Shannon diversity of AOA and AOB); (4) Soil microbial biomass (MBC, MBN); (5) Soil enzyme (sucrase, cellulase, urease activity); and (6) Production benefit (tobacco yield and production value). The PLS-PM was conducted using the “plspm” R package. All statistical analyses were performed using R software (version 4.4.1).

Two-way ANOVA showed that tillage treatments significantly affected SOC, TN, NO₃⁻-N, AP, and AK, whereas NH₄⁺-N was not significantly affected. Soil depth had significant effects on NO₃⁻-N, AP, and AK, and significant interactions between tillage and depth were observed for AP and AK. At the 0–20 cm depth, RT significantly increased SOC, TN, NH₄⁺-N, AP, and AK compared to DTM and DP (Figures 2A–C,E,F). However, NO₃⁻-N was highest under DP (Figure 2D), indicating enhanced nitrification. In the 20–40 cm layer, RT still maintained the highest SOC, TN, and AP levels (Figures 2A,B,E), while DP had the highest NO₃⁻-N and AK (Figures 2D,F).

Tillage also significantly influenced BD and SP, both of which also varied significantly with soil depth. The interaction effect was not significant for either parameter. Compared to DTM and DP, RT significantly increased BD in both layers, reaching 1.15 g cm⁻³ at 0–20 cm and 1.17 g cm⁻³ at 20–40 cm (Figure 2G). In contrast, DTM and DP showed higher SP than RT at both depths (Figure 2H). Soil pH values were stable across treatments, ranging from 7.99–8.19 (0–20 cm) and 8.06–8.10 (20–40 cm). No significant effects of tillage, depth, or interaction were observed on pH (Figure 2I).

Tillage treatments had significant effects on urease activity, but not on cellulase or sucrase. Urease also varied significantly with depth, and all three enzymes showed significant interaction effects between tillage and depth. At 0–20 cm, urease activity showed no significant differences among treatments (Figure 2J). However, DTM significantly reduced cellulase and sucrase activities compared to RT and DP (Figures 2K,L). At 20–40 cm, urease and cellulase activities were significantly higher under RT than DP (p < 0.05), while sucrase activity showed no significant differences among treatments (Figures 2J–L).

Tillage treatments significantly affected both MBC and MBN, and soil depth had significant effects on both. A significant interaction was observed only for MBC. Across both soil layers, RT consistently produced the highest MBC values (303.82 mg kg⁻¹ at 0–20 cm and 443.07 mg kg⁻¹ at 20–40 cm), while DTM showed the lowest (Figure 2M). For MBN, RT also showed the highest values at both depths (Figure 2N). At 0–20 cm, MBN under DTM was the lowest, while at 20–40 cm, MBN followed the trend: RT > DTM > DP (Figure 2N).

Norank_d_Bacteria, norank_p_Thaumarchacota and norank_p_Crenarchaeota were the top three taxa with the highest relative abundance among the AOA (Figure 3A). Specifically, norank_d_Bacteria and norank_p_Crenarchaeota were enriched in RT, while norank_p_Thaumarchacota was enriched in DTM and DP. Norank_d_Bacteria, Nitrosospira and Nitrosomonadaceae were the top three taxa with the highest relative abundance among AOB (Figure 3B). Specifically, norank_d_Bacteria were enriched in DTM, while Nitrosomonadaceae and Nitrosospira were enriched in DP. Mantel analysis explained the effects of environmental factors on the community composition of AOA and AOB (Figure 3C). The community structure of AOA was influenced by NO₃⁻-N, TN, AP, AK, BD, and SP. While NH₄⁺-N and Sucrase mainly regulated the community structure of AOB.

The ANOVA test for diversity revealed that the sobs indices of AOA (p = 0.12) and AOB (p = 0.14) communities were insensitive to the response of tillage treatments (Figures 4B,D). RT and DP significantly increased the Shannon index of the AOA community compared to DTM (p < 0.001; Figure 4A). For the AOB community, DP showed a significantly higher Shannon index than both RT and DTM (p < 0.001; Figure 4C). Pearson’s correlation analysis of AOA and AOB community diversity with soil properties revealed that Shannon’s index of AOB was significantly associated with soil nutrient content and soil physico-chemical properties (Figure 4E). Specifically, the Shannon index of AOB was significantly negatively correlated with SOC, TN, NH₄⁺-N, AP, AK, and BD while it was significantly positively correlated with NO₃⁻-N and SP. The Shannon index of the AOA community was significantly and positively correlated with NO₃⁻-N, soil enzyme (Sucrase, Cellulase) and MBN. The sobs index of the AOB community was positively correlated with NO₃⁻-N, whereas the sobs index of the AOA community was significantly negatively correlated with SOC, TN, NH₄⁺-N.

Microbial co-occurrence network analysis reveals changing patterns of interactions between AOA and AOB community under different tillage treatments (Figure 5). For the co-occurring network of AOA, the tillage treatments had a limited effect on the number of nodes incorporated into the network (20–22) but improved the number of links by 54.76 and 38.06% for DTM and DP, respectively, compared to RT (Figure 5A). This implies that DTM and DP improved the interaction between AOA. However, the highest number of links in RT’s network was in the AOB community. Compared to DTM, RT and DP improved 61.11 and 18.06% of links, respectively (Figure 5B). Further analysis of the co-occurrence network properties of AOA and AOB revealed distinct patterns across tillage treatments. In the AOA network, DTM exhibited the highest average degree and edge density, but the lowest relative modularity, whereas DP had the highest relative modularity (Supplementary Figure S1A). In the AOB network, RT showed the highest average degree, centralization, and modularity, while DP again exhibited the highest relative modularity (Supplementary Figure S1B). This implied that DTM and RT improved the interaction of AOA and AOB communities, respectively, but DP improved the cohesion of modules performing potential functions in both colonies.

Tillage treatments had no statistically significant effect on the overall tobacco yield or production value, although the highest numerical values for both were observed under the deep plowing (DP) treatment (2,825 kg hm⁻² and 88,959 CNY hm⁻², respectively) (Table 2). However, tillage significantly impacted the quality distribution of tobacco leaves. The proportion of superior-grade leaves was significantly higher under DP (77.15%) compared to rotary tillage (RT) (68.82%). Conversely, there were no significant differences among treatments for the proportions of medium-grade or inferior-grade leaves. The PLS-PM analysis indicated that soil physicochemical properties, nutrient content, and microbial characteristics collectively explained 81.7% of the variance in production benefits (R² = 0.817) (Figure 6). Specifically, physicochemical properties, including soil pH, porosity, and bulk density, affected production benefits indirectly via soil nutrients and soil microbial biomass. Soil nutrients and the diversity of AOA and AOB communities exerted a significant positive direct effect on soil microbial biomass. In turn, soil microbial biomass not only directly and significantly promoted soil enzyme activity (standardized path coefficient = 1.325, p < 0.05) but also demonstrated the strongest positive direct influence on production benefits (standardized path coefficient = 1.303, p < 0.10).

Our results showed that RT significantly increased SOC, TN, AP, and AK in the 0–20 cm soil layer compared to DTM and DP. The increase in SOC is primarily attributed to the reduced soil disturbance under RT, which minimizes soil aeration and slows the mineralization of organic matter, thereby facilitating the retention of crop residues in the surface layer (Lv et al., 2023; Shen et al., 2024). In contrast, the accumulation of TN, AP, and AK is largely driven by the shallow incorporation of basal fertilizers. Since rotary tillage restricts soil mixing to the top 15 cm, the applied nitrogen, phosphorus, and potassium remain concentrated in the surface soil. In the 20–40 cm layer, RT still maintained higher SOC, TN, and AP levels compared to DTM and DP, which contrasts with studies suggesting that deep tillage enhanced subsoil fertility by facilitating nutrient redistribution (Feng et al., 2020; Li J. et al., 2024; Qiao et al., 2025). This inconsistency may be explained by differences in experimental duration, residue input, and rooting depth of crops. While long-term deep tillage combined with deep-rooted crops and high biomass return generally improves subsoil nutrient availability, our relatively short (3-year) tobacco field experiment, characterized by low residue input and shallow root systems, likely limited nutrient translocation to deeper layers. Regarding potassium, RT increased AK in the 0–20 cm layer, whereas DP significantly elevated AK in the 20–40 cm layer, reflecting the enhanced mixing and downward movement of potassium under deep tillage as reported previously (Mallarino and Borges, 2006).

RT significantly increased BD in both the 0–20 cm and 20–40 cm layers, while significantly reducing SP, indicating that RT resulted in a less optimal physical environment with reduced air and water movement. The resulting higher bulk density and lower porosity reduced oxygen diffusion in the soil, suppressing nitrification, as further indicated by the lower NO₃⁻-N concentrations and leading to greater accumulation of NH₄⁺-N in the 0–20 cm surface layer (Pengthamkeerati et al., 2006). In contrast, DP significantly decreased bulk density and increased porosity in both soil layers, improving soil aeration and moisture conditions (Yang et al., 2021). This enhanced aeration likely accelerated the oxidative decomposition of soil organic carbon, contributing to the significantly lower SOC levels observed under DP. This improved soil physical environment favored nitrifying microbes, promoting the conversion of NH₄⁺-N to NO₃⁻-N (Pengthamkeerati et al., 2006), which resulted in significantly higher NO₃⁻-N concentrations in both the 0–20 cm and 20–40 cm layers compared to RT and DTM.

Soil enzyme activities and microbial biomass are key indicators of biological quality and nutrient cycling (Borny et al., 2024; Neemisha and Sharma, 2022; Nikitin et al., 2022; Yu et al., 2018). In our study, RT generally enhanced these properties compared to DTM and DP. At the 0–20 cm layer, urease activity did not differ significantly among treatments, likely reflecting the buffering capacity of established microbial communities and regulation by other factors such as pH (Fisher et al., 2017). However, cellulase and sucrase activities were significantly reduced under DTM compared to RT and DP, which may be attributed to the decreased availability of organic substrates, as indicated by the lower SOC in the DTM treatment. This limitation in substrate availability likely constrained enzyme activity. In the 20–40 cm subsoil, urease and cellulase activities were significantly higher under RT than DP, suggesting that RT preserved favorable conditions for microbial activity at depth, possibly by minimizing physical disturbance and maintaining substrate availability (Shen et al., 2024). DP, on the other hand, may have disrupted aggregate structures, impairing microbial habitat quality (Basic et al., 2004; Dai et al., 2021).

RT also produced the highest MBC and MBN across both depths, consistent with findings that reduced tillage promotes microbial biomass by preserving structure and reducing disturbance (Chen et al., 2019; Dhaliwal et al., 2024; Yu et al., 2025). The lowest MBC at both 0–20 cm and 20–40 cm, along with the lowest MBN at 0–20 cm under DTM, may indicate a “dual nutrient limitation” mechanism (Demoling et al., 2007). RT maintained the highest microbial biomass primarily due to the enrichment of organic substrates (SOC) in the surface layer. In contrast, DTM diluted these organic substrates through mixing, limiting heterotrophic growth. Furthermore, unlike DP which compensated for carbon loss with significantly enhanced nitrate (NO₃⁻-N) availability, DTM failed to generate comparable inorganic nitrogen enrichment. Consequently, the soil environment under DTM lacked both the concentrated labile carbon of RT and the high nitrate levels of DP, severely restricting overall microbial proliferation (Manzoni et al., 2012).

AOA and AOB play essential roles in soil N cycling, and their community composition is highly responsive to soil management practices such as tillage (Haiming et al., 2022; Li et al., 2015; Liu et al., 2018). In our study, different tillage treatments selectively enriched dominant ammonia-oxidizing microbial groups, reflecting the environmental selection imposed by tillage practices that shape microbial community composition (Khan et al., 2023). Mantel analysis further demonstrated that AOA communities were closely associated with NO₃⁻-N, TN, AP, AK, BD, and SP, while AOB communities were mainly influenced by NH₄⁺-N and sucrase activity. This indicates a functional niche differentiation between AOA and AOB, with AOA more influenced by soil nutrient status and physical conditions, whereas AOB respond mainly to immediate NH₄⁺ availability and labile carbon sources (Carey et al., 2016; Sterngren et al., 2015).

Although the sobs index of both AOA and AOB showed no significant response to tillage, the Shannon diversity index, which accounts for both richness and evenness, was significantly affected. Specifically, RT and DP significantly enhanced AOA diversity compared to DTM, while AOB diversity was highest under DP. These results suggest that tillage intensity influences microbial diversity patterns, with RT and DP promoting greater diversity likely by maintaining or improving habitat heterogeneity and resource availability (Yu et al., 2025; Zhang et al., 2022). This reduction in diversity under DTM can also be explained by the ‘dual nutrient limitation’ mechanism described above (Demoling et al., 2007). While RT supported diverse heterotrophic populations through abundant surface carbon, and DP fostered specific nitrifying communities via high nitrate availability, DTM offered neither resource advantage. This scarcity of both organic substrates and inorganic nitrogen likely created a stressful, homogeneous environment that intensified competitive exclusion among microbial taxa, thereby inhibiting the coexistence of diverse species and resulting in the lowest community diversity. This highlights the importance of considering both the nature and intensity of soil management practices on the complex responses of microbial communities (Chernov and Semenov, 2021). Correlation analysis revealed that AOA and AOB diversity were distinctly influenced by different soil nutrients, enzyme activities, and physical properties, suggesting that they occupy complementary niches that promote coexistence and stabilize soil N cycling under varying conditions.

It is important to note that while soil physicochemical properties were assessed at two depths (0–20 cm and 20–40 cm), the microbial community analysis (including DNA extraction and sequencing of AOA and AOB) in this study was restricted to the 0–20 cm soil layer. Therefore, the observed shifts in AOA and AOB diversity and community structure directly reflect the response of the topsoil community to tillage disturbances. Although the improved soil porosity and significant nitrate accumulation observed in the 20–40 cm subsoil under DP suggest that deep plowing may also favorably alter the subsoil microbial niche, direct molecular evidence for these deep-layer changes was not obtained in this study. Future research involving deep-profile sequencing is necessary to fully elucidate the vertical spatial dynamics of nitrogen-cycling microorganisms under different tillage depths.

Microbial co-occurrence networks offer valuable insights into ecological interactions, cooperation, and competition among microbial taxa in response to environmental changes (Gao et al., 2022; Hu X. et al., 2021; Li et al., 2021). In this study, different tillage practices altered the complexity and structural characteristics of both AOA and AOB networks, reflecting their distinct responses to soil disturbance. The number of AOA nodes remained stable across treatments, indicating a relatively consistent core community structure. However, both DTM and DP enhanced potential interactions among AOA taxa, suggesting that moderate disturbance and nutrient redistribution may promote network connectivity by creating more heterogeneous niches and increasing metabolic interdependence (Cornell et al., 2023; Song et al., 2023). The DTM network exhibited high connectivity but low modularity, implying greater ecological redundancy and cooperative potential among AOA taxa, whereas the high modularity observed under DP may result from microsites with distinct physicochemical conditions, driving niche differentiation and modular assembly (Guseva et al., 2022). In contrast, RT generated the most complex and centralized AOB network, indicating that low disturbance may maintain stable substrate supply and habitat structure, thereby facilitating cooperative interactions among functional groups (Hu et al., 2020). DP also showed relatively high modularity in the AOB network, possibly linked to improved oxygen availability and enhanced nitrification activity in deeper soil layers (Hu et al., 2024). Overall, these findings indicate that tillage practices not only influence the strength of interactions within AOA and AOB communities but also shape their organizational patterns, potentially affecting the stability and functional performance of soil N cycling.

Although tillage treatments did not produce statistically significant differences in total tobacco yield, DP demonstrated a clear advantage by yielding a significantly higher proportion of superior-grade leaves. This finding aligns with previous research suggesting that intensive tillage primarily benefits crop production by improving the soil’s physical environment and ecological functions (Huang et al., 2023). Practices like DP improve the soil’s physical structure (e.g., increasing porosity, reducing bulk density), which is a critical determinant of habitat for soil microorganisms (Six et al., 2004). These physical improvements indirectly enhance production benefits by positively influencing soil nutrient availability and the microbial community. PLS-PM shows that soil microbial biomass (MBC and MBN) was identified as the strongest direct positive driver of production benefits. As a key indicator of soil health (Wardle, 1992), microbial biomass directly boosted soil enzyme activity, accelerating nutrient turnover for plant uptake (Nannipieri et al., 2017) and promoting high-quality crop growth. However, it is important to note that while RT maintained the highest total microbial biomass, it did not yield the highest proportion of superior-grade leaves. This discrepancy implies that crop quality may not be driven solely by the total magnitude of biomass (MBC/MBN), but is likely the result of the synergistic interaction between microbial activity and an optimized soil environment (Bhandari et al., 2025). Unlike RT, DP improved soil aeration and porosity, which appeared to enhance the metabolic efficiency of the microbial community and favored specific functional groups (such as AOB) essential for nitrogen transformation. Thus, the coordination between microbial biomass and favorable physical conditions under DP may be more critical for determining leaf quality than the sheer biomass quantity observed under RT. Notably, while this accelerated nutrient turnover is crucial for current crop quality, it may also imply a potential risk of soil organic matter depletion over the long term. Therefore, integrating organic amendments with deep plowing strategies is recommended to compensate for carbon consumption and ensure sustainable agricultural production. This underscores the importance of cultivating a healthy soil ecosystem for high-value, sustainable production (Doran and Zeiss, 2000).