The main tobacco variety CB-1 in the tobacco-growing area of Fujian province was used as experimental material. The pot experiment was conducted from April to June 2024 at the Fujian Key Laboratory of Crop Breeding by Design, situated within the greenhouse of Fujian Agriculture and Forestry University. The greenhouse environment was maintained at a 16-h day temperature of 22 °C and an 8-h night temperature of 18 °C. The experimental design comprised a control group (denoted as CK) without the addition of either C₁₀H₅Cl₂NO₂ or K₂S₂O₈ to the soil. The treatment groups included soil added with 0.04 mg/kg C₁₀H₅Cl₂NO₂ (denoted as C), and soil amended with both 0.04 mg/kg C₁₀H₅Cl₂NO₂ and 100 mg/kg K₂S₂O₈ (denoted as CY). Each treatment involved five plants, with three biological replicates for a total of 45 pots.

The samples were collected at 45 days post-treatment. This included soil samples for soil characteristics, rhizospheric soil for microbial community analysis, and tobacco leaves for gene expression and metabolic profiling of the tobacco plants. For the microbial community analysis, soil samples within a range of 1–4 mm around the roots of the three treatment groups were collected, with approximately 100 g per treatment group. Additionally, 100 g of soil was collected from five pots per treatment group for soil characteristics analysis. For gene expression and metabolic profiling of tobacco, leaves in the 2nd–3rd positions from the top (counting from the uppermost leaf) were selected for this study. A total of ten leaves, sourced from five different plants, were collected as a sample. All samples were collected with three biological replicates. Following collection, the soil samples were stored at −20 °C, while the leaf samples were stored at −80 °C until analysis.

The collected soil samples were initially purified to remove impurities, and then passed through a 2 mm mesh for homogenization. The air-dried soil samples were subsequently analyzed for their properties and nutrient content. The organic matter was analyzed according to the NY/T1121.6-2006 method; pH was determined using the NY/T1377-2007 method; total nitrogen level was measured following the NY/T53-1987 method; total phosphorus content was assessed based on the NY/T88-1988 method; total potassium concentration was determined according to the NY/T87-1988 method; available nitrogen were evaluated using the method described in LY/T 1228-2015; available phosphorus levels were evaluated using the NY/T1121.7-2006 method; and available potassium concentration was measured following the NY/T889-2004 method. Soil particle size measurements were carried out according to the NY/T1121.3-2006 method.

The high-throughput metagenomic sequencing technology was used to investigate changes in soil microbial communities under three treatments. Initially, microbial DNA was extracted and purified from soil samples using a bacteria and fungi genomic extraction kit (Omega D3350-02; Solarbio D2300-100T). The resulting DNA fragments were generated through ultrasound treatment, followed by purification, end-repair, 3′-end adenylation, and ligation with sequencing adapters. Subsequently, agarose gel electrophoresis was employed to select appropriately-sized fragments for PCR amplification library construction. Metagenome sequencing was performed on the Illumina Hiseq2500 platform following standard protocols. After data processing and statistical analysis, including low-quality data filtering, output data generation, and quality control statistics, the metagenome assembly was carried out using MEGAHIT software while QUAST software (Tang and Borodovsky, 2014) evaluated the assembly results by removing contig sequences shorter than 300 bp. Additionally, MetaGeneMark software was used for coding region identification and removal of redundant data. Finally, prediction analysis of tobacco rhizosphere microbial community structure and alpha diversity under different treatments was conducted on BMK Cloud.¹

RNA-Seq was used for the treatments and their control to investigate the potential mechanism underlying K₂S₂O₈-mediated C₁₀H₅Cl₂NO₂ stress mitigation. Total RNA of samples (C, Y, CK) was extracted using the TRIzol reagent (Invitrogen, USA). RNA sequencing (RNA-Seq) and data processing were performed with the Illumina HiSeq platform at Biomarker Technologies Co., LTD. (Beijing, China) according to Cho et al. (2016). The RNA-Seq data have been submitted in the NCBI Sequence Read Archive (SRA) under the accession number PRJNA1221589.

After excluding reads containing adapter, poly-N, and low-quality sequences, the remaining clean reads were aligned to the reference genome in Sol Genomics Network database.² Subsequently, these aligned reads were assembled and quantitatively analyzed using StringTie software to determine the fragments per kilobase of exon per million fragments mapped (FPKM) values. DEGs were identified using a false discovery rate (FDR) ≤ 0.01 and |Fold change| ≥ 1.5 while calculating FDR and Fold change (FC) for all genes. Additionally, GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed for the three comparison groups: CK vs. C, CK vs. CY, and C vs. CY.

The metabolites of leaf samples at 45 days post-treatment (C, CY, and CK) were analyzed using widely targeted metabolomics methods. Freeze-dried leaves were homogenized using a mixer mill (MM 400, Retsch, Germany), and the leaf powder was pooled from each biological replicate sample, 100 mg of this powder was extracted overnight at 4 °C with 0.6 mL of 70% aqueous methanol. The extracts were subjected to analysis using ultra-performance liquid chromatography with electrospray ionization coupled to tandem mass spectrometry (UPLC–ESI–MS/MS) at Biomarker Technologies Co., LTD. (Beijing, China).

Metabolic data from each sample were analyzed using hierarchical cluster analysis (HCA), principal component analysis (PCA), and K-means clustering. HCA and PCA analyses were performed using Software R and GraphPad Prism v9.01 (GraphPad Software Inc., La Jolla, CA, USA), respectively. DEMs among samples from different groups were identified based on the following criteria: VIP ≥ 1, |Fold change| ≥ 1, and p value <0.01. The Venn diagram illustrates the quantitative relationship among different comparison groups. The Kyoto Encyclopedia of Genes and Genomes (KEGG) compound database³ was utilized for annotating the different metabolites which were then mapped onto the KEGG pathway database.⁴ Pathways containing significantly regulated metabolites underwent further analysis through metabolite sets enrichment analysis (MSEA). Significance assessment was conducted by calculating p-values obtained from hypergeometric tests.

The Spearman test method (Heinen and Valdesogo, 2020) was employed to conduct correlation analysis among metabolomics, transcriptomics, and microbiota. Results meeting the criteria of a p-value <0.05 and a Spearman correlation coefficient |r| > 0.8 were chosen for constructing a correlation network.

Total RNA was isolated from plantlets using TRIzol reagent (Invitrogen) according to the manufacturer’s protocol. The extracted RNA was then reverse-transcribed into complementary DNA (cDNA), which was used for quantitative real-time PCR (qRT-PCR) analysis with SYBR Premix ExTaq (Takara). The expression of the Actin gene was employed as an internal control. The experiment was conducted with three biological replicates, each comprising three individual plants, and each sample was analyzed in triplicate. The relative gene expression levels were determined using the 2^−ΔΔCt method (Livak and Schmittgen, 2001), and the primer sequences used for qRT-PCR are provided in Supplementary Table 1.

In comparison to the control (CK) (Figure 1A), tobacco seedlings exposed to C₁₀H₅Cl₂NO₂ herbicides exhibited leaf curling/narrowing and stunted growth (Figure 1B). Notably, K₂S₂O₈ treatment significantly alleviated these symptoms (Figure 1C), suggesting its potential role in mitigating herbicide-induced damage.

The physicochemical analysis of the soil (Figure 2) showed C₁₀H₅Cl₂NO₂ stress significantly reduced available nitrogen, phosphorus, and potassium (C vs. CK; Figures 2F–H). Importantly, K₂S₂O₈ application (CY) counteracted these reductions, restoring available nitrogen and potassium to near-CK levels (Figures 2F,H). Soil particle analysis (Figures 2I–J) further revealed that C₁₀H₅Cl₂NO₂ altered granular structure (>0.01 mm vs. <0.01 mm), while K₂S₂O₈ rehabilitated proportions to CK-equivalent states. These results demonstrate K₂S₂O₈‘s dual capacity to alleviate herbicide damage and restore soil functionality.

Metagenomic sequencing of root-associated communities yielded 364,770,378 clean reads from nine samples (treatments C, CY, and the control CK; Supplementary Table 2). Assembly generated 1,133,605 contigs (N50 > 680 bp), with open reading frame (ORF) prediction identifying 2,265,851 ORFs, confirming dataset robustness for further analysis.

Alpha diversity analysis using Shannon, Simpson, and Inverse-Simpson indices demonstrated that C₁₀H₅Cl₂NO₂ (C) exposure significantly reduced microbial diversity relative to CK. In contrast, K₂S₂O₈ (CY) amendment not only reversed this decline but enhanced diversity beyond control levels (Figures 3A–C), indicating effective mitigation of herbicide impacts on soil microbiota.

Taxonomic profiling (Supplementary Table 3) identified 4 kingdoms, 185 phyla, 305 classes, 495 orders, 928 families, 2,818 genera, and 12,641 species. Bacteria dominated microbial communities (96.05–97.99% relative abundance), with archaea constituting the remainder (Figure 3D). Five core bacterial phyla-Proteobacteria, Acidobacteria, Actinobacteria, Chloroflexi, and Gemmatimonadetes-collectively represented 83.42–87.29% of relative abundance across treatments (Figure 3E). Notably, C₁₀H₅Cl₂NO₂ reduced the abundance of beneficial genera including Acidobacteria, Chloroflexi, Gemmatimonadetes, Bradyrhizobium, Actinobacteria, Verrucomicrobia, Candidatus-Rokubacteria, and Sphingomonas, while K₂S₂O₈ treatment uniquely restored their prevalence (Figure 3F). These genus-specific shifts substantiate K₂S₂O₈’s capacity to rehabilitate functional soil microbiomes compromised by quinclorac stress.

To elucidate the mechanism by which K₂S₂O₈ alleviates C₁₀H₅Cl₂NO₂ stress, we integrated transcriptomic and metabolomic analyses of tobacco leaves under CK, C, and CY treatments. PCA distinguished the C group from CK and CY (Figure 4A). Notably, C₁₀H₅Cl₂NO₂ stress (C vs. CK) induced 3,019 down-regulated and 2,146 up-regulated genes, while K₂S₂O₈ supplementation (CY vs., C) reversed this trend, up-regulating 851 genes and down-regulating 627 genes (Figure 4B). Crucially, 71 DEGs were common across all comparisons (CK vs. C, CK vs. CY, C vs. CY; Figure 4C). Go enrichment confirmed that DEGs were primarily associated with stress response pathways, including cellular process (GO:0009987), metabolic process (GO:0008152), biological regulation (GO:0065007), localization (GO:0051179), response to stimulus (GO:0050896), signaling (GO:0023052). Critically, C₁₀H₅Cl₂NO₂ suppressed expression in these pathways (down-regulated > up-regulated in CK vs. C), while K₂S₂O₈ restored expression levels, directly supporting its role in mitigating phytotoxicity (Figure 4D).

Metabolite profiling revealed distinct clustering among CK, C, and CY groups (Figures 5A–C), validating data robustness. We identified 1,396 metabolites, dominated by terpenoids (17.9%), lipids (11.6%), organic acids (9.03%), sugars/alcohols (8.88%), and amino acids (8.88%) (Figure 5C).

OPLS-DA confirmed significant inter-group differences (Figures 6A–C), and volcano plots quantified metabolites change: CK vs. C had 62 increased and 71 decreased metabolites, while C vs. CY showed 36 increased and 26 decreased metabolites, indicating K₂S₂O₈’s normalization effect (Figures 6D–F). Among 203 differentially expressed metabolites (DEMs), seven were shared across all comparisons (Figure 6G). K-means clustering demonstrated that C₁₀H₅Cl₂NO₂ stress specifically depleted metabolites in cluster 3 and cluster 4 but elevated those in cluster 2 and cluster 5. Remarkably, K₂S₂O₈ supplementation (CY) restored these metabolites to near-control (CK) levels (Figure 6H), highlighting its efficacy in rescuing stress-disrupted metabolic pathways.

To reveal the core metabolic pathways modulated by K₂S₂O₈ in alleviating C₁₀H₅Cl₂NO₂ stress, we conducted an integrative analysis of metagenome and transcriptome. KEGG analysis identified DEGs and DEMs co-enriched in six key metabolic pathways, including lysine degradation, stilbenoid diarylheptanoid and gingerol biosynthesis, arginine and proline metabolism, phenylalanine biosynthesis, tyrosine metabolism, and flavonoid biosynthesis (Figures 7A–C). A total of 159 DEGs were identified in the six common metabolic pathways, which were classified into 6 clusters by K-means analysis based on their similar expression patterns. Critically, C₁₀H₅Cl₂NO₂ stress (C vs. CK) significantly suppressed gene expression in cluster 2 and cluster 4, while inducing expression in cluster 6. Strikingly, K₂S₂O₈ supplementation (CY) effectively reversed these stress-induced alterations, restoring expression levels in these clusters close to those observed in the control (CK) (Figure 7D). This restoration pattern strongly supports K₂S₂O₈’s role in counteracting C₁₀H₅Cl₂NO₂-induced dysregulation within these critical pathways. Additionally, compared to control (CK), genes in cluster 1 showed decreasing expression in both C and CY treatments, while genes in cluster 3 exhibited increasing expression in both treatments, though CY induced more pronounced changes than C alone. Among the seven DEMs identified in these six common metabolic pathways, most displayed significantly altered levels under C₁₀H₅Cl₂NO₂ stress (CK vs. C) but were restored towards control levels by K₂S₂O₈ supplementation (CK vs. C, CK vs. CY) (Figure 7E).

A multi-omics correlation network highlights K₂S₂O₈-mediated regulatory interactions. To elucidate the interplay among metabolomics, transcriptomics, and microbiota, we constructed a correlation network comprising 7 common DEMs, 159 DEGs, and microbial taxa with relative abundance greater than 1% (Figure 7F). The network consists of 79 nodes and 173 edges (96 positive, 77 negative). This included gene-microbe (50 nodes, 112 edges), gene-metabolite (18 nodes, 22 edges), and microbe-metabolite (2 nodes, 3 edges) interactions. Notably, a significant negative correlation was identified between the metabolite 5-Aminovaleric Acid and the microbe Sphingomonas regulated by seven genes (NewGene_3900, Nitab4.5_0000215g0020, Nitab4.5_0000791g0070, Nitab4.5_0000006g0050, Nitab4.5_0001220g0050, Nitab4.5_0000107g0090, and NewGene_21840). Furthermore, 4-Hydroxyphenyacetylgultamic Acid exhibited a significant negative correlation with Bradyrhizobium and Alphaproteobacteria regulated by Nitab_4.50002015 g00700 and Nitab450006992g00700, respectively. These specific regulatory axes underscore the complex interplay between the microbiome, gene expression, and metabolite levels potentially modulated by K₂S₂O₈ in mitigating stress.

To verify the reliability of the transcriptome data, we selected six genes that showed significant correlations with both the metabolome and microbiome for qPCR validation. Compared with the CK group, the C and CY groups showed consistent trends, with four genes up-regulated and two genes down-regulated. Notably, the gene expression profile of the CY group was more closely aligned with that of the CK group compared to the C group. This observation serves as further evidence of the efficacy of the K₂S₂O₈ treatment. The expression patterns of these six genes obtained through qRT-PCR were highly consistent with those from RNA-seq (Figure 8), confirming the reliability of the transcriptome-based differential gene expression analysis.