Integrative omics analysis of plant-microbe synergies in petroleum pollution remediation

Plant materials and growth conditions

Alfalfa seeds were used as the experimental material. The seeds were soaked in dH₂O overnight and then planted into soil (soil: vermiculite = 1:1). The cultivation conditions were as follows: ambient temperature 25 °C, relative humidity 60%, 16 h light (photon flux density 200 μmol/m²/s)/ 8 h dark cycle (Yang et al., 2022).

Simulated oil pollution

After growing for about 20 days, the alfalfa seedlings (Jin Huanghou) were transplanted into water cultivation containing a simulated saturated aliphatic hydrocarbon pollution mixed solution with a concentration of 1% for stress treatment (Dodecane: Cetane: Tetracosane=1:1:1). A half-strength Hoagland nutrient solution was added to support their growth (Lou et al., 2018). The alfalfa seedlings were collected at 0 h, 6 h, and 24 h after exposure to stress. The plants were ground into powder using liquid nitrogen and stored at −80 °C. Each treatment group was replicated three times.

RNA-sequencing analysis

Sequencing libraries and RNA sequencing were performed by the Meiji Biomedical Technology (Shanghai, China) using an Illumina Novaseq 6000 platform. To ensure the accuracy of subsequent bioinformatics analysis, fastp software was used to filter the raw sequencing data and obtain clean data. The quality-controlled data was aligned to the reference genome using HISAT2 software in order to obtain the mapped reads necessary for further analysis (Reference gene source: Medicago sativa; Reference genome version: Zhongmu No. 1) (Kim, Langmead & Salzberg, 2015). Mapped reads were assembled using Cufflinks software, and the results were compared with known transcripts for functional annotation (Trapnell et al., 2010).

Differential expression analysis of three groups (three biological replicates per condition) was performed using the DESeq2 package. The differentially expressed genes were identified with the criteria set as Padjust <0.05 and |log₂FC| ≥ 1. Annotation analysis of the transcripts of DEGs was performed in the GO and KEGG databases (Klopfenstein et al., 2018).

16S rRNA gene sequence analysis

Total DNA was extracted from 5 g of each soil sample using E.Z.N.A.^® soil DNA kit (Omega Bio-tek, Norcross, GA, US). The concentration and purity of DNA were assessed using a 1% agarose gel. Using the extracted DNA as a template, PCR amplification of the V3-V4 variable region of the 16S rRNA gene was performed with the upstream primer 338F: (5′-ACTCCTACGGGAGGCAGCAG-3′) and downstream primer 806R (5′-GGACTACHVGGG TWTCTAAT-3′) carrying barcode sequences (Liu et al., 2016). The PCR products from different samples were combined at an equal concentration and purified for sequencing on the Illumina PE300 platform. Chloroplast and mitochondrial sequences were removed from all samples using UPARSE v7.1 software (Edgar, 2013). The RDP classifier was used to compare the sequences with the Silva 16S rRNA gene database for taxonomic annotation of OTUs with a 50% confidence estimate (Wang et al., 2007). Additionally, functional prediction of 16S rRNA data was performed using PICRUSt2 (version 2.2.0) software (Douglas et al., 2020). A total of three groups of samples, with three biological replicates for each group. A: soil. B: A mixture of soil and oil. D-MX: A mixture of soil and oil used to grow alfalfa.

Metabolomics analysis

Metabolomics was performed as previously described with some modifications (Li et al., 2022). A total of 100 µL of concentrated hydroponic solution was transferred into a 1.5 mL centrifuge tube. Subsequently, 400 µL of extraction solution (acetonitrile: methanol = 1:1) containing 0.02 mg/mL of the internal standard L-2-chlorophenylalanine was added. After vortex mixing for 30 s, the sample was subjected to low-temperature ultrasonic extraction at 4 °C and 40 kHz for 30 min, followed by standing at −20 °C for an additional 30 min. The samples were then centrifuged at 4 °C, 13,000 g for 15 min. The supernatant was transferred and dried with nitrogen gas. The dried residues were dissolved in 100 µL solution (acetonitrile: water = 1:1) and then subjected to ultrasound extraction for 5 min at 4 °C and 40 kHz. After centrifugation at 4 °C and 13,000 g for 10 min, the supernatant was transferred to a sample vial for analysis by UHPLC -Q Exactive HF-X. The equal volumes of metabolites in all samples are mixed to prepare quality control (QC) samples, ensuring the reproducibility of the entire analysis process.

Samples were analyzed using electrospray ionization (ESI) in both positive and negative ion modes. The ESI source parameters were set as follows: Scan range (m/z), 70–1,050; Spray voltage (kV), 3.5 (ESI+) and 3.5 (ESI−); Heater temperature (^∘C), 425; Capillary temperature (^∘C), 325; Sheath gas flow rate (arb), 50; Auxiliary gas flow rate (arb), 13; Resolution (Full MS), 60,000; Resolution (MS/MS), 7,500; Normalized collision energy (eV), 20, 40, 60.

Data analysis

The raw data were processed using Progenesis QI v3.0 software (Waters Corporation, Milford, MA, USA), which performed baseline filtering, peak recognition, retention time correction, and peak alignment. This process generated a dataset containing retention time, mass-to-charge ratio (m/z), and peak intensity information. Metabolite identification was conducted by matching MS/MS spectral data in the HMDB database (http://www.hmdb.ca/) with an m/z tolerance of less than 10 ppm. The retention time and isotope similarity scores for all identified metabolites were provided in Table S7. Metabolic pathway enrichment analysis was performed using the KEGG database (http://www.genome.jp/kegg/), classifying the metabolites based on their involvement in specific pathways or biological functions.

Determination of total petroleum hydrocarbon degradation rate

A 4% mixture of petroleum hydrocarbons and soil was prepared, and the total petroleum hydrocarbon (TPH) content in the soil was assessed at various stages of alfalfa growth using a gravimetric method (Mishra et al., 2001). For each sample, 10 g of soil was subjected to sequential extraction with 100 mL each of hexane, dichloromethane, and chloroform. The combined extracts were evaporated to dryness under a nitrogen stream, and the residual TPH content was determined gravimetrically.

Statistical analyses

Data analysis and graph generation were performed using software program GraphPad Prism V. 7.0. One-way ANOVA was used for comparison of multiple data sets, with p < 0.05 as criterion for significant difference.

The selection of the differentially expressed genes

The statistical power of this experimental design, calculated in RNASeqPower was 0.748 (three biological replicates). Principal component analysis (PCA) demonstrated distinct differences between the control and treatment groups by eliminating irrelevant data. The results indicated clear variances between these groups, and the data within each group exhibited good repeatability (Fig. 1A). Differentially expressed genes (DEGs) were selected using criteria of Padjust <0.05 and |log₂FC| ≥ 1. After 6 h of simulated oil pollution stress, a total of 783 DEGs were identified, including 480 up-regulated and 303 down-regulated genes (Fig. 1B). After 24 h of stress, 582 DEGs were detected, with 434 up-regulated and 148 down-regulated genes (Fig. 1C). Venn diagram analysis revealed that 260 DEGs were common to both time points, with 523 unique DEGs after 6 h and 322 unique DEGs after 24 h (Fig. 1D). These results indicated that as the pollution time increased, the DEGs also differed.

GO and KEGG enrichment analysis

Several DEGs were implicated in the resistance to or degradation of petroleum hydrocarbon pollutants. To elucidate the signaling pathways associated with these DEGs, GO analysis and KEGG analysis were performed. The GO analysis for the WR-6 h vs CK comparison revealed significant enrichment of DEGs in pathways related to the biosynthetic and metabolic process of salicylic acid, and the oxidation–reduction process (Fig. 2A). In the WR-24 h vs CK group, DEGs were significantly enriched in pathways involved in the synthesis of secondary metabolites, secondary metabolic processes, and glutathione metabolism (Fig. 2B).

The KEGG pathway analysis revealed distinct enrichment patterns in different time points. In the WR-6 h vs CK group, DEGs were predominantly enriched in the signal pathways of plant-pathogen interactions, glutathione metabolism, plant hormone signal transduction, and MAPKs signaling pathway (Fig. 2C). In the WR-24 h vs CK group, DEGs were mainly enriched in the signal pathways of oxidation–reduction process, carbon fixation in photosynthetic process and antenna proteins in photosynthesis (Fig. 2D). Therefore, these findings suggested that alfalfa may respond to petroleum hydrocarbon pollutants through various biosynthesis pathways, abiotic stress signaling transduction, and various oxidative-reductive metabolic pathways.

Quantitative analysis of DEGs

Transcriptome sequencing identified a total of 1,431 differentially expressed genes. Among these, nine most significantly DEGs in the antioxidant-related signaling pathway (peroxidase P7, glutathione S-transferase, transcription factor WRKY7), abiotic stress-related signaling pathway, and growth and development-related signaling pathway were selected for validation. The qRT-PCR results showed that the antioxidant-related genes (P7, ascorbate oxidase gene, glutathione S-transferase gene, REDOX2, P450, WRKY7) were all up-regulated, and the abiotic stress-related transcription factors (NAC and MYB1) also showed up-regulation (Fig. 3). At the same time, the growth development-related gene GA20ox1 was up-regulated, while the flower development-related gene FT was down-regulated (Fig. 3). The expression patterns of these nine genes are consistent with the sequencing results, thus proving the accuracy and repeatability of the transcriptome data.

Microbial diversity analysis of different soil samples

Transcriptomic analysis indicated that alfalfa exposure to petroleum pollutants activates its antioxidant enzyme system. Simultaneously, 16S rRNA sequencing was conducted to analyze the diversity and composition of the microbial community in the soil. Alpha diversity analysis employed various diversity indices to evaluate the abundance and richness of microbial communities in samples. PCA analysis indicated distinct variations between the control group (A: soil, B: a mixture of soil and oil) and the treatment group (D-MX: a mixture of soil and oil used to grow alfalfa), and the data within each group shows good repeatability (Fig. 4A). The Chao and Shannon indices revealed that the microbial abundance and richness in the A group were lower compared to those in the B and D-MX groups, with no significant differences observed between the B and D-MX groups (Figs. 4B, 4C). Venn diagram analysis identified 133 unique species present in the D-MX group (Fig. 4D). These special species may need rhizosphere exudates as carbon source or nitrogen source to maintain their survival.

Differential analysis of microbial community structure

Differential analysis of microbial community structure revealed the species present in each sample at various taxonomic levels and their relative abundance. Compared with the A and B group, family-level analysis showed that Pseudomonadaceae, Nocardiaceae, Rhizobiaceae, and Caulobacteraceae were significantly enriched in D-MX group (Fig. 5A). Furthermore, genus-level analysis showed that Pseudomonas and Rhodococcus were significantly enriched in D-MX group (Fig. 5B). The analysis of the community heatmap further indicated that Pseudomonas, Rhodococcus and Brevundimonas were the three most abundant microorganisms in the D-MX group (Fig. 5C). Besides, a comparison of the abundance of 15 species with the greatest differences revealed a significant decrease in the abundance of streptomyces, while the abundance of Pseudomonas, Rhodococcus and Brevundimonas showed a significant increase (Fig. 5D). These results suggested that the three dominant species in the D-MX group may have a role in the degradation of petroleum pollutants.

Metabolomic analysis of root exudates

Root exudates play a crucial role in supporting microbial ecosystems in the rhizosphere, and they are also essential for material cycling within this environment. To investigate whether alfalfa root exudates can degrade petroleum pollutants, a metabolomic analysis was conducted on the hydroponic solution where alfalfa was grown.

The KEGG pathway analysis results indicated that these metabolites were involved in amino acid metabolism, lipid metabolism, biosynthesis of secondary metabolites and carbohydrate metabolism (Fig. 6A). Additionally, the HMDB compound classification revealed that the metabolites were predominantly categorized as lipids and lipid-like molecules, organic acids and derivatives, organoheterocyclic compounds and organic oxygen compounds (Fig. 6B). Subsequently, alfalfa was transplanted into soils contaminated with various hydrocarbons, and the oil content was measured at different time. When pantothenic acid, malic acid, and ascorbic acid selected from metabolomics were applied to the contaminated soils where alfalfa was grown, significant improvements in degradation rates were observed. Notably, pantothenic acid application increased the degradation rate by approximately 10% compared to other treatments (Fig. 6C).