Soybean-driven microbiome remodeling of soil metabolism enhances nano-ZnO-assisted remediation of TPH-contaminated soil

Sustainable strategies are needed to remediate total petroleum hydrocarbon (TPH)-contaminated soils. A 90-day pot experiment combined soybean cultivation with different concentrations of nano-zinc oxide (nano-ZnO), followed by soil physicochemical, amplicon sequencing, microbial network, metabolomic, and multi-omics analyses. Combined application of soybean and 60 mg/kg nano-zinc oxide (nano-ZnO) significantly enhanced soil fertility and enzyme activities, achieving a 68.7% total petroleum hydrocarbon (TPH) removal rate over 90 days. Amplicon sequencing revealed soybean cultivation primarily structured bacterial and fungal communities. The rhizosphere buffered microbial diversity against TPH and nanoparticle stress, and co-occurrence networks showed the combined treatment fostered the most complex, stable cross-kingdom interactions. Soil metabolomics indicated a targeted rhizosphere response to nano-ZnO, enriching flavonoids and biosurfactant precursors, contrasting with the non-specific stress response in bare soil. Multi-omics integration identified a functional guild (Pseudomonas, Sphingomonas, Rhodococcus, and Fusarium) associated with TPH dissipation. These results suggest that soybean remodels its rhizosphere microbiome structure and metabolic function to synergize with nano-ZnO, accelerating the recovery of petroleum-contaminated soils.

Soil is the fundamental substrate for terrestrial ecosystems and the cornerstone of global food security; however, it is increasingly threatened by anthropogenic pollution caused by intensifying industrial and agricultural activities (1). Among the myriad of pollutants, total petroleum hydrocarbons (TPH) present a particularly complex remediation challenge due to their hydrophobicity, recalcitrance, and tendency to persist in soil matrices, necessitating advanced sustainable remediation strategies (2, 3). Traditional physicochemical remediation approaches, such as soil washing and thermal desorption, are often economically unviable and ecologically disruptive, sparking a global shift towards eco-friendly biological solutions (3). Phytoremediation, the use of plants and their associated microbiomes to detoxify soils contaminated with petroleum hydrocarbons, has emerged as a viable green technology (3, 4). However, the efficiency of phytoremediation is frequently limited by the bioavailability of hydrophobic pollutants and the slow metabolic activity of native microbes, highlighting the need for innovative auxiliary technologies (4).

Soybean (Glycine max) is not only a globally vital crop but also serves as a powerful “ecosystem engineer” capable of improving soil health and fertility in degraded environments (5). Its extensive root system and nutrient-rich exudates shape a unique rhizosphere niche that recruits beneficial microorganisms capable of tolerating and degrading organic contaminants (5). The soybean rhizosphere microbiome plays a pivotal role in plant health by enhancing nutrient acquisition and mitigating abiotic stresses, including petroleum hydrocarbon toxicity (6, 7). Specifically, plant growth-promoting rhizobacteria (PGPR) within this niche can fix nitrogen, solubilize phosphate, and produce phytohormones, which are essential for sustaining plant growth in TPH-contaminated soils (8, 9). Recent reviews have emphasized that harnessing these microbial partners is key to the commercialization of biostimulants and the advancement of sustainable agriculture (10, 11).

Nanotechnology has recently intersected with agriculture to offer novel solutions for stress management and soil remediation (12, 13). Nano-zinc oxide (nano-ZnO) is of particular interest; at optimal concentrations, it acts as an essential micronutrient source that can regulate plant primary metabolites and antioxidant activities, thereby improving plant stress responses under contaminated-soil conditions (14, 15). Studies have shown that nano-ZnO can reduce contaminant uptake and improve plant physiological status under stress conditions (16, 17). However, the interaction between nanomaterials and the soil microbiome is complex and concentration-dependent. While nanoparticles can stimulate specific microbial defense mechanisms at low doses (18), they also pose potential antimicrobial risks at higher concentrations that must be carefully managed to avoid disrupting ecosystem functions (14, 18).

A critical gap remains in understanding the synergistic interplay between plants, nanomaterials, and specific functional microbial guilds in the context of petroleum-contaminated soil restoration (3, 19). Successful bioremediation of TPH-contaminated soils relies on microbial taxa with versatile catabolic capabilities. Among bacteria, Pseudomonas species are renowned for their diverse catabolic plasmids, production of biosurfactants, and ability to initiate oxidative ring-cleavage of aromatic hydrocarbons via dioxygenase pathways (20). Sphingomonas species, characterized by their unique glycosphingolipid-containing outer membrane, exhibit tolerance to hydrophobic organic contaminants and are frequently enriched in petroleum-contaminated soils (2, 20). Rhodococcus species are notable for their broad catabolic versatility, possession of multiple alkane hydroxylase gene copies, and capacity to form biofilms that enhance contact with hydrophobic substrates (21). In the fungal domain, Fusarium species may contribute to TPH dissipation through the “fungal highway” mechanism, whereby fungal hyphae facilitate bacterial dispersal and provide structural conduits for pollutant transport across otherwise disconnected soil pores (22, 23). Together, these taxa may form a degradation-associated consortium that supports the dissipation of complex petroleum hydrocarbon mixtures (20, 22).

Therefore, this study tests the hypothesis that soybean cultivation can actively “remodel” the soil microbiome to unlock a synergistic interaction with nano-ZnO, converting the nanoparticle from a potential stressor into a functional input for the remediation of TPH-contaminated soil. By integrating soil metabolomics with microbial network analysis, we aim to decipher how this tripartite system (Soybean–Microbe–Nano-ZnO) enhances soil fertility (24) and accelerates pollutant dissipation. We specifically focus on identifying the key functional guild, including Pseudomonas, Sphingomonas, Rhodococcus, and Fusarium, may contribute to this remediation process that may contribute to TPH dissipation (25, 26), and on elucidating the metabolic reprogramming through which the soybean rhizosphere orchestrates this response.

Soil preparation and experimental setup

The soil used in this experiment was collected from the surface layer (0–20 cm) of an agricultural field located in the suburbs of Chengdu, Sichuan Province, China. The soil is pedologically classified as typical purple soil (Eutric Regosol in FAO taxonomy), which is widely distributed in the Sichuan Basin and known for its rapid weathering and high mineral content (27). After removing plant residues and stones, the soil was air-dried naturally, ground, and passed through a 2-mm nylon sieve. To simulate petroleum hydrocarbon contamination, crude oil obtained from a local refinery was dissolved in acetone and thoroughly mixed with the prepared soil to achieve a final total petroleum hydrocarbon (TPH) concentration of 10,000 mg kg -1. The spiked soil was then aged in the dark at room temperature for 60 days to allow for the volatilization of toxic light components and the establishment of an adsorption-desorption equilibrium, with periodic mixing to ensure homogeneity.

The experiment followed a randomized complete block design with a full factorial arrangement of two factors: soybean cultivation (bare soil vs. soybean planted) and nano-ZnO application at four levels (0, 30, 60, and 90 mg kg -1). This resulted in eight treatment groups: CK, ZnO30, ZnO60, ZnO90, CKSoy, ZnO30Soy, ZnO60Soy, and ZnO90Soy. The 60 mg kg-1 dose was selected for detailed multi-omics interpretation because it produced the highest TPH removal in the dose-response experiment, whereas the 90 mg kg-1 treatment showed a slight decline in removal efficiency. Commercially available spherical nano-ZnO particles ([brand/supplier to be added]; purity > 99.5%, average particle size < 100 nm) were suspended in deionized water and ultrasonicated for 30 min to ensure uniform dispersion before being thoroughly mixed with the aged petroleum-contaminated soil. Each plastic pot (20 cm diameter, 18 cm height) was filled with 5 kg of the treated soil. For soybean treatments, seeds of Glycine max (L.) Merr. cultivar “Nandou 12,” a variety widely cultivated in Southwest China for its stress tolerance, were surface-sterilized with 70% ethanol for 1 min and 3% H 2 O 2 for 5 min, followed by rinsing five times with sterile distilled water. Four seeds were sown in each pot and thinned to two uniform seedlings after one week of emergence. The pots were maintained in a greenhouse under controlled conditions with a 14/10 h light/dark photoperiod, a temperature regime of 28/22 °C (day/night), and a relative humidity of approximately 70%. Soil moisture was maintained at 60% of the field water holding capacity by weighing and watering daily with deionized water. There were six biological replicates for each treatment. After a 90-day growth period, rhizosphere soil samples were carefully collected by gently shaking the roots and collecting the tightly adhering soil, immediately frozen in liquid nitrogen, and stored at −80 °C for molecular and metabolomic analyses, while a portion was air-dried for physicochemical characterization.

The experiment followed a randomized complete block design with a full factorial arrangement of two factors: soybean cultivation (bare soil vs. soybean planted) and nano-ZnO application at four levels (0, 30, 60, and 90 mg kg -1). This resulted in eight treatment groups: CK, ZnO30, ZnO60, ZnO90, CKSoy, ZnO30Soy, ZnO60Soy, and ZnO90Soy. Commercially available spherical nano-ZnO particles (purity > 99.5%, average particle size < 100 nm) were suspended in deionized water and ultrasonicated for 30 min to ensure uniform dispersion before being thoroughly mixed with the aged petroleum-contaminated soil. Each plastic pot (20 cm diameter, 18 cm height) was filled with 5 kg of the treated soil. For soybean treatments, seeds of Glycine max (L.) Merr. cultivar “Nandou 12,” a variety widely cultivated in Southwest China for its stress tolerance, were surface-sterilized with 70% ethanol for 1 min and 3% H 2 O 2 for 5 min, followed by rinsing five times with sterile distilled water. Four seeds were sown in each pot and thinned to two uniform seedlings after one week of emergence. The pots were maintained in a greenhouse under controlled conditions with a 14/10 h light/dark photoperiod, a temperature regime of 28/22 °C (day/night), and a relative humidity of approximately 70%. Soil moisture was maintained at 60% of the field water holding capacity by weighing and watering daily with deionized water. There were six biological replicates for each treatment. After a 90-day growth period, rhizosphere soil samples were carefully collected by gently shaking the roots and collecting the tightly adhering soil, immediately frozen in liquid nitrogen, and stored at −80 °C for molecular and metabolomic analyses, while a portion was air-dried for physicochemical characterization.

Analysis of soil physicochemical properties and TPH degradation

Soil pH was measured in a soil-water suspension (1:2.5, w/v) using a digital pH meter. Soil total carbon (TC) and total nitrogen (TN) were quantified using an elemental analyzer. Total potassium (TK) was determined by flame photometry following digestion with aqua regia. Total phosphorus (TP) was measured using the molybdenum blue colorimetric method after perchloric acid digestion. Available phosphorus (AP) was extracted with 0.5 M NaHCO 3 and measured using the molybdenum blue colorimetric method. Soil nitrate-nitrogen (NO 3 - -N) and ammonium-nitrogen (NH 4 + -N) were extracted with 2 M KCl and analyzed using a continuous flow analyzer.

To evaluate the remediation efficiency, the residual concentration of total petroleum hydrocarbons (TPH) in the soil was determined by ultrasonic extraction followed by gas chromatography, according to the standard method USEPA 8015D (28). Briefly, 2 g of freeze-dried soil was mixed with anhydrous sodium sulfate and extracted with 20 mL of dichloromethane/hexane (1:1, v/v) in an ultrasonic bath for 30 min. The extract was filtered, concentrated using a rotary evaporator, and purified through a silica gel column. The TPH concentration was quantified using a gas chromatograph equipped with a flame ionization detector (GC-FID).

Soil enzyme activities, which serve as sensitive bioindicators of soil health and metabolic capability, were assayed using standard colorimetric methods. Urease activity was determined by the phenol-sodium hypochlorite method after incubating soil with urea at 37 °C for 24 h and expressed as μg NH 3 -N g -1 (dw) h -1 (29). Catalase (CAT) activity was assessed by potassium permanganate titration based on the decomposition rate of H 2 O 2 and expressed as mg O 2 g -1 (dw) s -1 (29). Amylase activity was determined using starch as a substrate, with the reducing sugars released quantified by the 3,5-dinitrosalicylic acid (DNS) colorimetric method and expressed as mg C 6 H 12 O 6 mg -1 (dw) (29). Acid phosphatase (ACP) activity was measured using p-nitrophenyl phosphate as a substrate, with the released p-nitrophenol quantified colorimetrically at 405 nm and expressed as μg PNP g -1 (dw) h -1 (30).

Soil DNA extraction, PCR amplification, and amplicon sequencing

Total genomic DNA was extracted from 0.5 g of fresh rhizosphere soil using the DNeasy PowerSoil Pro Kit (QIAGEN, Germany) according to the manufacturer’s detailed instructions. The quality and concentration of extracted DNA were verified using a NanoDrop 2000 spectrophotometer and 1% agarose gel electrophoresis. For bacterial community analysis, the V3–V4 hypervariable region of the 16S rRNA gene was amplified using the primers 341F and 806R (31). For fungal community analysis, the ITS1 region was amplified using primers ITS1F and ITS2R (32). All PCR reactions were carried out using Phusion ® High-Fidelity PCR Master Mix. PCR products were purified using the QIAquick Gel Extraction Kit and quantified using a Qubit 2.0 Fluorometer. Sequencing libraries were generated using the TruSeq ® DNA PCR-Free Sample Preparation Kit and sequenced on an Illumina NovaSeq 6000 platform to generate 250 bp paired-end reads.

Bioinformatics analysis

Raw sequencing data were processed using the QIIME 2 pipeline (version 2023.2) (33). Primers were trimmed, and paired-end reads were merged, denoised, and chimera-filtered using the DADA2 plugin to generate Amplicon Sequence Variants (ASVs) (34). Taxonomy was assigned to ASVs using the SILVA database (release 138) for bacteria (35) and the UNITE database (release 8.2) for fungi (36). To normalize sequencing depth, samples were rarefied to the lowest sequence count across all samples prior to diversity analysis. Alpha diversity indices (Shannon and Chao1) were calculated to assess community richness and diversity. Beta diversity was evaluated based on Bray-Curtis distances and visualized using Principal Coordinate Analysis (PCoA). Permutational multivariate analysis of variance (PERMANOVA) was performed using the vegan package (37) in R to test for significant differences in community structures among treatments.

Co-occurrence network analysis

To decipher microbial interactions under petroleum stress and nano-ZnO treatments, co-occurrence networks were constructed for four representative treatment groups: CK, CKSoy, ZnO60, and ZnO60Soy. To reduce noise, only ASVs present in at least 80% of samples within each group were retained. Pairwise Spearman’s rank correlations were calculated between all retained ASV pairs, and a robust correlation was defined as having a Spearman’s correlation coefficient |ρ| > 0.6 and a statistically significant p-value < 0.01, adjusted for false discovery rate using the Benjamini-Hochberg procedure (38). The resulting networks were visualized using Gephi software (version 0.9.2) (39), with nodes colored by phylum-level taxonomy. Topological properties, including average degree, modularity, average path length, and graph density, were calculated to characterize network complexity and stability. Keystone taxa were identified based on their high degree and betweenness centrality scores within the networks (40).

Soil metabolomics analysis

Soil metabolites were extracted from 100 mg of freeze-dried soil using 1 mL of ice-cold methanol/water (3:1, v/v) containing L-2-chlorophenylalanine as an internal standard. Metabolite profiling was performed using a UHPLC system coupled with a Q Exactive High-Resolution Mass Spectrometer (Thermo Fisher Scientific), following the protocols described by Dunn et al. (41). Chromatographic separation was carried out on an ACQUITY UPLC HSS T3 column using water and acetonitrile containing 0.1% formic acid as mobile phases. The mass spectrometer was operated in both positive and negative electrospray ionization (ESI) modes. Raw data were processed using Compound Discoverer 3.1 for peak alignment, normalization, and identification against the KEGG and mzCloud databases. Differential metabolites were identified based on a Variable Importance in Projection (VIP) score > 1.0 from a Partial Least Squares-Discriminant Analysis (PLS-DA) model and a p-value < 0.05 from a Student’s t-test. The metabolomics comparisons were conducted between: (i) ZnO60 vs. CK (bulk soil response to nano-ZnO), and (ii) ZnO60Soy vs. CKSoy (rhizosphere response to nano-ZnO), to delineate the divergent metabolic trajectories between planted and unplanted systems.

Statistical and multi-omics correlation analysis

Statistical analyses were performed using SPSS 26.0 (IBM Corp., Armonk, NY, USA) and R software (v.4.2.0). One-way analysis of variance (ANOVA) followed by Tukey’s HSD post-hoc test was used to compare soil physicochemical properties, TPH degradation rates, and alpha diversity indices among treatments. The normality of data was checked using the Shapiro-Wilk test prior to ANOVA. To link microbial community structure with soil function, Spearman’s rank correlation analyses were conducted between the relative abundance of dominant microbial genera and key soil environmental variables, including TPH concentration, enzyme activities, and fertility indicators (TC, TN, AP). Mantel tests were performed to assess the overall relationship between microbial community composition and TPH degradation rates. Correlations between the relative abundance of key microbial OTUs and the intensity of differential metabolites were calculated and visualized as heatmaps using the pheatmap package in R, separately for bacterial (Figure 1) and fungal (Figure 2) communities.

Synergistic enhancement of TPH degradation and soil enzymatic activities by soybean and nano-ZnO

The remediation efficiency of petroleum-contaminated purple soil was significantly influenced by the interaction between soybean cultivation and nano-ZnO application. As shown in Figure 3A, the residual concentration of TPH decreased significantly across all treatments compared to the initial spiked concentration of 10,000 mg kg -1, but the extent of degradation varied markedly among experimental groups. The bare soil control (CK) exhibited limited natural attenuation, with a TPH removal rate of only 18.5% over the 90-day period, reflecting the constraints of indigenous microbial activity under petroleum stress. In contrast, soybean cultivation alone (CKSoy) significantly elevated the degradation rate to 42.3%, confirming the “rhizosphere effect” whereby root exudates stimulate microbial metabolism and enhance pollutant bioavailability. Most notably, the co-application of nano-ZnO further accelerated this process in a concentration-dependent manner. The moderate concentration treatment (ZnO60Soy) achieved the highest TPH removal rate of 68.7%, significantly outperforming both the bare soil and the soybean-only treatments, while the highest nano-ZnO dose (ZnO90Soy) yielded a slightly lower removal rate, suggesting that excessive nano-ZnO may exert inhibitory effects on the rhizosphere microbiome. These results indicate that at an optimal concentration, nano-ZnO acts as a stimulant rather than a toxicant, enhancing the phytoremediation potential of soybean. This dose-response pattern provided the basis for focusing the subsequent network and metabolomic interpretation on the ZnO60Soy treatment.

Concurrently, soil enzymatic activities mirrored the TPH degradation trends and served as sensitive bioindicators of the functional status of the soil ecosystem (Figures 3B–E). Urease activity, which reflects the nitrogen cycling capacity of the soil microbial community, was significantly elevated in the combined soybean and nano-ZnO treatments compared to bare soil controls, with the ZnO60Soy and ZnO90Soy groups showing the highest values (Figure 3B). Catalase (CAT) activity, an indicator of the soil’s capacity to detoxify reactive oxygen species generated during hydrocarbon degradation, followed a similar trend, with the ZnO90Soy treatment achieving the highest activity, approximately 2.2-fold that of the CK control (Figure 3C). Amylase activity, reflecting the carbon cycling potential of the microbial community, was significantly enhanced in the soybean-planted treatments, with values in the ZnO60Soy and ZnO90Soy groups significantly exceeding all other treatments (Figure 3D). Acid phosphatase (ACP) activity, a key indicator of phosphorus mobilization, showed the most pronounced increase under combined treatment, with the ZnO90Soy group achieving the highest activity of approximately 105 μg PNP g -1 (dw) h -1, more than twice that of the CK control (Figure 3E). Furthermore, soil fertility metrics, including total carbon, total nitrogen, available phosphorus, and inorganic nitrogen species, were significantly improved in the combined treatments (Table 1), likely reflecting the enhanced mineralization of organic matter facilitated by the stimulated microbial community.

Soybean cultivation dominates the reshaping of soil microbial community structure

To elucidate the biological correlates of the enhanced TPH degradation, high-throughput amplicon sequencing was employed to characterize both bacterial and fungal communities across all treatments. Principal Coordinate Analysis (PCoA) based on Bray-Curtis dissimilarities revealed distinct clustering patterns driven by the experimental factors (Figures 4A, B). For bacterial communities, the first principal coordinate (PCo1, explaining 37.02% of the total variation) clearly separated samples based on the presence or absence of soybean, while nano-ZnO treatments formed sub-clusters along PCo2 (22.27%) within each main group (Figure 4A). A similar pattern was observed for fungal communities, with PCo1 (27.2%) separating planted from unplanted treatments and PCo2 (15.94%) capturing the secondary effects of nano-ZnO (Figure 4B). These results indicate that the rhizosphere effect was the dominant force structuring the soil microbiome, overriding the effects of nano-ZnO concentration and petroleum stress. PERMANOVA confirmed that soybean cultivation accounted for a significantly greater proportion of community variation than nano-ZnO treatment for both bacteria and fungi (p < 0.001).

Analysis of alpha diversity provided further insights into community responses to the experimental treatments (Figures 4C, D). For bacterial communities, the Shannon diversity index was highest in the CK and low nano-ZnO bare soil treatments and declined progressively with increasing nano-ZnO concentration in unplanted soils, with the ZnO90 treatment showing the lowest diversity (Figure 4C). This reduction reflects the known antimicrobial effects of metal oxide nanoparticles at high concentrations. Strikingly, this negative effect was significantly mitigated in the soybean-planted treatments; the Shannon index of the ZnO90Soy group was significantly higher than that of ZnO90, demonstrating that the rhizosphere acts as a “microbial refuge” buffering the community against the combined stress of petroleum hydrocarbons and nanoparticles. For fungal communities, the ZnO30Soy treatment exhibited the lowest Shannon diversity, while the other soybean-planted groups maintained relatively stable diversity indices that were not significantly different from each other (Figure 4D), suggesting a more complex and non-linear fungal response to the combined treatments.

Enhanced network complexity and stability in the rhizosphere under moderate nano-ZnO exposure

Beyond simple diversity metrics, the ecological interactions among microbial taxa were deciphered through co-occurrence network analysis of four representative treatment groups (CK, CKSoy, ZnO60, ZnO60Soy). The topological properties of the bacterial networks revealed that soybean cultivation significantly increased network complexity compared to bulk soil (Figures 5A, B). The CKSoy network featured a greater number of nodes and edges, as well as higher average degree and graph density relative to CK, indicating that the rhizosphere environment fosters a more interconnected and potentially cooperative bacterial community structure. The nano-ZnO-only treatment (ZnO60) produced a distinct network architecture characterized by the emergence of clearly separated modules dominated by Proteobacteria and Actinomycetota, suggesting a community-level reorganization in response to nanoparticle stress (Figure 5C). Strikingly, the highest network complexity among all treatments was observed in the combined treatment (ZnO60Soy), where the bacterial network displayed the greatest number of inter-module connections and a higher proportion of positive correlations (Figure 5D). This suggests that under moderate nano-ZnO pressure, the rhizosphere microbiome responded by strengthening cooperative interactions, such as cross-feeding and syntrophy, to efficiently process the complex hydrocarbon mixture.

Analogous patterns were observed in the fungal co-occurrence networks (Figures 5E–H). The fungal network of the CK treatment was dominated by Ascomycota with relatively dense intra-phylum connections (Figure 5E). Soybean cultivation (CKSoy) introduced greater inter-phylum connectivity, with Basidiomycota and Mucoromycota nodes forming new linkages with the dominant Ascomycota cluster (Figure 5F). The ZnO60 treatment resulted in a sparser fungal network with fewer edges compared to the control (Figure 5G), consistent with the known sensitivity of certain fungal taxa to metal oxide nanoparticles. In the ZnO60Soy treatment, however, the fungal network recovered substantially in terms of connectivity, with an increased proportion of Ascomycota nodes exhibiting high betweenness centrality, indicative of their role as keystone taxa bridging otherwise disconnected fungal modules (Figure 5H). Taken together, these network analyses suggest that the combined application of soybean and moderate nano-ZnO promotes the formation of specialized and stable functional consortia across both bacterial and fungal kingdoms, consistent with the “Intermediate Disturbance Hypothesis”.

Divergent metabolic trajectories in rhizosphere and bulk soils in response to nano-ZnO

Untargeted metabolomics analysis was conducted to link the structural shifts in the microbiome to functional metabolic outputs. The metabolomics data are presented across three analytical layers, each addressing a distinct comparison (Figure 6). The first row (Figures 6A–C) compares the overall metabolic profiles between ZnO60 vs. CK and ZnO60Soy vs. CKSoy, revealing that nano-ZnO induced far greater metabolic perturbation in bare soil (552 differential metabolites in total) than in the soybean rhizosphere (67 differential metabolites), with only 2 metabolites shared between the two systems (Figure 6C). The second row (Figures 6D–F) and third row (Figures 6G–I) present analogous comparisons for distinct metabolite subsets, consistently showing that the bare soil system exhibited a broader, less specific response (83 and 43 unique differential metabolites, respectively) compared to the more targeted rhizosphere response (35 and 28 unique differential metabolites, respectively).

In bare soil, the metabolic response to nano-ZnO was characterized by a broad upregulation of stress-related compounds, consistent with a generalized cellular stress response to nanoparticle toxicity (Figure 6A). In stark contrast, the metabolic shift in the rhizosphere was substantially more targeted and specific (Figure 6B), with pathway enrichment analysis indicating selective enrichment of flavonoid-related compounds and specific amino acids implicated in biosurfactant synthesis. These soybean-derived flavonoids are well-documented chemotactic signals for the recruitment of beneficial rhizosphere microorganisms and are known to enhance the bioavailability of hydrophobic organic contaminants such as TPH. The remarkably low overlap between the bare soil and rhizosphere metabolic responses (only 2 shared differential metabolites across all comparisons; Figures 6C, F, I) confirms that the soybean rhizosphere microbiome did not simply resist nano-ZnO stress but actively reprogrammed its metabolic output in response to the nanoparticle signal, a fundamentally different ecological strategy from the non-specific stress response observed in bare soil.

Identification of keystone taxa linking microbiome structure to remediation function

To identify microbial taxa associated with the observed TPH remediation, a multi-omics correlation analysis was performed integrating OTU relative abundance data with differential soil metabolite profiles. The resulting correlation heatmaps are presented separately for bacterial OTUs (Figure 1) and fungal OTUs (Figure 2). Hierarchical clustering of the bacterial heatmap (Figure 1) revealed two major OTU clusters with contrasting correlation patterns. The first cluster, comprising OTUs affiliated with Pseudomonas and Sphingomonas, exhibited strong positive correlations (r > 0.8, p < 0.01) with a suite of metabolites associated with aromatic hydrocarbon catabolism and biosurfactant precursor synthesis, and strong negative correlations with residual TPH concentration. The second cluster included OTUs most negatively correlated with the pro-degradation metabolite suite, likely representing taxa suppressed under the combined soybean and nano-ZnO treatment. Within the bacterial network, OTUs assigned to Rhodococcus occupied central positions with high betweenness centrality, consistent with their role as metabolic hubs connecting aliphatic and aromatic hydrocarbon degradation pathways.

In the fungal correlation heatmap (Figure 2), OTUs affiliated with Fusarium formed a distinct cluster characterized by strong positive correlations with organic acid metabolites and amino acid biosurfactant precursors, and positive co-occurrence with the Pseudomonas - and Sphingomonas -affiliated bacterial OTUs identified above. This cross-kingdom co-occurrence pattern is consistent with the “fungal highway” mechanism, whereby fungal hyphae serve as physical conduits facilitating bacterial dispersal and the transport of hydrophobic substrates across discontinuous soil pore networks. A subset of Ascomycota OTUs, exhibiting high degree centrality in the ZnO60Soy fungal network (Figure 5H), also showed significant positive correlations with flavonoid-class metabolites in Figure 2, suggesting that these fungi may respond to plant-derived signals to co-assemble with the bacterial degraders. Collectively, the multi-omics correlation analysis provides correlation-based evidence that the synergistic application of soybean and nano-ZnO reshapes the soil microbiome to favor a cross-kingdom functional guild—comprising Pseudomonas, Sphingomonas, Rhodococcus, and Fusarium —that may contribute to TPH dissipation through complementary catabolic and structural roles.

The remediation of petroleum-contaminated soil represents a complex ecological challenge where the interactions between plants, microbes, and soil amendments collectively determine the restoration trajectory. Our findings suggest that the co-application of soybean and moderate concentrations of nano-ZnO (60 mg kg -1) yields a significant synergistic effect on TPH degradation, achieving a removal rate of 68.7% over 90 days and substantially surpassing the efficacy of either treatment applied in isolation. This synergy is underpinned by the concurrent stimulation of multiple soil enzyme activities, including urease, catalase, amylase, and acid phosphatase, which collectively reflect an enhanced capacity for nutrient cycling, oxidative stress detoxification, and organic matter mineralization in the rhizosphere (42, 43). The rhizosphere effect observed here aligns with the “priming effect” theory, whereby root exudates provide labile carbon and nitrogen sources that stimulate the microbial production of extracellular enzymes, thereby co-metabolizing recalcitrant petroleum hydrocarbons (44). The significant elevation of soil fertility indicators—particularly total carbon, available phosphorus, and inorganic nitrogen—in the combined treatments (Table 1) further supports the interpretation that the soybean rhizosphere actively drives nutrient remobilization as a by-product of enhanced hydrocarbon mineralization.

Crucially, our data suggest that nano-ZnO, at an optimal dose of 60 mg kg -1, may function as a micronutrient-related cofactor rather than solely as a stressor within the soybean rhizosphere. Zinc is an essential structural and catalytic component of numerous metalloenzymes, including superoxide dismutase, carbonic anhydrase, and dehydrogenases involved in hydrocarbon catabolism. The elevated catalase and urease activities observed in the ZnO60Soy treatment imply that the controlled bioavailability of Zn² + derived from the partial dissolution of nanoparticles may have relieved a micronutrient limitation in the rhizosphere, thereby enhancing the metabolic efficiency of the hydrocarbon-degrading microbial guild (45). The decline in remediation efficiency at the highest nano-ZnO dose (ZnO90Soy), despite its continued stimulatory effect on certain enzyme activities, suggests that the relationship between nano-ZnO concentration and ecosystem function follows a hormetic dose-response curve, where benefits are realized only within a narrow optimal window (46). This observation challenges the prevailing view of engineered nanomaterials solely as environmental contaminants, suggesting their potential utility as precision “nano-fertilizers” in bioremediation strategies when deployed within this hormetic window.

The divergent responses of the soil microbiome to nano-ZnO in the presence and absence of soybean highlight the critical protective role of the rhizosphere environment. In bulk soil, high concentrations of nano-ZnO induced a significant decline in bacterial Shannon diversity (Figure 4C), consistent with the known antimicrobial mechanisms of metal oxide nanoparticles, including membrane disruption, generation of reactive oxygen species, and interference with cellular zinc homeostasis (14, 18). However, the soybean rhizosphere successfully buffered these negative impacts, maintaining relatively high bacterial diversity even under the combined stress of petroleum contamination and nano-ZnO exposure. This “microbial refuge” phenomenon can be attributed to the physicochemical complexity of the rhizosphere environment: soybean root exudates, rich in organic acids such as citrate and malate, likely chelated free Zn² + ions, reducing their bioavailability to toxic levels while simultaneously providing abundant carbon substrates to sustain diverse microbial populations (47). The more complex and non-linear response of fungal communities to the combined treatments (Figure 4D) suggests that fungi occupy distinct ecological niches in the rhizosphere and respond to nano-ZnO through different physiological mechanisms than bacteria, a distinction that warrants further investigation.

The co-occurrence network analysis revealed that the combined treatment (ZnO60Soy) produced the most complex and interconnected microbial networks for both bacteria (Figure 5D) and fungi (Figure 5H), consistent with the “Intermediate Disturbance Hypothesis,” which posits that moderate environmental stress prevents competitive exclusion and promotes functional diversity (48). The increased proportion of positive correlations and higher inter-module connectivity in the ZnO60Soy bacterial network (Figure 5D) suggest the emergence of cooperative functional consortia, where different microbial taxa engage in syntrophic interactions to collectively degrade the complex mixture of aliphatic and aromatic hydrocarbons present in petroleum-contaminated soil (23, 38). The reorganization of the ZnO60 bacterial network into distinct Proteobacteria- and Actinomycetota-dominated modules (Figure 5C) is consistent with the known differential sensitivity of these phyla to metal oxide nanoparticles and their contrasting roles in aromatic versus aliphatic hydrocarbon degradation, respectively. The recovery of fungal network connectivity in the ZnO60Soy treatment relative to ZnO60 alone (Figures 5G vs. 4H) is consistent with the possibility that soybean shields the fungal community from nanoparticle toxicity, preserving the structural integrity of the cross-kingdom microbial network that may be relevant to TPH dissipation.

The metabolomics analysis provided a molecular window into the active recruitment strategies employed by soybean under combined petroleum and nano-ZnO stress and yielded an integrative set of findings of this study. The strikingly low overlap between the metabolic responses of bare soil and rhizosphere soil to nano-ZnO—with only 2 shared differential metabolites across all three analytical comparisons (Figures 6C, F, I)—indicates that the soybean-microbiome holobiont did not simply buffer nanoparticle stress but fundamentally reprogrammed its metabolic output in response to the nano-ZnO signal. In bare soil, the metabolic response was characterized by a broad upregulation of non-specific stress compounds (Figure 6A), consistent with a community-wide emergency response to acute nanoparticle toxicity. In contrast, the rhizosphere metabolic response was highly targeted, featuring selective enrichment of flavonoid-class compounds and amino acid-derived biosurfactant precursors (Figure 6B). Flavonoids such as daidzein and genistein are well-documented chemotactic signals in the soybean rhizosphere that recruit nitrogen-fixing rhizobia and other beneficial PGPR (49). In the context of petroleum contamination, these isoflavonoids likely serve an additional function as structuring signals that selectively assemble a TPH-degrading microbiome and simultaneously act as antioxidants to mitigate ROS generated by hydrocarbon toxicity (50). The detection of biosurfactant precursor metabolites in the ZnO60Soy treatment further indicates a targeted metabolic investment in enhancing TPH bioavailability, a rate-limiting step in the biodegradation of hydrophobic petroleum components.

The multi-omics correlation analysis (Figures 1, 2) nominated a cross-kingdom functional guild associated with TPH remediation in the combined treatment. Among bacteria, Pseudomonas and Sphingomonas OTUs exhibited the strongest positive correlations with the pro-degradation metabolite suite and the strongest negative correlations with residual TPH concentration (Figure 1). Pseudomonas species are among the most metabolically versatile hydrocarbon degraders, possessing diverse catabolic plasmids encoding dioxygenase and monooxygenase enzymes capable of initiating the oxidative ring-cleavage of aromatic hydrocarbons, as well as the capacity to produce rhamnolipid biosurfactants that enhance TPH bioavailability (20). Sphingomonas species, characterized by their unique glycosphingolipid-containing cell envelopes, are exceptionally tolerant of hydrophobic organic contaminants and have been frequently identified as keystone taxa in petroleum-contaminated soils (2, 20). Rhodococcus OTUs, identified as metabolic hubs in the bacterial co-occurrence network, are notable for their broad catabolic versatility and possession of multiple alkane hydroxylase gene copies, positioning them as critical linkers between aliphatic and aromatic degradation pathways (21). In the fungal domain, Fusarium OTUs showed strong positive correlations with organic acid metabolites and co-occurred consistently with the bacterial degraders (Figure 2), supporting the “fungal highway” mechanism whereby fungal hyphae serve as physical conduits for bacterial dispersal and substrate transport across disconnected soil pore networks (22, 23). The strong correlation between plant-derived flavonoid metabolites and the abundance of Fusarium -affiliated OTUs suggests that soybean actively recruits fungal partners through chemical signaling, extending the plant’s microbiome engineering capacity to the fungal kingdom. Taken together, these findings support a working model in which soybean uses targeted metabolic signaling to assemble a cooperative bacterial-fungal consortium that collectively converts nano-ZnO from a potential toxicant into a functional micronutrient input, thereby co-associated with TPH dissipation.

Despite these promising findings, several limitations must be acknowledged. First, the experiment was conducted in a controlled pot system using artificially spiked and aged soil, which, while ensuring experimental homogeneity, does not fully replicate the spatial heterogeneity, weathering history, and indigenous microbial diversity of actual petroleum spill sites. The behavior of nano-ZnO and its interactions with the soil matrix may differ substantially under fluctuating field conditions involving rainfall leaching, temperature extremes, and competing organic matter inputs (14, 51). Second, the study employed a single particle size and morphology of nano-ZnO; given that the reactivity and toxicity of nanomaterials are shape- and size-dependent, future studies should systematically evaluate different aspect ratios and surface coatings to optimize remediation outcomes (51). Third, while the present study identified key functional taxa and metabolic pathways through correlation-based approaches, the specific genetic mechanisms underpinning the observed degradation activity—such as the expression levels of alkB, nahH, or rhdA genes encoding alkane hydroxylases and ring-hydroxylating dioxygenases—were not directly quantified. Integration of metatranscriptomic or qPCR-based functional gene profiling in future work would provide definitive molecular evidence for the proposed degradation pathways (21). Finally, the long-term ecological risks of nano-ZnO and petroleum-residue transfer into the food chain remain unassessed in the present study. We did not measure available Zn, plant biomass, Zn accumulation in soybean tissues, residual TPH or compound-specific PAHs in plant tissues, or longer-term post-remediation soil-health indicators beyond the measured physicochemical and enzyme endpoints. A longitudinal field trial incorporating nano-ZnO fate tracking, plant Zn analysis, TPH/PAH residue assessment in soybean tissues, and crop safety evaluation is essential to verify that the observed remediation gains do not incur unacceptable ecological or food safety costs (13, 52).

Despite these promising findings, several limitations must be acknowledged. First, the experiment was conducted in a controlled pot system using artificially spiked and aged soil, which, while ensuring experimental homogeneity, does not fully replicate the spatial heterogeneity, weathering history, and indigenous microbial diversity of actual petroleum spill sites. The behavior of nano-ZnO and its interactions with the soil matrix may differ substantially under fluctuating field conditions involving rainfall leaching, temperature extremes, and competing organic matter inputs (14, 51). Second, the study employed a single particle size and morphology of nano-ZnO; given that the reactivity and toxicity of nanomaterials are shape- and size-dependent, future studies should systematically evaluate different aspect ratios and surface coatings to optimize remediation outcomes (51). Third, while the present study identified key functional taxa and metabolic pathways through correlation-based approaches, the specific genetic mechanisms underpinning the observed degradation activity—such as the expression levels of alkB, nahH, or rhdA genes encoding alkane hydroxylases and ring-hydroxylating dioxygenases—were not directly quantified. Integration of metatranscriptomic or qPCR-based functional gene profiling in future work would provide definitive molecular evidence for the proposed degradation pathways (21). Finally, the long-term ecological risks of nano-ZnO accumulation in the food chain, particularly its potential trophic transfer into soybean seeds and subsequent implications for food safety, remain unassessed in the present study. A longitudinal field trial incorporating nano-ZnO fate tracking and crop safety evaluation is essential. A longitudinal field trial incorporating nano-ZnO fate tracking, plant Zn analysis, TPH/PAH residue assessment in soybean tissues, and crop safety evaluation is essential to verify that the observed remediation gains do not incur unacceptable ecological or food safety costs (13, 52).

This study demonstrates a novel synergistic mechanism through which soybean cultivation and moderate nano-zinc oxide application at 60 mg per kilogram co-drive the restoration of petroleum-contaminated soil. This study supports a synergistic working model in which soybean cultivation and moderate nano-zinc oxide application at 60 mg per kilogram enhance the remediation of petroleum-contaminated soil. Soybean functions as an active ecosystem engineer, shaping the structure and metabolic function of the rhizosphere microbiome through targeted exudate-mediated signaling. This plant-driven microbiome remodeling converts the potential toxicity of nano-zinc oxide into a functional micronutrient input, promoting the assembly of a complex, cooperative, and highly efficient cross-kingdom degradation consortium candidate cross-kingdom degradation-associated consortium comprising Pseudomonas, Sphingomonas, Rhodococcus, and Fusarium. The concurrent enhancement of soil fertility indicators and enzyme activities confirms the concurrent enhancement of soil fertility indicators and enzyme activities suggests that pollutant removal and ecosystem recovery proceed as coupled processes within this tripartite system. Soil metabolomics further establishes. Soil metabolomics further indicates that the soybean rhizosphere does not passively tolerate nanoparticle stress but actively reprograms its metabolic output to favor targeted hydrocarbon catabolism over generalized stress responses. This plant-assisted nano-bioremediation strategy offers a promising approach for the restoration of petroleum-contaminated agricultural soils, provided that the long-term fate of nano-zinc oxide and potential TPH or PAH transfer into soybean tissues are carefully evaluated before food-chain application.