Transcriptome and metabolome response of eggplant against Ralstonia solanacearum infection

Plant material

Eggplant inbred lines “NY-1” (R genotype, highly resistant to bacterial wilt) and “KY” (S genotype, highly susceptible to bacterial wilt) were obtained from the South Subtropical Crop Research Institute Chinese Academy of Tropical Agricultural Sciences. Seeds were sown in 15 cm diameter pots. The growing material was placed in pots and composed of sterile vermiculite and clay mixed in a 3:1 volume/volume ratio. Seedlings were grown under 28 °C/25 °C day/night with a 16 h light/8 h dark photoperiod condition. After 4 weeks of culture, when seedlings were at the 4-leaf stage, the culture was incubated with R. solanacearum.

Bacterial strain and inoculation

The R. solanacearum strain GMI1000-tac-EGFP was grown overnight on 2, 3, 5-triphenyl tetrazolium chloride medium at 28 °C and suspended in sterile distilled water (Xiao et al., 2021). The suspension was adjusted to 0.12 (10⁸ colony-forming units/ml) at 600 nm. The roots of eggplants were cut at 0–1 cm from the apex and then inoculated in 50 ml suspended R. solanacearum. After inoculation with R. solanacearum, plants were grown under 30 °C/32 °C day/night with a 16 h light/8 h dark photoperiod condition. Disease was rated on a scale of 0 to 4: 0 = no symptoms, 1 = 0–25% leaves wilted, 2 = 25–50% leaves wilted, 3 = 50–75% leaves wilted, 4 = 75–100% wilted and plant dead. The disease index (DI) was calculated using the formula: DI = ((N0 × 0 + N1 × 1 + N2 × 2 + N3 × 3 + N4 × 4)/(total number of plants)). N0 to N4 were the number of plants with disease rating scale values of 0 to 4, respectively. The EGFP fluorescence of R. solanacearum was detected using LUYOR-3415.

RNA-seq

At 0, 24, and 48 hpi, the roots of 10 eggplants were collected, mixed, immediately frozen in liquid nitrogen, and stored at −80 °C. Three biological replicates were established for each treatment group. Total RNA was extracted using the Spin Column Plant total RNA Purification Kit (Sangon Biotech, Shanghai, China) following the manufacturer’s protocol. RNA quantification was performed using the Qubit RNA Assay Kit in Qubit 2.0 Fluorometer. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system. After the Illumina sequencing libraries were established, cDNA libraries were sequenced on the Illumina HiSeq platform. The mRNA-Seq was assembled and analyzed by the Guangzhou Gene Denovo Biotechnology Corporation (Guangzhou, China). The statistical power of this experimental design is 0.9; values for alpha and CV were 0.05, which were calculated in RNASeqPower. The effect parameter was 2. The sample size results at 6× sequencing depths were 7.55.

The low-quality reads were filtered using fastp (v0.19.3) with the following standard: the N content in any sequencing reads exceeded 10% of the base number of the read and the low quality (Q<=20) bases contained in reads exceeded 50%. The cleaned reads were aligned to the eggplant reference genome (eggplant genome consortium V3) (Barchi et al., 2019) using HISAT2 (version 2.1.0; http://daehwankimlab.github.io/hisat2/). The reads count of each gene was calculated by featureCounts (v1.6.2; Liao, Smyth & Shi, 2014). The differential expression analysis of two groups was performed using the DESeq2 R package (V1.22.1; Love, Huber & Anders, 2014). For identifying DEGs, absolute fold change ≥ 2 and false discovery rate (FDR) < 0.01 were used as screening criteria. The expression patterns of DEGs were analyzed using the Mfuzz R package (Kumar & Futschik, 2007). The heatmap was analyzed by the pheatmap R package.

Metabolite profiling using UPLC-MS/MS

The freeze-dried eggplant root was crushed using a mixer mill (MM 400; Retsch, Haan, Germany) with a zirconia bead for 1.5 min at 30 Hz. Approximately 100 mg lyophilized powder was dissolved 1.2 ml of 70% methanol solution. After the solution was vortexed for 30 s every 30 min six times, the solution was placed in a refrigerator at 4 °C overnight. Then the solutions were centrifuged at 12,000 rpm for 10 min and the supernatant fluid was filtered (SCAA-104, 0.22 μm pore size) before UPLC-MS/MS analysis.

UPLC conditions

Sample extracted supernatant fluid was separated using an UPLC system (UPLC, SHIMADZU Nexera X2). The chromatographic column was Agilent SB-C18 (1.8 µm, 2.1 mm × 100 mm) and the column temperature was at 40 °C. The flow rate was 0.35 ml/min with a total injection volume of 4 μl. The mobile phase A was pure water with 0.1% formic acid and the mobile phase B was acetonitrile with 0.1% formic acid. Within 0–9 min, the ratio of A to B was 95:5 → 5:95 and kept for 1 min. Within 10.00–11.10 min, the ratio of A to B was 5:95 → 95:5 and kept for 3 min. The effluent was alternatively connected to an ESI-triple quadrupole-linear ion trap (QTRAP)-MS.

ESI-Q TRAP-MS/MS

The effluent was analyzed by an AB4500 Q TRAP UPLC/MS/MS System which equipped with an ESI Turbo Ion-Spray interface. The ESI Turbo Ion-Spray interface was controlled by the Analyst 1.6.3 software (AB Sciex, Framingham, MA, USA). The ESI source operation parameters were as follows: the ion source was turbo spray; positive ions mode was 5,500 V and the negative ion mode was −4,500 V; the ion source gas I, gas II, and curtain gas were set to 50, 60, and 25.0 psi, respectively. Ten and 100 μmol/l polypropylene were used for instrument tuning and mass calibration for triple quadrupole and LIT modes, respectively. Triple quadrupole analysis was performed under a multiple reaction monitoring model with the collision gas (nitrogen) set to medium. After the declustering potential and collision energy optimization, the declustering potential and collision energy for individual multiple reaction monitoring transitions were obtained. Based on the metabolites eluted within each period, a specific set of multiple reaction monitoring transitions was monitored.

Data analysis

Quality control (QC) samples were prepared by mixing sample extracts. A QC sample was inserted into each of the 10 detected samples during the stability evaluation of analytical conditions. The metabolites was identified in accordance with the metware database and quantified by multiple reaction monitoring.

Significantly regulated metabolites between groups were determined by VIP ≥ 1 and absolute Log₂FC ≥ 1. VIP values were extracted from OPLS-DA results, which also contained score and permutation plots generated using the R package MetaboAnalystR (Chong & Xia, 2018). Data were subjected to log transformation (log2) and mean centering before OPLS-DA. A permutation test (200 permutations) was performed to avoid overfitting. The gene-metabolite correlation networks were analyzed by differential genes and differential metabolites with Pearson correlation coefficient > 0.80 and p-value < 0.05 in each pathway.

KEGG annotation and enrichment analyses

KEGG annotation of the identified metabolites was performed using the KEGG compound database (http://www.kegg.jp/kegg/compound/). Then, the annotated metabolites were mapped to the KEGG pathway database. The differential metabolites were classified according to the types of pathways in KEGG according to the annotation of KEGG.

Analysis of the bacterial wilt resistance of eggplant

After the R and S genotypes of eggplant material were inoculated with the GMI1000-tac-EGFP strain, the DI and EGFP fluorescence were analyzed. The results showed that at 10 days after inoculation with R. solanacearum, the DIs of R and S were 0 and 2.48, respectively (Fig. 1A). R genotypes showed normal phenotype, whereas S genotypes showed wilt (Fig. 1B). EGFP fluorescence was detected in S stem and root but was not observed in the R stem and root (Figs. 1C and 1D). This result showed that R genotypes were highly resistant to GMI1000.

Induced responses to bacterial wilt in global transcriptome changes of eggplant

The transcriptome was compared after inoculation of GMI1000-tac-EGFP to understand the mechanism of the bacterial wilt defense response of eggplant. Three time points (i.e., 0, 24, and 48 h postinoculation (hpi)) were analyzed.

Approximately, 390.98 million cleaned reads were generated for nine samples (Table 1). About 85% cleaned reads were aligned to the eggplant reference genome. Transcriptomic sequences were deposited in the NCBI Sequence Read Archive under accession number PRJNA837016.

PCA showed that the first two PCAs explained 59.48% of the total variation (Fig. 2A). The heatmap of DEGs showed a significant difference in gene expression level after the inoculation of R. solanacearum (Fig. 2B). After filtration by FDR < 0.01 and absolute Log2 (fold change (FC)) ≥ 1, 1,831 (799 upregulated and 1,032 downregulated), 1,416 (708 upregulated and 708 downregulated), and 1,032 (538 upregulated and 494 downregulated) DEGs were identified in R-0h_vs_R-24h, R-0h_vs_R-48h, and R-24h_vs_R-48h, respectively (Fig. 2C). These DEGs are listed in Tables S1–S3. After the DEGs were selected as absolute Log₂ > 5, 27 DEGs were selected. In the R-0h_Vs_R-24h group, eight genes were unregulated and nine genes were downregulated. In R-0h_Vs_ R-48h, the 11 DEGs were downregulated. But in the R-24h_Vs_ R48-h, there were only one downregulated gene (Table S4).

A total of 485 and 751 genes were upregulated and downregulated, respectively, in at least one time point. Four genes were common in the R-0h_vs_R-24hUp and R-0h_vs_R-48hDn sets, and a gene was common in the R-0h_vs_R-24hDn and R-0h_vs_R-48hUp sets (Fig. 3).

KEGG and KOG classification of DEGs

DEGs were mapped to the KEGG pathway to understand DEGs function in the eggplant defense response. The top 20 pathways included metabolic pathways, biosynthesis of secondary metabolites, plant hormone signal transduction, MAPK signaling pathway, plant–pathogen interaction, and flavonoid biosynthesis pathway (Fig. S1).

Expression pattern analysis of DEGs

A set of genes with similar expression patterns was functionally correlated. Six expression patterns were obtained in accordance with the DEGs data. The expression pattern of the centers of cluster 1 was upregulated at 24 and 48 hpi. The expression pattern of the centers of cluster 2 was downregulated at 24 hpi and then upregulated at 48 hpi. Ultimately, the expression level at 0 and 48 hpi were not significantly different. The expression pattern of the centers of cluster 3 was upregulated at 24 hpi and then downregulated at 48 hpi. Ultimately, the expression level at 0 and 48 hpi was not ignificantly different. The expression pattern of the centers of cluster 4 was upregulated at 24 hpi and then downregulated at 48 hpi. However, the expression level at 48 hpi was significantly increased compared with that at 0 hpi. The expression pattern of the centers of cluster 5 was downregulated at 24 hpi and then upregulated at 48 hpi. However, the expression level at 48 hpi was significantly decreased compared with that at 0 hpi. The expression pattern of the centers of cluster 6 was downregulated at 24 and 48 hpi (Fig. 4 and Table S5). A total of 488 (cluster 1) and 443 (cluster 6) DEGs were maintained to be upregulated and downregulated at 24 and 48 hpi. In cluster 1, there were seven WRKY transcriptional factors, four TIFY transcriptional factors and three bHLH transcriptional factors. Meanwhile, there were seven resistance genes were identified (Table 2).

Widely targeted metabolome analysis

On the basis of UPLC-MS/MS and the metware metabolite database, 661 metabolites were detected (Table S6). The PCA results showed that the first two components could explain 42.23% of the dataset variation (Fig. 5A). Cross-validation indicated the first two components relevant for the classification of variation, which illustrated different directions of response to R. solanacearum. The heatmap showed different expression profiles after inoculation with R. solanacearum (Fig. 5B). A total of 44 (11 upregulated and 33 upregulated), 25 (six upregulated and 19 upregulated), and 24 (14 upregulated and 10 upregulated) differential metabolites were identified in R-0h_vs_R-24h, R-0h_vs_R-48h, and R-24h_vs_R-48h, respectively (Fig. 5C). Differential metabolites are listed in Tables S7–S9.

A total of 15 and 41 metabolites were upregulated and downregulated, respectively, in at least one time point. A total of two and 11 metabolites were upregulated and downregulated, respectively, at R-24h and R-48h (Fig. 6).

KEGG of differential metabolites

KEGG classification results showed that the top pathways were metabolic pathways and biosynthesis of secondary metabolites (Fig. S2). Several metabolites were involved in phenylpropanoid biosynthesis, flavone and flavonol biosynthesis, and plant hormone signal transduction pathway.

Expression pattern analysis of differential metabolites

A set of metabolites with similar expression patterns was functionally correlated. Six expression patterns were obtained in accordance with differential metabolites data. The expression pattern of the centers of cluster 1 was downregulated at 24 hpi and then upregulated at 48 hpi. However, the expression level at 48 hpi was significantly decreased compared with that at 0 hpi. The expression pattern of the centers of cluster 2 was upregulated at 24 and 48 hpi. The expression pattern of centers of cluster 3 was upregulated at 24 hpi and then significantly downregulated at 48 hpi. The expression level at 48 hpi was significantly decreased compared with that at 0 hpi. The expression pattern of the centers of cluster 4 was downregulated at 24 and 48 hpi. The expression pattern of the centers of cluster 5 was downregulated at 24 hpi and then upregulated at 48 hpi. Ultimately, the expression level at 48 hpi was significantly increased compared than at 0 hpi. The expression pattern of the centers of cluster 6 was downregulated at 24 hpi and at 48 hpi. Ultimately, the expression level at 24 and 48 hpi were not significantly different (Fig. 7 and Table S10). The (−)-Jasmonoyl-L-Isoleucine, the 5′-Glucosyloxyjasmanic acid and 1-O-Salicyl-D-glucose were upregulated at R-24h and R-48h. The salicylic acid was downregulated at R-24h and R-48h.

Integration analysis of transcriptomic and metabolic datasets

DEGs and differential metabolites were simultaneously assigned to KEGG pathways (p < 0.05) to understand the eggplant response to bacterial wilt. The results showed that only alpha-linolenic acid metabolism (ko00592) and plant hormone signal transduction pathway (ko04075) were significantly enriched in the R-0h_vs_R-24h group (Fig. S3). The metabolites involved in these two pathways were JA, ABA, jasmonate, and 9-Hydroxy-12-oxo- 15(Z)-octadecenoic acids. However, 69 genes were involved in these two pathways (Table S9). The ko00592 pathway was the precursor of the JA biosynthesis pathway which ultimately biosynthesizes methyl-jasmonate.

Transcriptome and metabolome data were also compared by the Pearson correlation analysis (Pearson correlation coefficient > 0.8). Gene–metabolite correlation networks were also constructed (Fig. 8). In the R-0h_vs_R-24h group, SMEL_008g298210.1,SMEL_011g363600.1, and SMEL_005g223990.1 were highly correlated with JA. SMEL_003g183930.1 and SMEL_006g251030.1 were highly correlated with ABA. In the R-0h_vs_R-48h group, 13 genes were highly correlated with (−)-Jasmonoyl-L-Isoleucine. These results suggested that JA might regulate gene expression in the bacterial wilt resistance defense of eggplant.