Transcriptomic and metabolomic differences between banana varieties which are resistant or susceptible to Fusarium wilt

Plant material cultivation and Foc4 treatment

In this study, Guijiao 1 (G1, Fusarium wilt susceptible variety) and Guijiao 9 (G9, Fusarium wilt moderately resistant banana variety) were used to performed the transcriptomics and metabolisms experiments. We took the banana seedlings out of the nutrient cups and rinsed them three to five times with tap water to wash the culture substrate, then three to four times with distilled water, and finally three to four times with ultrapure water. We soaked the seedlings in thymol solution (5 mg/mL) for 3 min and then rinsed them with ultrapure water three times.

We placed the banana seedlings in a container that had been sterilized in advance, added 800 mL of Hoagland liquid medium and 50 mL of 5 mmol/L calcium chloride solution, and then added an appropriate amount of ceramist to cover the bare roots. All the plants were placed in a greenhouse and allowed to grow under natural light conditions. After 1 day of culture, the first sample was taken as the control group.

The banana seedlings were removed from Hoagland liquid medium, immersed in 800 mL of Foc4 spore suspension (Hoagland liquid medium) containing spores at 1 × 10⁷·ml⁻¹ concentration for 2–2.5 h, and then returned to the original Hoagland liquid medium (containing ceramsite). The second sample was taken on the third day after inoculation with Foc4. The third sample was taken on the seventh days after inoculation.

Root tissues (about 1 cm of the white root tip) were collected and washed quickly with PBS buffer to remove surface impurities of the root, flash-frozen with liquid nitrogen, placed in 1.5 mL tubes precooled in liquid nitrogen, and stored at −80 °C in a refrigerator for RNA extraction, library building, and transcriptome sequencing. The sample was ground in a homogenizer and then added to Trizol reagent (Invitrogen, Waltham, MA, USA). After vortex mixing, chloroform was added and centrifuged at low temperature. The RNA was precipitated with ethanol and dissolved in Rnase free water. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). RNA concentration was measured using Qubit^® RNA Assay Kit in Qubit^® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). The sample size estimation was estimated by using RnaSeqSampleSize (Zhao et al., 2018).

The 20 mL root culture medium was placed in a pre-cooled 50 mL centrifuge tube, flash-frozen in liquid nitrogen, and then stored in a refrigerator at −80 °C. It was used to identify the metabolome of root exudates. After the samples were slowly thawed at 4 °C, appropriate amount of samples were added to pre-cooled methanol/acetonitrile/aqueous solution (2:2: 1, v/v), vortex mixing, ultrasound at low temperature for 30 min, stand at −20 °C for 10 min, centrifuge at 14,000 g and 4 °C for 20 min, vacuum drying with supernatant, adding 100 μL acetonitrile solution during mass spectrometry (acetonitrile: Water = 1:1, v/v) dissolved, swirled, centrifuged 14,000 g at °C for 15 min, and the supernatant was taken for analysis. Transcriptome and metabolome samples were taken at the same time. Three biological replicates were collected for each treatment.

Sequencing and transcriptome data analysis

The RNA-Seq libraries were prepared using the Illumina TruSeq RNA Sample Pre-Kit following the manufacturer’s protocols. The sequencing had been performed in paired-end reads (2 × 151 bp) using the Illumina HiSeq sequencing platform (Beijing Biomics Biotech Co. Ltd., Beijing, China). The RNA-Seq data have been deposited at GenBank under the accession PRJNA916078.

Trimmomatic (version 0.38) was used to remove the adaptor sequences and the low-quality (<Q20) bases at the 5′ and 3′ ends (Bolger, Marc & Bjoern, 2014), and reads longer than 70 bp were used for further experiment. The reads were mapped to the banana genome (D’Hont et al., 2012) using bowtie2 (version 2.1.0) (Langmead & Salzberg, 2012) with default parameters after preprocessing of RNA-Seq data. TopHat2 (version 2.2.1) was also used to perform a sequence comparison between clean reads and the reference genome to obtain position information on the reference genome or gene and sequence characteristic information unique to the sequencing sample. Gene expression levels were presented as FPKM (fragments per kilobase of transcript per million fragments mapped) values calculated by using DESeq2 (version 3.17). Genes with expression levels >1 FPKM were retained for further analysis. We used FDR < 0.05 and more than two-fold change as the criteria to classify DEGs between two compared pairs after the pairs were processed by DEseq2 (version 3.17, www.bioconductor.org). The resulting p values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate.

AmiGO (version 2.0) with the default parameters was used to obtain the gene ontology terms of each gene and analyze GO (Gene Ontology) functional enrichment by using hypergeometric tests with FDR (false discovery rate) correction to obtain an adjusted p < 0.01 between test gene groups and the whole annotation data set, respectively. The differentially expressed genes in the KEGG pathway were analyzed using Cytoscape (version 3.10.0) (Ideker, 2011) with the ClueGO plugin (Wegdam et al., 2014).

Identification and data analysis of untargeted metabolomics

Samples were separated on an ultra-high-performance liquid chromatography system (UHPLC, Agilent 1290 Infinity LC) HILIC column. The column temperature was 25 °C, the flow rate was 0.5 mL/min, and the injection volume was 2 μl. The mobile phase consisted of two phases: A (water + 25 mM ammonium acetate + 25 mM ammonia water) and B (acetonitrile). The gradient elution procedure was as follows: 0–0.5 min, 95% B; 0.5–7min, B changes linearly from 95% to 65%; 7–8 min, B changes linearly from 65% to 40%; 8–9 min, B maintained at 40%; 9–9.1 min, B changes linearly from 40% to 95%; 9.1–12 min, B maintained at 95%. During the entire analysis, the samples were placed in an autosampler at 4 °C.

To avoid the influence of the fluctuating instrument detection signal, the samples were analyzed continuously in random order. Quality control (QC) samples were inserted into the sample cohort to monitor and evaluate the stability of the system and the reliability of the experimental data.

The AB Triple TOF 6,600 Mass Spectrometer (AB SCIEX) was used to perform mass spectrometry analysis. The ESI Source conditions after HILIC chromatographic separation were as follows: Ion Source Gas1, Gas2 and Curtain Gas (CUR) was of 60, 60 and 30, respectively; source temperature was of 600 °C; IonSapary Voltage Floating (ISVF) was of ± 5,500 V (both positive and negative modes), TOFMS scan M/Z range from 60 to1,000 Da; product ION scan M/Z range from 25 to1,000 Da; TOF MS scan retention time was of 0.20 s/spectra, and product ion scan retention time was 0.05 s/spectra. To acquire high resolution mass spectra, the secondary mass spectra data obtained were obtained using information-dependent acquisition (IDA) of ESI positive ion mode and the ESI negative ion mode. The declustering potential were following as in situ: ±60 v. Frontiers: 35 ± 15 eV. The IDA was set as follows: ex situ within 4 Da and candidate ions to monitor per cycle: 10.

The original data was converted into MzXML format by ProteoWizard, and then XCMS software was used for peak alignment, retention time correction and peak area extraction. Firstly, the metabolite structure identification and data preprocessing were carried out for the data extracted by XCMS, and then the experimental data quality evaluation was carried out, and finally the data analysis was carried out. Peak alignment, retention time correction, and peak area extraction were performed by MSDAIL software. Metabolite structure identification, data preprocessing, experimental quality evaluation, and data analysis were carried out for data extracted from MSDAIL. Principal component analysis (PCA) was performed in SIMCA-P 14.1 software to investigate the differences between biological replicates of different samples. Orthogonal partial least squares discriminant analysis (OPLS-DA) was used to measure the differences between the treatment groups. The differential metabolites were screened according to the contribution value (VIP ≥ 1), the significance of the change between groups (p < 0.05), and the fold change (fold change ≥ 2).

The data reported in this article have been deposited in the OMIX, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (https://ngdc.cncb.ac.cn/omix: accession no. OMIX004812).

Identification of Foc4-responsive genes

In the Foc4-susceptible variety (G1), total of 76 genes showed decreased expression and 12 genes showed increased expression after 3 days of Foc4 inoculation, and 71 genes showed decreased expression and 17 genes showed increased expression on seventh day after inoculation. In all, 11 genes were up-regulated and 11 genes were down-regulated on seventh day compared with those on third day after inoculation (Fig. 1A, S2, and Tables S5–S7).

For the Foc4-resistant variety (G9), the expression of 39 genes decreased and 98 genes increased 3 days after inoculation with Foc4 and the expression of 290 genes decreased and 129 genes increased on seventh day after inoculation. However, 1,282 genes were down-regulated and only 268 genes were up-regulated on seventh day compared with those on the 3rd day after inoculation (Figs. 1B, S3, and Tables S8–S10).

The expression levels of 16 genes in the Foc4-resistant variety (G9) were lower than those in the Foc4-susceptible variety (G1) even before pathogen treatment. On the third day after inoculation, the expression levels of 47 genes in G9 were lower than those in G1 and the expression levels of 10 genes were higher than those in G1 (Figs. 1C, 2A, and Table S11). On the seventh day after inoculation with Foc4, the expression levels of 644 genes in G9 were significantly lower than those in G1 and only the expression levels of 99 genes were higher than those in G1 (Figs. 1C, 2B, and Tables S12 and S13).

Function annotation of Foc4-responsive genes

To fully analyze the functions of DEGs, we performed GO and KEGG functional enrichment analysis (Figs. 3, 4, S4–S7, and Tables S5–S13). In G1, 3 days after inoculation (Fig. S4A), the DEGs were mainly enriched in the oxidation-reduction process (eight), heme binding (four), the extracellular region (four), the hydrogen peroxide catabolic process (three), and the response to oxidative stress (three) and after 7 days of inoculation (Fig. S4B), the DEGs were mainly enriched in the extracellular region (six), the integral component of the membrane (10), and the cell wall (three) in GO terms.

In G9, 3 days after inoculation (Fig. S5A), the DEGs were mainly enriched in cell wall organization (four), the xyloglucan metabolic process (three), and cell wall biogenesis (three) in GO terms. Compared with the control, 7 days after inoculation, in G9 (Fig. S5B), the DEGs were mainly enriched in the oxidation-reduction process (45), the negative regulation of endopeptidase activity (six), cell wall organization (seven), oxidoreductase activity (14), monooxygenase activity (eight), iron ion binding (12), hydrolase activity (seven), heme binding (12), iron ion binding (12), the extracellular region (14), and the cell wall (seven) in GO terms.

However, compared with three and 7 days after inoculation (Fig. S5C), the DEGs were mainly enriched in the oxidation-reduction process (72), the metabolic process (19), nucleosome assembly (eight), carbohydrate transport (seven), the integral component of the membrane (155), oxidoreductase activity (23), potassium ion transmembrane transporter activity (five), nucleosome (11), and the glucose metabolic process (four). Under normal growth conditions, in G1 and G9, the DEGs were enriched in the membrane (two), methylation (one), the metabolic process (one), and terpene synthase activity (one) (Fig. 3A). In GO terms of 3 days after inoculation, DEGs were enriched in the terpenoid biosynthetic process (two), cell wall organization (two), translation initiation factor activity (two), hydrolase activity (two), the apoplast (two), the membrane (three), the cell wall (two), and catalytic activity (two) (Fig. 3B). After 7 days of inoculation, DEGs were mainly enriched in the extracellular region (25), kinase activity (six), L-lactate dehydrogenase activity (two), methionine adenosyl transferase activity (three), hydrolase activity, acting on ester bonds (eight), metal ion binding (37), phenylalanine ammonia-lyase activity (four), peroxidase activity (13), heme binding (23), the steroid biosynthetic process (three), the S-adenosylmethionine biosynthetic process (three), plant-type cell wall organization (five), the cinnamic acid biosynthetic process (four), the oxidation-reduction process (48), the L-phenylalanine catabolic process (five), cellular oxidant detoxification (13), the response to oxidative stress (13), the hydrogen peroxide catabolic process (13), and the cell wall (nine) (Fig. 3C).

In G1 of 3 days after inoculation (Fig. S6A), the main DEG-enriched KEGG pathways identified were phenylpropanoid biosynthesis; nitrogen metabolism; alanine, aspartate, and glutamate metabolism; and circadian rhythm. After 7 days, phenylpropanoid biosynthesis and plant-pathogen interaction were the mainly enriched KEGG pathways (Fig. S6B).

In G9 of 3 days after inoculation (Fig. S7A), the main DEG-enriched KEGG pathways identified were ubiquitin-mediated proteolysis, the Wnt signaling pathway, plant hormone signal transduction, nitrogen metabolism, and MAPK signaling. After 7 days of inoculation (Fig. S7B), the main DEG-enriched KEGG pathways identified were glycolysis/gluconeogenesis, nitrogen metabolism, starch and sucrose metabolism, and plant hormone signal transduction.

Under normal growth conditions, DEGs between G1 and G9 (Fig. 4A) were enriched in glycerophospholipid metabolism, glycerolipid metabolism and cutin, and suberine and wax biosynthesis KEGG pathways. After 3 days of inoculation, the main DEG-enriched KEGG pathways were glycolysis/gluconeogenesis and terpenoid backbone biosynthesis. After 7 days of inoculation, the main DEG-enriched KEGG pathways were plant-pathogen interaction and phenylpropanoid biosynthesis (Fig. 4B).

Protein-protein interaction network analysis of Foc4-responsive genes

We constructed an interaction network based on the DEGs retrieved from known PPI databases (Figs. 5, S8 and S9). Regarding G1, when comparing the control and the inoculated G1 variety on third day after inoculation, six DEGs could be linked via known nitrogen metabolism relationships. Significantly, the Macma4_09_g15700 (4 link and mainly down-regulated) encoding nitrate reductase was proposed as a potential hub (Fig. S8).

As for G9, when comparing the control and the inoculated variety on third day after inoculation, eight DEGs could be linked via known nitrogen metabolism relationships. Macma4_09_g15700 (six link and down-regulated) was also proposed as a potential hub (Fig. S9A). On comparing the control and the inoculated variety 7 days after inoculation (Fig. S9B), 43 DEGs could be linked via energy production and conversion, inorganic ion transport and metabolism, and other relationships. Macma4_09_g15510 (nine link and down-regulated) and Macma4_09_g15700 both encoding nitrate reductase were also proposed as potential hubs. In addition, Macma4_10_g32300 (six link and down-regulated) encoding glucose-6-phosphate isomerase, Macma4_07_g06890 (nine link and down-regulated) encoding L-lactate dehydrogenase, and Macma4_09_g25660 pyruvate kinase (nine link and down-regulated) were proposed as the key hub proteins.

When comparing the inoculated varieties on third day after inoculation with those on seventh day after inoculation, 410 DEGs could be linked via carbohydrate transport and metabolism, energy production and conversion, and other relationships. Among them, Macma4_11_g21890 (11 link and up-regulated) encoding enolase, Macma4_07_g06890 (10 link and down-regulated), Macma4_02_g23780 (10 link and down-regulated) encoding pyruvate kinase, and Macma4_07_g06880 (10 link and down-regulated) encoding L-lactate dehydrogenase were key potential hubs (Fig. S9C).

Comparing G1 and G9 after inoculation for 7 days, total of 65 DEGs could be linked via energy production and conversion, oxidation-reduction reactions, secondary metabolites biosynthesis, and lipid metabolism relationships (Fig. 5).

Metabolite identification results

In this study, the differences between the metabolites in the two varieties after Foc4 inoculation were analyzed using UPLC-QTOF-MS. The results showed that 400 metabolites were matched with the database in the positive ion mode (Table S14), 237 metabolites were matched in the negative ion mode (Table S15), and 611 metabolites were combined. There are 12 super classes (Table S16) and 82 classes (Table S17).

The peaks extracted from all the experimental samples and QC samples were analyzed by PCA, as shown in Fig. S10. The results showed that QC samples in positive and negative ion modes are closely clustered together, indicating that the experiment has good repeatability. In the OPLS-DA displacement test, the samples of the comparison groups showed an obvious trend of separation, indicated that the model had a good predictive ability and was effective (Fig. S11).

Differential metabolite screening

According to metabolomics identification, without Foc4 inoculation, there were eight differential metabolites between G9 and G1 in the positive ion mode and two differential metabolites between G9 and G1 in the negative ion mode (Fig. 6). A total of 3 days after Foc4 inoculation, there were four differential metabolites between G9 and G1 in the positive ion mode and five differential metabolites between G9 and G1 in the negative ion mode. On the seventh day after Foc4 inoculation, there were eight differential metabolites between G9 and G1 in the positive ion mode and nine differential metabolites between G9 and G1 in the negative ion mode.

All the metabolites detected in positive and negative ion modes (including those not identified) were analyzed for differences based on univariate analysis. When the difference between groups was FC > 1.5 or FC < 0.67 and the p value < 0.05, the metabolites were defined as differential metabolites and super class classification statistics were performed for the differential metabolites (Figs. 7–9).

In G9 and G1, without Foc4 inoculation, the up-regulated metabolites identified in the negative ion mode mainly belonged to lipids and lipid-like molecules, phenylpropanoids and polyketides, and organic acids and derivatives (Fig. 7). The metabolites identified in the positive ion mode mainly belonged to alkaloids and derivatives, benzenoids, lipids and lipid-like molecules, organic acids and derivatives, organic nitrogen compounds, and phenylpropanoids and polyketides. The up-regulated metabolites were mainly of the organ heterocyclic compounds.

After 3 days of Foc4 inoculation, the differential metabolites identified in positive and negative ion modes were similar in classification and the up-regulated metabolites mainly belonged to benzenoids, lipids and lipid-like molecules, organic acids and derivatives, the organ heterocyclic compounds, and phenylpropanoids and polyketides (Fig. 8).

After 7 days of Foc4 inoculation, the up-regulated metabolites identified in the negative ion mode were mainly benzenoids, lipids and lipid-like molecules, organic oxygen compounds, phenylpropanoids and polyketides (Fig. 9). The metabolites identified in the positive ion mode mainly belonged to lipids and lipid-like molecules, organic oxygen compounds, and phenylpropanoids and polyketides. The up-regulated metabolites mainly belonged to benzenoids and nucleosides, nucleotides, and analogues.