Sequencing overview
A total of 45.88 GB clean data were generated by sequencing platform, and clean data number of each sample was more than 6.23 GB. GC content ranged from 47.16% to 49.09%, and Q30 of each sample was above 92.92% (Additional file: Table S1). The sequencing results showed that sequencing fragments had high randomness and reliability (Additional file: Figure S1A). After transcript de novo assembly, 60324 Unigenes in total were obtained, and the N50 was 2409 kb. Furthermore, 19670 (32.61%) of them were over 1 kb in length (Additional file: Figure S1B). All these indicative data display high assembly integrity.
Functional annotation and differential expression analysis
DEGs annotation and function classification
A total of 5140 DEGs were annotated to seven databases. nr (RefSeq non-redundant proteins) had the largest number of annotated DEGs (5066), while KEGG had the least (1745) (Fig.1a). The venn diagram of function annotation in different database as Fig.1b. It was learned that those DEGs between healthy and fungal diseased samples chiefly classified into “signal transduction mechanisms”, “carbohydrate transport and metabolism”, “defense mechanisms”, “energy production and conversion”, “general function prediction only”, “posttranslation modification, protein turnover, chaperones”, “translation, ribosomal structure and biogenesis” (Fig.1c, d). We can believe that G. elata Bl. f. glauca possibly respond to fungal disease by enhancing or weakening these physiological or biological activities.
GO enrichment and KEGG enrichment analysis
Using GO database, 2482 DEGs were enriched into 3958 GO terms. GO terms are usually classified into 3 categories: biological process (BP), cellular component (CC), molecular function (MF). Here, 2363 (59.70%) of these GO terms fell into BP, 509 (1.49%) belong to CC, and 1086 (27.44%) were part of MF. Noteworthily, 36 GO terms were about signal transduction, and 24 GO terms were relevant to hormone. By Kolmogorov-Smirnov test, 421 GO terms were significantly enriched (P-value<0.05). Part of them were showed in Additional file: Table S2 and top 30 were displayed as Fig.2a.
Using KEGG database, 122 pathways were enriched and top 50 was showed as Figure 2c. The enrichment degree was based on the P-value and enrichment factor (Fig. 2b). 9 pathways were significantly enriched (p<0.05), and they fell into 3 pathway categories: metabolism, environmental information processing, organismal systems. Specific pathway names displayed in Table 1.
Table 1 KEGG pathway enrichment analysis (p<0.05)
Pathway category
|
Pathway description
|
Specific pathway
|
ko ID
|
DEG
|
All Unigene
|
P-value
|
Metabolism
|
Carbohydrate metabolism
|
Starch and sucrose metabolism
|
ko00500
|
49
|
169
|
0.041
|
Metabolism of cofactors and vitamins
|
Ubiquinone and other terpenoid-quinone biosynthesis
|
ko00130
|
12
|
32
|
0.047
|
Metabolism of terpenoids and polyketides
|
Brassinosteroid biosynthesis
|
ko00905
|
7
|
11
|
0.005
|
Diterpenoid biosynthesis
|
ko00904
|
7
|
12
|
0.009
|
Biosynthesis of other secondary metabolites
|
Flavone and flavonol biosynthesis
|
ko00944
|
5
|
5
|
0.001
|
Phenylpropanoid biosynthesis
|
ko00940
|
36
|
96
|
0.001
|
Flavonoid biosynthesis
|
ko00941
|
14
|
32
|
0.008
|
Environmental Information Processing
|
Signal transduction
|
Plant hormone signal transduction
|
ko04075
|
44
|
136
|
0.008
|
Organismal Systems
|
Environmental adaptation
|
Plant-pathogen interaction
|
ko04626
|
38
|
122
|
0.023
|
Differential expression analysis
A total of 7540 DEGs were identified. 4326 of these DEGs were up-regulated in diseased group, and 3214 were down-regulated (Fig.3a, b). In addition, 40440 Unigenes did not demonstrate significantly differential expression. In other words, DEGs between healthy and diseased samples accounted for 15.71% of all Unigenes.
Transcription factor prediction
By the standard of FDR<0.01 and FC>2, 1295 DEGs were identified as transcription factors with transcription factor prediction tool (Fig.4). Here, transcription factor family covers transcription factor (TF), transcription regulator (TR), protein kinases (PK). It could be clear to see that many DEGs were the members of transcription factor families MYB, ERF, C2H2, NAC, bHLH, C3H, WRKY, bZIP, GRAS, PHD, SNF2, SET. Coincidently, most of those transcription factor families have been proved that their expression or regulation can directly or indirectly affect plant disease resistance[39-58]. Exceptionally, present reports about C3H are mainly related to cold resistance, rather than disease resistance[59, 60].
KEGG pathways analysis
So far, it has been proved that plant disease resistance is relative to plant-pathogen interaction, plant hormone signal transduction, and other pathways about certain secondary metabolite biosynthesis or metabolism[56, 61-64]. Consistently, we got similar results in this study (Fig.5-7, Table 1).
In plant-pathogen interaction (Fig.5), genes except WRKY1/2 were all up-regulated. They were CDPK (calcium-dependent protein kinase), Rboh (respiratory burst oxidase homolog), CNGC (cyclic nucleotide gated channel), calcium-binding protein CML (calmodulin-like protein), LRR (leucine-rich repeat) receptor-like serine/threonine-protein kinase FLS2, MEKK1 (mitogen-activated protein kinase kinase kinase 1), MKK4/5 (mitogen-activated protein kinase kinase 4/5), WRKY transcription factor 33, WRKY transcription factor 22, RIN4 (RPM1-interacting protein 4), serine/threonine-protein kinase PBS 1, molecular chaperone HtpG. Biological processes these up-regulated genes principally involved were hypersensitive response (HR), cell wall reinforcement, defense-related gene induction, phytoalexin accumulation and miRNA production. Some of these genes were involved in PAMP-triggered immunity. Only WRKY transcription factor 2 displayed down-regulated expression, and it was connected with HR, defense-related gene induction and programmed cell death.
Further to analyze the map of plant hormone signal transduction (Fig.6), we learned that GH3 (auxin responsive glycoside hydrolase 3 gene family), AHP (histidine-containing phosphotransfer protein), ARR-B (two-component response regulator ARR-B family), PIF4 (phytochrome-interacting factor 4), ERF1 (ethylene-responsive transcription factor 1), JAZ (jasmonate ZIM domain-containing protein) were up-regulated. In the same pathway, AUX1 (auxin influx carrier), ARF (auxin response factor), CRE1 (cytokinin receptor enzyme), DELLA protein, PP2C (protein phosphatase 2C), EIN2 (ethylene-insensitive protein 2), BZR1/2 (brassinosteroid resistant 1/2), JAR1 (jasmonic acid-amino synthetase), COI1 (coronatine-insensitive protein 1), transcription factor TGA showed down-regulated. As it described, transcription factor TGA is connected with disease resistance. DEGs in this metabolic pathway involved many biological processes, such as cell enlargement, plant growth, cell division, shoot initiation, stem growth, stomatal closure, seed dormancy, fruit ripening, senescence, monoterpenoid biosynthesis, indole alkaloid biosynthesis, cell elongation, of course, disease resistance as well. Interestingly, above biological processes usually accompanied by phosphorylation (+p), dephosphorylation (-p), ubiquitination (+u). Phosphorylation and ubiquitination are common post-translational modification of proteins. They play an important role in pattern-triggered immunity (PTI), and simultaneously be necessary to receptor complex activation signals and cell homeostasis. Noteworthily, phytohormone played a vital role in this pathway. specifically, they included jasmonic acid (JA), salicylic acid (SA), ethylene (ET), brassinosteroid (BR), auxin, cytokinine, gibberellin, abscisic acid.
Furthermore, numerous DEGs regulating secondary metabolites biosynthesis, like ubiquinone, brassinosteroid, diterpenoid, phenylpropanoid, flavone and flavonol, revealed active state. In addition, genes related metabolic process of secondary metabolites also revealed significant differential expression. In starch and sucrose metabolism pathway, DEGs related to fructose and glucose synthesis is up-regulated, while DEGs linked to starch and glycogen production showed down-regulated expression.
In particular, brassinosteroid is one of crucial phytohormone closely related to plant growth and stress response. In brassinosteroid biosynthesis map, CYP90D2 (steroid 3-oxidase) showed up-regulated expression; CYP90A1 (cytochrome P450 family 90 subfamily A polypeptide 1) displayed down-regulated expression; CYP734A1/BAS1 (PHYB activation tagged suppressor 1) was mix-regulated, with two genes up-regulated and one gene down-regulated (Fig.7).
In summary, fungal disease response is a complex process involving multiple biological processes. And meanwhile, there is more than one gene participated in this process. Interestingly, one gene was not appeared in single pathway. That is to say, one gene may perform more than one function simultaneously. Anyway, these significantly enriched pathways would well reveal the underlying response mechanism of fungal disease in G. elata Bl. f. glauca.
Potential immune response mechanism of fungal disease in G. elata Bl. f. glauca
Actually, in this study, many genes related to stress response and disease resistance demonstrated high expression and significant difference. Moreover, They were members of certain transcription factor families, like WRKY, GH3, JAZ, CML, ERF, TGA. Furthermore, these genes were closely connected with derivatives of jasmonic acid, salicylic acid, brassinosteroid, ethylene and auxin. By BLAST (blast.ncbi.nlm.nih.gov/Blast.cgi), it is revealed that amino acid sequences of four JAZ genes in G. elata family were highly similar to certain sequences in Dendrobium catenatum, Phalaenopsis equestris, Apostasia shenzhenica (Figure 8). Noteworthily, they were all belong to TIFY10 family.
Finally, we found 10 candidate genes responding to fungal disease in G. elata Bl. f. glauca (Fig.8; Table 2). Seven of them participated in plant hormone signal transduction (ko04075) pathway, and three were part of plant-pathogen interaction (ko04626) pathway.
Table 2 Information of disease resistance genes. ko04626: plant-pathogen interaction; ko04075: plant hormone signal transduction. K13425: WRKY22; K14487: GH3; K13464: JAZ; K13448: CML; K14516: ERF1; K13424: WRKY33; K14431: TGA.
gene ID
|
FDR
|
log2FC
|
regulated
|
pathway
|
KEGG entry
|
nr annotation
|
c65017.graph_c0
|
8.65E-170
|
10.60551
|
up
|
ko04626
|
K13425
|
probable WRKY transcription factor 27, partial [Phalaenopsis equestris]
|
c32310.graph_c0
|
7.67E-97
|
9.531941
|
up
|
ko04075
|
K14487
|
probable indole-3-acetic acid-amido synthetase GH3.1 [Phalaenopsis equestris]
|
c75818.graph_c1
|
5.48E-19
|
5.547811
|
up
|
ko04075
|
K13464
|
protein TIFY 10a-like [Dendrobium catenatum]
|
c76234.graph_c1
|
1.23E-15
|
5.462372
|
up
|
ko04075
|
K13464
|
protein TIFY 10c-like [Dendrobium catenatum]
|
c60520.graph_c0
|
1.11E-77
|
4.956140
|
up
|
ko04075
|
K13464
|
protein TIFY 10a-like [Dendrobium catenatum]
|
c75190.graph_c0
|
3.26E-110
|
4.430928
|
up
|
ko04626
|
K13448
|
probable calcium-binding protein CML18 [Phalaenopsis equestris]
|
c78203.graph_c0
|
7.49E-63
|
4.276427
|
up
|
ko04075
|
K14516
|
ethylene-responsive transcription factor 1B-like [Dendrobium catenatum]
|
c78388.graph_c1
|
1.21E-46
|
3.499401
|
up
|
ko04075
|
K13464
|
protein TIFY 10a-like [Dendrobium catenatum]
|
c71906.graph_c0
|
6.55E-51
|
2.652510
|
up
|
ko04626
|
K13424
|
WRKY transcription factor WRKY24-like isoform X1 [Dendrobium catenatum]
|
c74033.graph_c1
|
8.03E-27
|
-2.15672
|
down
|
ko04075
|
K14431
|
transcription factor TGA1-like [Dendrobium catenatum]
|