Control effect of CEF-082 on Verticillium wilt of cotton and H2O2 content
The disease index was 18.61 in the control group (water+ V. dahliae) and 7.62 in the treatment group (CEF-082+ V. dahliae) 14 d after V. dahliae inoculation (Fig. 1A). The results showed that CEF-082 enhanced the resistance of cotton to Verticillium wilt, and the biocontrol effect was 59.1% (Fig. 1C).
The H2O2 content in the treatment group was higher than that in the control group throughout the majority of the duration of the experiment and lower than that in the control group at 5 dpi with V. dahliae. The H2O2 content in the treatment group was highest at 2 dpi (12.80 μmol/g), while the H2O2 content in the control group was highest at 1 dpi (10.38 μmol/g). The changes in the two groups were similar and were stable 5 d later (Fig. 1B).
Verification of RNA-Seq analysis by qRT-PCR
Twelve DEGs were randomly selected. The gene expression levels in the control and treatment groups were compared by qRT-PCR. The RNA-seq data showed that the expression of the 12 genes was upregulated at 0 h, 12 h or 48 h. The qRT-PCR results showed that the expression of nine of the 12 genes was upregulated, which was consistent with the results of their upregulated expression in the transcriptome; however, the expression of three genes was downregulated, which was inconsistent with their expression in the transcriptome, namely, Gh_D12G2793, Gh_D08G2484 and Gh_D05G3615 (Fig. 2). In addition, the levels of upregulation of 5 genes in the qRT-PCR data were lower than those in the RNA-seq data. The qRT-PCR data were up to 75% consistent with the transcriptome data.
Functional annotation and enrichment analysis of the DEGs
The minimum correlation between the three replicates was 95.5% (Fig. S1). Principal component analysis (PCA) of 18 arrays (Fig. S2) was also used to compare the samples and to explore the dynamic changes in the cotton transcriptome after treatment with CEF-082 and V. dahliae.
The average clean reads of the 18 samples was 62.08 M. The lowest Q20 value of the clean reads was 97.93, and the lowest Q30 value was 90.06 (Table S2). A total of 47183 new transcripts were found, of which 7288 belonged to new protein-coding genes (Table S3).
There were 3480 upregulated and 2158 downregulated DEGs at 0 h, 1716 upregulated and 1205 downregulated DEGs at 12 h, and 1524 upregulated and 629 downregulated DEGs at 48 h. The greatest number of DEGs was identified after inoculation with CEF-082 for 24 h. After inoculation with V. dahliae, the number of DEGs gradually decreased.
DEGs induced by CEF-082 at 0 h without V. dahliae inoculation
After inoculation with CEF-082 for 24 h (0 h), 5638 DEGs were identified, and KEGG pathway enrichment analysis revealed 15 significantly enriched pathways, including plant-pathogen interaction, MAPK signalling pathway-plant, flavonoid biosynthesis, phenylpropanoid biosynthesis, galactose metabolism, arachidonic acid metabolism, carotenoid biosynthesis, glutathione metabolism, sesquiterpenoid and triterpenoid biosynthesis, linoleic acid metabolism, other glycan degradation, glycosphingolipid biosynthesis - ganglio series, brassinosteroid biosynthesis, diterpenoid biosynthesis and sphingolipid metabolism (Q-value <0.05) (Table 1). In the plant-pathogen interaction pathway, there were 106 FLS2 genes, 88 upregulated and 18 downregulated; 7 Rboh genes, 5 upregulated and 2 downregulated; 5 upregulated calcium-dependent protein kinase (CDPK) genes; 5 CNGC genes, 3 upregulated and 2 downregulated; and 57 glutathione S-transferase (GST) genes in the glutathione metabolism pathway, of which 49 were upregulated and 8 downregulated (Fig. 3). These genes were related to the metabolism of reactive oxygen species (ROS) and Ca2+. In the MAPK signalling pathway-plant pathway, 304 DEGs regulated 30 crucial points related to ROS, Ca2+, abscisic acid (ABA), ethylene (ET), jasmonic acid (JA), H2O2 and FLS2. In the flavonoid biosynthesis pathway, the genes encoding chalcone synthase (CHS) and ferulate-5-hydroxylase (F5H) were induced. In the phenylpropanoid biosynthesis pathway, the key genes PAL and 4CL were also induced.
The GO enrichment analysis revealed that the 5638 genes were mainly enriched in 86 terms, including the intrinsic component of membrane, integral component of membrane, membrane part, membrane, catalytic activity, response to biotic stimulus, cell wall, oxidoreductase activity, defense response, response to stimulus, response to stress, and response to fungus (Q-value <0.001), and the first 15 terms are listed in Table 2. Of the 16 genes in the response to fungus term, 15 were upregulated and 1 was downregulated. The GO classification showed that there were 18, 14 and 12 terms in biological process, cellular component and molecular function, respectively, and the KEGG classification indicated that the DEGs mainly belonged to the metabolism pathway (2856 DEGs).
DEGs coinduced by CEF-082 and V. dahliae
The 463 shared DEGs at 12 h and 48 h were significantly enriched in 6 KEGG pathways (Table 3). In the plant-pathogen interaction pathway, 29 DEGs regulated 8 crucial points, including CNGCs, calmodulin (CaM), FLS2, disease resistance protein RPS2 (RPS2), heat shock protein 90 kDa (HSP90), pto-interacting protein 1 (Pti1), disease resistance protein RPM1 (RPM1), and EIX receptor 1/2 (EIX1/2). In the phenylpropanoid biosynthesis pathway, 23 DEGs regulated 9 crucial points. In the flavonoid biosynthesis pathway, 12 DEGs regulated 8 crucial points. The enriched GO terms included terpenoid metabolic process, oxidoreductase activity, defense response, H2O2 metabolic process and ROS metabolic process terms.
DEGs induced only in cotton inoculated with V. dahliae in the presence of CEF-082
A total of 1209 specific DEGs were identified at 12 h and 48 h, which were induced only in cotton plants inoculated with V. dahliae in the presence of CEF-082, but not when cotton plants were inoculated with V. dahliae only. The cluster thermogram showed the expression patterns of these genes at different stages (Fig. S3). KEGG classification showed that these DEGs mainly belonged to metabolism (672 DEGs) and were significantly enriched in 5 KEGG pathways, including flavonoid biosynthesis, indole alkaloid biosynthesis, MAPK signalling pathway-plant, plant-pathogen interaction, and phenylpropanoid biosynthesis (Table 4). GO classification showed that there were 14, 12 and 9 terms in the biological process, cellular component and molecular function, respectively. GO enrichment indicated that these DEGs were enriched in ROS metabolic process (14 DEGs), H2O2 metabolic process (12 DEGs), H2O2 catabolic process (12 DEGs), defense response (31 DEGs), superoxide dismutase activity (5 DEGs), antioxidant activity (19 DEGs), oxidoreductase activity, acting on superoxide radicals as acceptor (5 DEGs), cofactor binding (75 DEGs) and DNA binding (121 DEGs) (Fig. S4).
At 12 h and 48 h, 96 shared DEGs were obtained, which were induced only in cotton plants inoculated with V. dahliae in the presence of CEF-082, but not when cotton plants were inoculated with V. dahliae only (Fig. S5). KEGG analysis of the 96 DEGs indicated that they were mainly enriched in glutathione metabolism and flavonoid biosynthesis (Table 5). GO analysis showed that the DEGs were enriched in the terms superoxide dismutase activity, oxidoreductase activity, acting on superoxide radicals as acceptors, and antioxidant activity. Of the 96 DEGs, 9 encoded TFs and 20 encoded predicted PRGs (Table S4).
A protein-protein interaction network (Fig. S6) was constructed via the 96 DEGs shared between 12 h and 48 h and other genes interacting with them in cotton. Six hub genes were obtained, Gh_A05G1020, Gh_D09G0858, BGI_novel_G004376, Gh_A08G0125, Gh_D07G1197 and Gh_A05G3508. Among them, Gh_D07G1197 was annotated in the flavonoid biosynthesis pathway.
Putative R genes and TFs involved in resistance to Verticillium wilt
On the basis of the transcriptome analysis, a total of 65 candidate genes that may be related to the resistance of cotton to Verticillium wilt were identified, including 5 CNLs (whose members contain an NB-ARC domain), 3 CNs (members of the U-box domain-containing protein kinase family protein), 5 NLs (whose members contain an NBS-LRR domain), 7 RLPs (whose members contain an eLRR-TM-S/TPK domain), 7 Ns (whose members contain an NBS domain only), 9 TNLs (members of the TIR-NBS-LRR class), 6 Ts (members of NAC domain containing protein 17), 1 Mlo-like (a member of the Mlo-like resistant proteins) and 2 other types (which have resistance functions but do not fit the known classes). These genes mainly included a disease resistance protein, 2 probable calcium-binding protein (CML45), 3 ethylene-responsive transcription factor (ERF), 2 cyclic nucleotide-gated ion channel 2 (CNGC2), 5 MYB TFs and 2 GST (Table 6-1, Table 6-2, and Table 6-3). A clustering thermogram of 65 genes (Fig. 4) showed that certain genes were upregulated at 0, 12 and 48 h; certain genes were downregulated at 0 h and upregulated at 12 and 48 h; and certain genes were downregulated at 0, 12 and 48 h.