Statistical analysis of RNA-seq results
A. alternata morphology,symptom changes in inoculated chrysanthemum leaves, and the dual RNA-seq analysis process are shown in Figure 1. Three samples sets, each with three biological replicates, were subjected to dual RNA-seq at each time point, and 27 cDNA libraries were generated: CK1h_1, CK1h_2, CK1h_3, CK12h_1, CK12h_2, CK12h_3, CK24h_1, CK24h_2, CK24h_3, Aa1h_1, Aa1h_2, Aa1h_3, Aa12h_1, Aa12h_2, Aa12h_3, Aa24h_1, Aa24h_2, Aa24h_3, In1h_1, In1h_2, In1h_3, In12h_1, In12h_2, In12h_3, In24h_1, In24h_2, In24h_3 (CK: control treatment, only chrysanthemum leaves; In: inoculation treatment, contain chrysanthemum leaves and A. alternata mycelium; Aa: only A. alternata mycelium). Table S2 shows the summary statistics of original reads and filtered clean reads of three replicates at each time point, for chrysanthemum. Inoculated and control samples generated on a average 41.20 Mb and 108.76 Mb clean reads, respectively, all with 100% read ratio. Furthermore, at each tome point, inoculated and control samples contained average 4.12 Gb clean bases and 10.88 Gb of clean bases, respectively. Table S3 lists the summary statistics of original reads and filtered clean reads obtained from three replicates at each time point, for A. alternata. The average clean reads of the inoculated and control samples generated on average 108.93 Mb and 108.64 Mb clean reads, respectively, with a read ratio ≥ 92.54%. Moreover, at each time point, inoculated and control samples contained an average of 10.89 Gb and 10.86 Gb of clean bases, respectively.
The clustered quality indicators of chrysanthemum are shown in Table S4. Infected and control chrysanthemum samples contained an average of 35,843 and 54,560 unigenes, respectively. The total length of chrysanthemum library transcripts was ≥ 14,571,366, the average length of the library was ≥ 642, N50, N70, and N90 ≥ 865, 556, 296, respectively. The GC ratio was ≥ 40.45%. Comparison of all unigenes to the seven major functional databases for annotation, generated the following numerical data: 89,889 (NR: 72.62%), 55,679 (NT: 44.98%), 61,156 (SwissProt: 49.41%), 64,694 (KOG: 52.26%), 64,705 ( KEGG: 52.27%), 68,727 (GO: 55.52%), and 60,671 (Pfam: 49.01%) (Table S5). The average total mapping percentage of A. alternata at each time point was higher than 56.9%, and that of the control group was higher than 86.53% (Table S6).
Identification of DEGs
Comparison of gene expression between the ‘In’ and ‘CK’ sample series detected 27,029 DEGs (21,216 up-regulated and 5,813 down-regulated) for In1h vs. CK1h, 76,932 DEGs (18,446 up-regulated and 58,486 down-regulated) for In12h vs. CK12h, and 49,571 DEGs (29,642 up-regulated and 19,929 down-regulated) for In24h vs. CK24h (Figure 2A). Illustration of these results as a Venn diagram clearly showed that both unique and shared DEGs were identified between, and among times points (Figure 2B). For example, 18,318, 20,696, and 37,618 shared DEGs were detected in the 1 HPI vs. 12 HPI, 1 HPI vs. 24 HPI, and 12 HPI vs. 24 HPI comparisons, respectively, while 15,960 DEGs were found in the 1 HPI vs. 12 HPI vs. 24 HPI comparison (Figure 2B). These results suggested that, as pathogen infection progressed, an increasing number of genes became involved in defense responses.
The degree of GO term enrichment was similar for the three inoculation time points, and DEGs were divided into 54 functional categories according to biological processes (25), cellular components (16), and molecular functions (13). The most significantly enriched biological processes were “regulation of transcription, DNA-templated”, “carbohydrate metabolic process”, and “translation”; the most significantly enriched cellular components were “cytoplasm”, “ribosome”, and “chloroplast”, while “protein serine/threonine kinase activity”, “nucleic acid binding” and “oxidoreductase activity” occupied the important positions in molecular functions (Figure 2D). Moreover, a total of 30 KEGG pathways were significantly enriched at 1, 12, and 24 HPI, each with a varying number of DEGs (Table S7). Maps with the highest DEG representation were those for ‘plant-pathogen interactions’ (ko 04626), followed by those for ‘plant hormone signal transduction’ (ko 04075), ‘MAPK signaling pathway-plant’ (ko 04016), ‘carbon metabolism’ (ko 01200) ‘protein processing in endoplasmic reticulum’ (ko 04141) and ‘biosynthesis of amino acids’ (ko01230). The above results indicated that chrysanthemum infected with A. alternata involved a series of defense strategies interacting with multiple pathways to jointly regulate and respond to pathogenic stress. These strategies dominated at different infection stages.
Gene expression comparison between the ‘In’ and ‘Aa’ sample series found 4,027 DEGs (2,729 up-regulated and 1,298 down-regulated) for In1h vs. Aa1h, 5364 DEGs (3697 up-regulated and 1667 down-regulated) for In12h vs. Aa12h, and 5,541 DEGs (3572 up-regulated and 1969 down-regulated) for In24h vs. Aa24h (Figure 2A). Moreover, analysis of Venn diagram showed 2,881, 3,066, 4,176, and 2,584 shared DEGs in the 1 HPI vs. 12 HPI, 1 HPI vs. 24 HPI, 12 HPI vs. 24 HPI, and 1 HPI vs. 12 HPI vs. 24 HPI comparison (Figure 2C).
The degree of GO term enrichment was similar among the three stages of A. alternata infection and DEGs were divided into 38 functional categories, according to biological processes (16), cellular components (12), and molecular functions (10). “Metabolic process”, “organic substance metabolic process” and “cellular process” were the most significantly enriched biological processes; the most significantly enriched cellular components were detected in “cell”, “cell part”, and “intracellular”; while “hydrolase activity”, “organic cyclic compound binding” and “heterocyclic compound binding” occupied the important positions in molecular functions (Figure 2E). Moreover, a total of 20 KEGG pathways were significantly enriched at the three stages, but with varying numbers of DEGs (Table S8). Maps with the highest DEGs representation were for ‘biosynthesis of antibiotics’ (ko 01130), followed by ‘MAPK signaling pathway-yeast’ (ko 04011), ‘amino sugar and nucleotide sugar metabolism’ (ko 00520) and ‘glycine, serine and threonine metabolism’ (ko 00260). The above results indicated that A. alternata induced a variety of metabolic activities during chrysanthemum infection, that generated energy and toxic metabolites to attack host cells. These metabolic processes played a key role in the interaction between chrysanthemum and A. alternata.
DEGs involved in phytohormone signaling
Phytohormones, such as salicylic acid (SA), ethylene (ET), jasmonic acid (JA), brassionosteroid (BR), auxin (AUX), and abscisic acid (ABA), are widely involved and play critical regulatory roles in plant-pathogen interactions [16]. The related DEGs of several hormone signaling pathways in infected chrysanthemum leaves were analyzed. Several DEGs involved in SA biosynthesis and signaling were differentially expressed, e.g., three DEGs of NPR1 (non-expressors of disease-related genes) and TGA (TGACG motif binding factor) were down-regulated at 1HPI but significantly up-regulated at 24 HPI; two DEGs related to PR1 (proteins related to disease-course) were up-regulated at 1 HPI and one of them was up-regulated at 24 HPI. All DEGs from JAZ (jasmonate ZIM-domain) were up-regulated during the whole process, and those MYC2 homologous were significantly up-regulated at 24 HPI. Several genes known to be ET-responsive were up-regulated, including EBF 1/2 (F-box protein 1/2 bound to EIN3) at 24 HPI, EIN3 (ethylene-insensitive 3) at 1 HPI and 24 HPI, and ERF1/2 (ethylene response factor 1/2) which exhibited change by a higher multiple. Most DEGs in AUX signaling, such as AUX / IAA (auxin / indole-3-acetic acid), SAUR (small auxin-up RNA) and auxin-responsive GH3 (Gretchen Hagen3 gene) also showed notably up-regulated expression. Previous studies have also shown that BRs comprise a unique class of growth-promoting steroid hormones, known to be key regulators of plant immunity [17]. DEGs encoding BR signaling cascades included BAK1 (BRI-related receptor kinase 1), BSK (brassinosteroid steroid signal transduction kinase), TCH4 (xyloglucan endotransglucosylase transglucosylase, also named Touch 4) and BZR1/2 (an anti-azole transcription factor). Except for DEGs encoding BZR1/2, that responded to A. alternata at 1 HPI, but were down-regulated by a high multiple (7) at 12HPI, the remaining DEGs belonging to the BR signaling cascades, expression level gradually increased at three infection stages. Finally, DGEs involved in the ABA signaling pathway, such as PYR/PYL (pyrabactin resistance1/PYR1-like), PP2C (protein phosphatase 2C), and SnRK2 (sucrose non-fermenting 1-related protein kinase2), were all up-regulated at 24 HPI; ABF (ABA binding factor) was up-regulated at 1 HPI and 24 HPI, but down-regulated at 12 HPI, similar to BZR1/2. The schematic diagram of the relevant hormone pathways is shown in Figure 3A.
DEGs involved in plant-fungal interaction
During biotic stress, chrysanthemum DEGs encoding CDPK (calcium-dependent protein kinase) and Rbohs (respiratory explosive oxidase homologs) were significantly up-regulated over time and were involved in hypersensitive reaction (hereafter named HR) and cell wall reinforcement. DEGs encoding PR-proteins, such as chitinase (CHI, PR3), were generally up-regulated, while the DEGs encoding b-1,3-glucanase (PR2) were specially down-regulated at 1 HPI, but significantly up-regulated at 12 HPI, with a high multiple (> 6) notably induced at 24 HPI. In addition, most DEGs encoding potential CNGCs (cyclic nucleotide gated channels) were up-regulated at 24 HPI, and those encoding CaM/CMLs (crassulacean acid metabolism/calmodulin-like protein) showed similar trends. Furthermore, several downstream defense-related PR proteins, such as PR9 (peroxidase), PR10 (ribonuclease), and PR14 (lipid-transfer protein), were induced and significantly up-regulated at all three time points (Figure 3B).
DEGs related to virulence in A. alternata
Alternaria spp. produce a variety of secondary metabolites during the pathogenic process, and more than 70 compounds with significant toxicity have been isolated [18], with important roles in fungal virulence. Most of these toxins are versatile compounds of polyketides and non-ribosomal peptides, which are usually generated by non-ribosomal peptide synthase (NRPS) and polyketide synthase (PKS), respectively [19]. We identified three NRPS genes and seven PKS genes in A. alternata, all showing up-regulated expression at the three inoculation stages. DEGs of NRPS (CC77DRAFT_1065195) presented a significantly higher multiple (> 6) at 24 HPI (Figure 4). A previous study have shown that pksJ and pksH were correlated with the production of alternariol (AOH) and alternariol-9-methyl ether (AME) [20]. We also identified the pksJ homolog (CC77DRAFT_1058721) and pksH homolog (CC77DRAFT_976935), both up-regulated during A. alternata infection (Figure 4).
Furthermore, fungal cell wall degrading enzymes (CAZymes) can promote degradation of the plant cell wall, penetration into the host tissue, and adhesion layer formation [21]. CAZymes consist of four functional classes: glycoside hydrolases (GHs), glycosyl transferases (GTs), polysaccharide lyases (PLs), and carbohydrate esterases (CEs), classified according to their catalytic modules or functional domains [21]. Expression levels were also investigated during the interaction between A. alternata and chrysanthemum. There were thirteen DEGs of GHs, with most of them significantly up-regulated at the three time points; two DEGs of GTs, one up-regulated and the other down-regulated. Only one DEG of PLs was found, and its expression showed an obviously upward trend. Similar to PLs, only one CE displayed higher expression (Figure 4).
Weighted gene co-expression network analysis (WGCNA)
Weighted gene co-expression network analysis (WGCNA) was carried out to identify genes related to phenotypes and investigated the co-expression networks to elucidate the interaction network between C. morifolium and A. alternata. Ultimately, 17 and 29 gene co-expression modules were discovered in C. morifolium and A. alternata, respectively, shown in Figure 5A and B.
Genes from the ‘Cm_brown’, ‘Cm_midnightblue’, ‘Cm_salmon’, ‘Cm_greenyellow’, ‘Cm_lightcyan’, ‘Aa_magenta’, ‘Aa_yellow’, ‘Aa_brown’, ‘Aa_skyblue’ and ‘Aa_black’ modules were highly correlated with the traits observed at the three infection stages (Figure 5A, B). KEGG annotation analyses were performed to further explore what pathways the genes from the modules above were involved in. Plant cell walls can act as a natural physical barrier against pathogens [21]. The cuticle is the first layer of the cell wall that prevents pathogens from invading the cells [22], and usually consists of a horny and a waxy protective film. Its biosynthesis involves several genes, including wax-ester synthase/diacylglycerol O-acyltransferase (WSD) [23], and fatty acid omega-hydroxy dehydrogenase (HTH) [24]. In the ‘Cm_turquoise’ module, WSD and HTH homologs in C. morifolium were significantly up-regulated, and several of these gene homologs (e.g., Unigene36512_All, Unigene36513_All, and CL1059.Contig1_All) displayed an expression fold-change > 10 (Figure S1). The above analysis showed that the cuticle played a positive role in the chrysanthemum defense against A. alternata. As reported, pectin lyase, pectate lyase, and xylanase can break down the pectin and xylans present in the plant cell wall. In the ‘Aa_black’ and ‘Aa_green’ modules, lyase homologs (e.g., CC77DRAFT_1043109, CC77DRAFT_1048882, and CC77DRAFT_167134) exhibited a high expression level in A. alternata during infection (Figure S1).
Highly correlated modules and key genes identification
The relationship between module and trait allowed us to evaluate the correlation coefficient between modules from C. morifolium and A. alternata. A network of C. morifolium and A. alternata modules were shown in Figure 6A, and highly correlated modules (r ≥ 0.8 and p-value < 0.05) were linked by a line (Figure 6A). In the early infection stage, three C. morifolium gene modules (‘Cm_brown’, ‘Cm_midnightblue’, and ‘Cm_salmon’,) and five A. alternata modules (‘Aa_magenta’, ‘Aa_yellow’, ‘Aa_brown’, ‘Aa_darkorange’, and ‘Aa_pink’) were highly correlated; in the middle infection stage, one C. morifolium gene modules (‘Cm_greenyellow’) and three A. alternata modules (‘Aa_royalblue’, ‘Aa_black’, and ‘Aa_green’) were highly correlated; in the late infection stage, two C. morifolium gene modules (‘Cm_turquoise’ and ‘Cm_lightcyan’) and four A. alternata modules (‘Aa_steelblue’, ‘Aa_grey60’, ‘Aa_black’, and ‘Aa_green’) were highly correlated (Figure 6B). Gene significance (value ≥ 0.8) and connectivity (top 20%) were used together to identify key genes in each of the modules above (Figure 6C, D).
Based on gene transcription levels, we performed correlation coefficient analyses between genes from highly correlated modules in C. morifolium and A. alternata to identify the interplay genes. From modules that highly correlated with the late infection stage, a number of genes were identified and a network of highly correlated genes (r ≥ 0.8 and p-value < 0.05) were linked by a line, as shown in Figure 7A. Most of these genes were up-regulated with the spread of A. alternata (Figure 7B, C, D, E).
Regulatory network of key genes annotated in PPI and PRG databases
The Plant Resistance Genes database (PRGdb; http://prgdb.org) is a bioinformatics platform for plant resistance gene analysis [25]. The pathogen–host interactions database (PHI-base: www.phi-base.org) contains molecular and biological information on genes which have been proven to affect the outcome of host-pathogen interactions [26]. Key genes from the highly correlated modules were examined for further analyses. Seventy-five key genes from A. alternata were annotated by PHI-base and twelve key genes were annotated by PRGdb. The regulatory network of these eighty-seven key genes were shown in Figure 8. Notably, two RGA1-like disease resistance protein homologs (CL2806.Contig1_All and CL2806.Contig4_All) were identified in the late infection stage. Two transcription factors (CL14283.Contig1_All and Unigene25854_All) were also identified and may play important roles in response to A. alternata infection in C. morifolium (Figure 8). Additionally, ACL2 homolog (CC77DRAFT_784023), ACL1 homolog (CC77DRAFT_986135), and BUF1 homolog (CC77DRAFT_528893) were identified, which were predicted to influence the virulence of A. alternata (Figure 8).
Validation of RNA-seq data by qRT-PCR
To confirm the reliability of the generated dual RNA-seq data, the expression of 12 DEGs were analyzed using qRT-PCR assays, of which six were derived from chrysanthemum, (CL11098.Contig2_All, CL1653.Contig1_All, CL5572.Contig1_All, Unigene47090_All, CL3907.Contig2_All and CL11265.Contig3_All); Figure 9A), and six were from A. alternata (CC77DRAFT_945175, CC77DRAFT_1044312, CC77DRAFT_1036704, CC77DRAFT_598231, CC77DRAFT_779096 and CC77DRAFT_950634; Figure 9B). qRT-PCR results and RNA-seq data showed similar up-regulation or down-regulation expression patterns. The correlation coefficients between qRT-PCR and RNA-seq of the 12 DEGs were all ≥ 0.85. Minor discrepancies regarding the expression levels might suggest a difference in sensitivity between the two methods. These results highlighted the reliability of the RNA-seq data.