Transcriptome analysis reveals pathways facilitating the growth of tobacco powdery mildew in Arabidopsis

Abstract

Powdery mildew (PM) fungi are biotrophic pathogens that rely on living hosts to survive and thrive. However, their colonization is restricted by host defenses at both the penetration and post-penetration stages. The tobacco PM strain Golovinomyces cichoracearum (Gc) SICAU1 has overcome penetration resistance of Arabidopsis but its growth is arrested by post-penetration resistance. While Gc SICAU1 only poorly grows in Arabidopsis Col-0 wild-type plants, it can sustainably grow for more than 20 days on the same infected leaves of the double mutant pad4–1 sid2–1 that is defective in both the synthesis and signaling of salicylic acid (SA). To understand the underlying molecular mechanisms, we conducted a comparative transcriptome analysis between Col-0 and pad4–1 sid2–1 in response to Gc SICAU1. We found that 4811 genes were differentially expressed more than four-fold between any two of the measured seven time points (0, 1, 3, 6, 8, 10 and 12 days post-inoculation). Gene expression pattern analysis suggests that differential expression of 348 genes and 190 genes may explain resistance in Col-0 and susceptibility in pad4–1 sid2–1, respectively. Gene Ontology (GO) analysis suggests that Gc SICAU1 might be arrested in Col-0 by both pattern-triggered immunity and SA-dependent defense. By contrast, its sustained growth in pad4–1 sid2–1 may be attributable to the activation of a detoxification pathway that is normally repressed by the SA-signaling pathway. Taken together, our results suggest that multiple distinct, yet interconnected pathways control the growth of tobacco powdery mildew in Arabidopsis.

Background

Powdery mildew (PM) diseases epidemically occur in nearly 10,000 plant species including many economically and agriculturally important crops (Kuhn et al. 2016). Successful colonization of a powdery mildew begins with germination of a conidium upon landing on the surface of a host plant, which is followed by the formation of an appressorium that further develops a penetration peg to break the host cell wall. At the tip of the penetration peg, it forms a feeding structure called the haustorium to steal photosynthates for its epiphytic growth (Kuhn et al. 2016). During the infection processes, the host plant mounts spatiotemporally distinct defenses to stop the invasion of the pathogen. The first line of defense is penetration resistance, which is essentially mounted by the rigid host cell wall as a physical barrier; attempted fungal penetration induces the formation of the papilla, i.e. cell wall apposition, enriched for callose at the penetration site (Collins et al. 2003; Assaad et al. 2004). Using Arabidopsis-powdery mildew interaction as a model phytopathosystem, it has been demonstrated that papilla formation may be part of pattern-triggered immunity (PTI), which is activated upon recognition of chitin by the pattern-recognition receptor CERK1 (CHITIN ELICITOR RECEPTOR-LIKE KINASE 1) together with LYK5 (LYSIN-MOTIF RECEPTOR LIKE KINASE 5) (Miya et al. 2007; Cao et al. 2014). Consistently, the cerk1 mutant is more susceptible to the adapted PM Golovinomyces cichoracearum than the wild type (Wan et al. 2008). Application of chitin results in callose deposition at the treated tissue (Underwood and Somerville 2013). Conversely, loss-of-function mutations in PENETRATION 1 (PEN1) lead to delayed formation of the papilla and increased penetration rate of the non-host powdery mildew pathogen Blumeria graminis f. sp. hordei in Arabidopsis (Collins et al. 2003; Assaad et al. 2004). Accumulation of callose as a major constituent of papillae requires the callose synthase PMR4/GSL5 (POWDERY MILDEW RESISTANCE 4 / GLUCAN SYNTHASE LIKE 5) (Nishimura et al. 2003). Although loss-of-function in PMR4 leads to enhanced resistance to adapted PM, which may be due to compensatory activation of SA-dependent defense, overexpression of PMR4 or the GTPase RabA4c enhances callose deposition to the papilla, leading to complete penetration resistance to PM (Ellinger et al. 2013; Eggert et al. 2014; Ellinger et al. 2014). Indeed, the formation of a callose encasement of the haustorium complex (EHC) was observed in some invaded epidermal cells of Col-0, an Arabidopsis accession susceptible to an adapted PM (Wang et al. 2009), but the frequency of EHC formation is much higher in plants challenged with a non-adapted PM (Wen et al. 2011), implying that adapted PM pathogen can suppress callose deposition and more effectively break penetration resistance to establish colonization. Not surprisingly, the frequency of EHC is significantly increased in plants expressing RESISTANCE TO POWDERY MILDEW8 genes (including RPW8.1 and RPW8.2). Moreover, the expression of RPW8.2 can be induced to mediate cell death and the RPW8.2 protein is specifically targeted to the extra-haustorial membrane (EHM) encasing the haustorium (Wang et al. 2009), which nicely explains why RPW8 mediates broad-spectrum resistance to PM. Apparently, RPW8.2 activates EHM-based, post-penetration resistance including enhanced EHC formation against adapted PM (Kuhn et al. 2016).

Certain host genes are required for successful pathogenesis of PM and thus act as susceptibility factors. Loss-of-function in such genes can result in enhanced disease resistance to PM. For example, some plant mutants display resistance to adapted PM pathogens that were impaired for potential PM susceptibility. Such mutants include edr (enhanced disease resistance), pmr (powdery mildew resistance) and mlo (mildew locus O) mutants of Arabidopsis, barley or wheat (Buschges et al. 1997; Frye et al. 2001; Vogel et al. 2002; Nishimura et al. 2003; Tang et al. 2005; Consonni et al. 2006; Tang et al. 2006; Wang et al. 2014; Wu et al. 2015). Therefore, it is feasible to engineer disease resistance in crops by identifying and editing susceptibility factors.

The tobacco powdery mildew strain Gc SICAU1 is a poorly-adapted pathogen to Arabidopsis. Although it can finish its life cycle in the Arabidopsis wild type Col-0, its growth is arrested because of post-penetration defense (Zhang et al. 2015). Nevertheless, the double mutant pad4–1 sid2–1 can support its growth to more than 20 days which is unusually long because vigorous PM growth in a compatible host normally lasts for 12–15 days before the collapse of the infected leaves, indicating that the mutant not only lacks the post-penetration defense, but some pathways are activated to support the sustained growth of the PM strain (Zhang et al. 2015). Thus, PAD4 (PHYTOALEXIN DEFICIENT 4) and SID2 (SALICYLIC ACID INDUCTION DEFICIENT 2) are important for post-penetration resistance to tobacco powdery mildew in Arabidopsis. PAD4 interacts directly with EDS1 (ENHANCED DISEASE SUSCEPTIBILITY 1) to mediate basal defense and SA-signaling-dependent disease resistance (Zhou et al. 1998; Rusterucci et al. 2001). SID2 encodes the isochorismate synthase acting in SA biosynthesis (Nawrath and Metraux 1999). Thus, Arabidopsis-Gc SICAU1 interactions could be an ideal pathosystem to dissect the mechanism of sustained growth of PM in Arabidopsis.

In the present study, to understand why loss of PAD4 and SID2 results in unusually long sustained growth of tobacco powdery mildew in pad4–1 sid2–1, we performed a comparative transcriptome analysis between Col-0 and pad4–1 sid2–1 upon Gc SICAU1 inoculation. Our data confirmed that basal defense and SA-signaling-associated pathways contribute to post-penetration defense against tobacco powdery mildew. Moreover, a detoxification pathway contributes to the sustained growth of tobacco PM in Arabidopsis.

Results

The number of genes responsive to tobacco powdery mildew is quite different between Col-0 and pad4–1 sid2–1

Previously, we isolated a tobacco powdery mildew strain, Gc SICAU1, which poorly sporulated in Arabidopsis ecotype Col-0 (Zhang et al. 2015). However, it can grow in the double mutant pad4–1 sid2–1 abundantly as in its host tobacco plants (Fig. 1). To understand the underlying mechanism, we performed a time course transcriptome analysis by collecting samples at 0, 1, 3, 6, 8, 10, and 12 days post inoculation (dpi) with two biological replicates for RNA-seq. Data from the two biological replicates were highly consistent as demonstrated by the correlation coefficient matrix between samples (Additional file 1: Figure S1a). Clean reads were mapped to TAIR10 reference genome. Both unique mapped rate and total mapped rate exhibited the following features: in the samples collected from before and at 3 dpi, the unique mapped rates for both Col-0 and pad4–1 sid2–1 were quite similar, being of about ~ 78% and the total mapped rates of about ~ 93% (Additional file 2: Table S1 and Additional file 1: Figure S1b, c). However, at and after 6 dpi in the mutant, the unique mapped rates and total mapped rates reduced over the time course and were lower than that in Col-0 (Additional file 2: Table S1 and Additional file 1: Figure S1b, c). The reduced percentage of mapped reads in pad4–1 sid2–1 could be due to the increase of the pathogen at the later time points.

Fig. 1

Powdery mildew disease phenotypes in tobacco and Arabidopsis. Powdery mildew disease phenotypes in a representative tobacco plant (a), leaf (b), indicated Arabidopsis plants (c) and leaves (d) at 12 dpi of the tobacco powdery mildew isolate Golovinomyces cichoracearum SICAU1

To identify the genes showing the most differential response to Gc SICAU1, we selected the differentially expressed genes (DEGs) at variations of absolute Log₂ fold change ≥2 and a false discovery rate ≤ 0.01 between any two time points. Subsequently, there were 4811 DEGs identified (Additional file 3: Table S2). Then, the DEGs were further screened by comparison between different time points and 0 dpi in Col-0 and in pad4–1 sid2–1, respectively. Compared to 0 dpi, the number of up-regulated genes increased from 255 genes at 1 dpi to 1101 genes at 12 dpi in Col-0, and from 176 genes at 1 dpi to 741 genes at 12 dpi in pad4–1 sid2–1 (Fig. 2a, b). Conversely, the number of down-regulated genes reduced from 494 genes at 1 dpi to 344 genes at 8 dpi, and then increased to 650 genes at 12 dpi in Col-0; whereas, the number of down-regulated genes reduced from 592 genes at 1 dpi to 433 genes at 8 dpi, and then increased to 630 genes at 12 dpi in pad4–1 sid2–1 (Fig. 2a, b). Obviously, there were more up-regulated genes in Col-0 than in pad4–1 sid2–1 at 1, 3, 10 and 12 dpi, but more up-regulated genes in pad4–1 sid2–1 than in Col-0 at 6 and 8 dpi; whereas, there were more down-regulated genes in pad4–1 sid2–1 than in Col-0 at 1–10 dpi, but similar number at 12 dpi (Fig. 2a, b).

Fig. 2

Number of differentially expressed genes over the time course of infection and between Col-0 and pad4–1 sid2–1. a Number of up-regulated and down-regulated genes in Col-0 at the indicated time points from 1 dpi (T1) to 12 dpi (T12) in comparison with un-inoculated plants (T0). b Number of up-regulated and down-regulated genes in pad4–1sid2–1 at the indicated time points (from T1 to T12) in comparison with un-inoculated plants (T0). c Number of up-regulated and down-regulated genes in pad4–1sid2–1 at the indicated dpi in comparison with Col-0 at the same time points. d Venn graph shows the gene numbers up-regulated in Col-0 (compared to 0 dpi, designated class A, blue cycle), down-regulated in pad4–1sid2–1 (compared to Col-0, designated class B, yellow cycle) and down-regulated in pad4–1sid2–1 (compared to 0 dpi, designated class C, green cycle). e Venn graph shows the gene numbers down-regulated in Col-0 (compared to 0 dpi, designated class a, blue cycle), up-regulated in pad4–1sid2–1 (compared to Col-0, designated class b, yellow cycle) and up-regulated in pad4–1 sid2–1 (compared to 0 dpi, designated class c, green cycle)

Because the DEGs obtained during the time course could be due to either the aging of the plants or the infection of the pathogen, we also compared the number of DEGs between Col-0 and pad4–1 sid2–1 at each time point. At 0 dpi, there were 26 genes expressed higher and 62 genes expressed lower in pad4–1 sid2–1 than in Col-0 (Fig. 2c). Over the time course of infection, the number of higher- and lower-expressed genes was increased, particularly, a sudden increase occurred to the number of higher-expressed genes in pad4–1 sid2–1 at 6 dpi and to the number of lower-expressed genes at 3 dpi (Fig. 2c). The higher-expressed genes in pad4–1 sid2–1 could contribute to the sustained growth of Gc SICAU1, whereas, the lower-expressed genes in pad4–1 sid2–1 could reflect the lack of post-penetration defense in pad4–1 sid2–1 and thus contribute to post-penetration defense in Col-0. Therefore, these data indicate that 3 dpi may be critical for transcriptome-reprogramming to defend against the pathogen and 6 dpi may be critical to establish the sustained growth of the pathogen.

The number of genes contributes to post-penetration defense in Col-0 and susceptibility in pad4–1 sid2–1

To identify the genes most likely associated with post-penetration resistance in Col-0 and susceptibility in pad4–1 sid2–1, we further performed comparison between the DEGs by making two Venn diagrams. One Venn diagram was made by comparing the genes exhibiting up-regulation in Col-0 over the time course versus 0 dpi (Class A) with the genes showing down-regulation in pad4–1 sid2–1 versus Col-0 at each time point (Class B), and the genes exhibiting down-regulation in pad4–1 sid2–1 over the time course versus 0 dpi (Class C) (Fig. 2d). Another diagram was made by comparing the genes down-regulated in Col-0 over the time course versus 0 dpi (Class a) with the genes up-regulated in pad4–1sid2–1 at each time point versus Col-0 (Class b), and the genes up-regulated in pad4–1 sid2–1 over the time course versus 0 dpi (Class c) (Fig. 2e). The genes consensually up-regulated in Col-0 and down-regulated in pad4–1 sid2–1 could be those contributing to post-penetration defense, whereas, the genes consensually upregulated in pad4–1 sid2–1 could contribute to susceptibility in the mutant. The data showed that there were 16 genes consensually detected in the Classes A, B and C, whereas,256 genes in the Classes A and B, 30 genes in the Classes A and C, and 46 genes in the Classes B and C (Fig. 2d and Additional file 4: Table S3). These genes could be the ones contributing to the defense in Col-0. On the other hand, there were 142 genes consensually detected in the Classes b and c, whereas, 35 genes in the Classes a and b, 7 genes in the Classes a, b and c, and 6 genes in the Classes a and c (Fig. 2e and Additional file 5: Table S4). These genes could be the ones contributing to the sustained growth of the pathogen in pad4–1 sid2–1.

Taken together, there are 348 genes that might contribute to the activation of post-penetration defense against the tobacco PM strain Gc SICAU1 in Col-0. Whereas, 190 genes may contribute to the susceptibility in pad4–1 sid2–1.

Categories of biological processes from GO analysis on DEGs contribute to post-penetration defense in Col-0 and susceptibility in pad4–1 sid2–1

In order to identify the biological functions contributing the most to post-penetration defense in Col-0 and susceptibility in pad4–1 sid2–1, we performed Gene Ontology (GO) analyses on the DEGs over the time course of infection of Gc SICAU1. Our results showed that the genes expressed higher in Col-0 than in pad4–1 sid2–1 were mainly enriched in 25 biological processes (Table 1). Whereas, the genes expressed higher in pad4–1 sid2–1 than in Col-0 were enriched in 20 biological processes (Table 1).

Table 1 GO enrichment of genes expressed higher in Col-0 or in pad4–1 sid2–1 over the time course

The 25 biological processes that may contribute to post-penetration defense in Col-0 can be further classified into three classes. Class I contained three categories of GO biological process terms, including receptor protein kinase signaling pathway, MAPK cascade and protein kinase activity, which mainly belong to PTI (Table 1 and Additional file 6: Table S5). In this class, we detected high induction of some receptor genes such as LYK5 (AT2G33580), WAK1 (AT1G21250), RLK5 (AT4G28490), RLP23 (AT2G32680), together with 16 cysteine-rich receptor-like protein kinases (CRKs), 11 receptor-like proteins (RLPs), three receptor kinases (RKs), two receptor-like kinases (RLKs) and one L-type Lectin Receptor Kinase (LECRK) (Additional file 6: Table S5). Further screening of DEGs identified 26 CRKs whose expression changed more than 4-fold between any 2 time points (Additional file 7: Table S8). Class II included 11 GO biological process terms (Table 1). In this class, four terms related to “SA signaling pathway”, including salicylic acid biosynthetic process, systemic acquired resistance, SA-mediated signaling pathway and response to SA, and there were 46 genes in this category with higher expression in Col-0 than in pad4–1 sid2–1 (Additional file 6: Table S6). Six genes from SA signaling pathway was also enriched in the JA-signaling pathway (Additional file 6: Table S6). In addition, there were eight more genes enriched in JA-signaling pathway (Additional file 6: Table S6). We used the category “defense response” to include all genes enriched in the biologic processes of the GO terms defense response, defense response to fungus and defense response in incompatible interaction. There were 30 genes in this category, of which 14 genes expressed higher in Col-0 than in pad4–1 sid2–1 at all tested time points and nine genes were at the highest expression at 12 dpi (Additional file 6: Table S6). Class III contained the rest of the 11 biological processes (Table1). Taken together, genes involved in PTI- and SA-associated pathways were expressed higher in Col-0 than in pad4–1 sid2–1.

The 20 biological processes that may contribute to the sustained growth of the tobacco powdery mildew in pad4–1 sid2–1 were further classified into “detoxification” and “secondary metabolism” (Table 1 and Additional file 6: Table S7). The “detoxification” class contained genes associated with the detoxification pathway, including 13 genes in the term of toxin catabolic process/response to toxic substance, of which four encode GSTUs (Glutathione S-Transferase), two encode UGTs (UDP-Glucose Transferase) and one encodes CYP (Cytochrome P450 monooxygenase); seven genes in the term proline transport/amino acid transport/glutathione transferase activity and nine genes involved in flavonoid biosynthetic process/flavonoid glucuronidation/flavonol biosynthetic process (Additional file 6: Table S7). Totally, seven UGT genes were differentially induced over the time course of infection of Gc SICAU1 with expression higher in pad4–1 sid2–1 than in Col-0 (Additional file 6: Table S7). Intriguingly, while UGT76E11 (AT3G46670) was expressed higher in pad4–1 sid2–1 than in Col-0 over the time course of infection, the other six UGT genes, i.e. UGT73B3 (AT4G34131), UGT73B4 (AT2G15490), UGT73B5 (AT2G15480), UGT73C1 (AT2G36750), UGT74E2 (AT1G05680) and UGT75B1 (AT1G05560), were expressed higher in Col-0 than in pad4–1 sid2–1 at one to three time points before 3 dpi, but higher in pad4–1 sid2–1 than in Col-0 at 6–12 dpi (Additional file 6: Table S7). In addition, we also detected the other 27 UGTs that were differentially expressed in Col-0 and pad4–1 sid2–1 upon infection of Gc SICAU1, most of them expressed higher in pad4–1 sid2–1 than in Col-0 (Additional file 8: Table S9). There were seven genes encoding apoplast proteins, eleven genes responsive to chitin/chitin binding or involved in chitin catabolic process, and two genes functioning in indole glucosinolate metabolic process (Additional file 6: Table S7). These genes provide the opportunity to further dissection of the mechanism in plant-powdery mildew interactions.

Expression patterns over the time course of infection classify genes contributing to post-penetration defense in Col-0 and susceptibility in pad4–1 sid2–1

To be more accurate in demonstrating the expression patterns of the genes involved in post-penetration defense in Col-0 and susceptibility in pad4–1 sid2–1, we performed analysis on the expression pattern of genes for the consensually up-regulated or down-regulated genes in Col-0 and pad4–1 sid2–1. First, we performed clustering analysis on all the DEGs that were classified into 119 clusters (Additional file 9: Figure S2). Then we further analyzed the clusters containing the consensually up-regulated or down-regulated genes in Col-0 and pad4–1 sid2–1. Subsequently, the clusters contributing to resistance were classified into four groups, i.e. Resistance group A (R_A), R_B, R_C and R_D, each represented by 8 clusters (Fig. 3). R_A included clusters 13, 48, 53, 57, 73, 75, 84 and 110 with an obviously up-regulated expression in Col-0 but not much change in pad4–1 sid2–1 over the time course (Fig. 3a). R_B included clusters 17, 30, 70, 72, 74, 82, 85 and 104 with a pattern of up-regulation over the time course but higher in Col-0 than in pad4–1 sid2–1 (Fig. 3b). R_C included clusters 7, 8, 10, 28, 4, 83, 96 and 98 with patterns of slightly up-regulation over the time course in both Col-0 and pad4–1 sid2–1 with expression higher in Col-0 than in pad4–1 sid2–1 (Fig. 3c). R_D included clusters 16, 39, 60, 61, 69, 86, 113 and 117 with patterns that the expression were decreased over the time course of infection, but the expression were higher in Col-0 than in pad4–1 sid2–1 (Fig. 3d).

Fig. 3

Four types of clustered DEGs involved in post-penetration defense. a Resistance type A clusters with DEGs up-regulated in Col-0 but no obvious change in pad4–1 sid2–1. b Resistance type B clusters with DEGs up-regulated in Col-0 but delayed up-regulation in pad4–1 sid2–1. c Resistance type C clusters with DEGs up-regulated in both Col-0 and pad4–1sid2–1 but the up-regulation in pad4–1sid2–1 are lower than that in Col-0. d Resistance type D clusters with DEGs expressed highly in Col-0 and down-regulated in pad4–1sid2–1

Meanwhile, we also classified the DEGs up-regulated in pad4–1 sid2–1 into five susceptibility groups. Susceptibility group A (S_A) was represented by clusters 5, 47, 92, 97, 102, 106 and 107, whose expression was obviously up-regulated in pad4–1 sid2–1, particularly at and after 6 dpi of Gc SICAU1 but not much change in Col-0 (Fig. 4a). S_B was represented by clusters 9, 22, 31, 34, 88, 91, 108 and 114 with expression increased over the time course of infection but higher in pad4–1 sid2–1 than in Col-0 (Fig. 4b). S_C was represented by clusters 21, 35, 37 and 101 with patterns of expression increased over the time course but higher in pad4–1 sid2–1 than in Col-0 at and after 6 dpi of Gc SICAU1 (Fig. 4c). S_D was represented by clusters 52, 58, 93 and 111 with patterns of expression decreased upon inoculation of Gc SICAU1, but the expression was higher in pad4–1 sid2–1 than in Col-0 (Fig. 4d). S_E was represented by clusters 1, 6, 42 and 44 with patterns of expression increased at 1–6 dpi and then dropped back at later time points but the expression was higher in pad4–1 sid2–1 than in Col-0 (Fig. 4e).

Fig. 4

Five types of clustered DEGs involved in susceptibility to powdery mildew. a Susceptibility type A clusters of DEGs dramatically up-regulated in pad4–1 sid2–1 particularly at and after 6 dpi but no obvious change or only slightly up-regulated in Col-0. b Susceptibility type B clusters of DEGs up-regulated in both Col-0 and pad4–1 sid2–1 but the up-regulation in pad4–1sid2–1 are higher than that in Col-0. c Susceptibility type C clusters of DEGs up-regulated in both Col-0 and pad4–1 sid2–1 but dramatically up-regulated in pad4–1sid2–1 particularly at and after 6 dpi. d Susceptibility type D clusters of DEGs expressed higher in pad4–1sid2–1 than in Col-0 but down-regulated upon powdery mildew infection. e Susceptibility type E clusters of DEGs up-regulated in both Col-0 and pad4–1 sid2–1 at 1~6 dpi and then dropped back at 8~12 dpi with expression higher in pad4–1sid2–1 than in Col-0

Reverse-transcription quantitative PCR confirmed the RNA-seq data

To validate the RNA-seq data, we conducted an independent time course experiment and collected samples for reverse-transcription quantitative PCR (RT-qPCR) on 19 genes selected from some resistance and susceptibility groups. Four genes from R_A, including PR1 (AT2G14610), BGL2 (AT3G57260), PR5 (AT1G75040) and DLO1 (AT4G10500), were highly up-regulated in Col-0 upon infection of Gc SICAU1 but only slightly up-regulated at later time points in pad4–1 sid2–1 (Fig. 5a). Four genes from R_B, including EFR72 (AT3G16770), PR4 (AT3G04720), NATA1 (AT2G39030) and ST2A (AT5G07010), were up-regulated in both Col-0 and pad4–1 sid2–1 upon infection of Gc SICAU1 with expression higher in Col-0 than in pad4–1 sid2–1 at most of the tested time points (Fig. 5b). Five genes from R_C were examined by RT-qPCR, of which four genes, including HCHIB (AT3G12500), ORA59 (AT1G06160), SEN4 (AT4G30270) and NAC6 (AT5G39610), were up-regulated in both Col-0 and pad4–1 sid2–1 upon infection of Gc SICAU1 with expression higher in Col-0 than in pad4–1 sid2–1 in most of the tested time points (Fig. 5c). Nevertheless, PDF1.2 (AT5G44420) from R_C was up-regulated higher in Col-0 than in pad4–1 sid2–1 at 1 and 3 dpi, but higher in pad4–1 sid2–1 than in Col-0 at and after 6 dpi (Fig. 5d). AT3G05730, a gene encoding a defensin-like family protein, was not classified into the R or S group, exhibiting an pattern of low expression at earlier time points but up-regulated in Col-0 and pad4–1 sid2–1 at and after 6 dpi with expression higher in Col-0 than in pad4–1 sid2–1 at 6, 8 and 10 dpi, higher in pad4–1 sid2–1 at 12 dpi (Fig. 5d). Two genes from S_C were examined of which WRKY6 (AT1G62300) was up-regulated in both Col-0 and pad4–1 sid2–1 with expression higher in pad4–1 sid2–1 than in Col-0 at 3 and 8 dpi, whereas, the expression of SEN1 (AT4G35770) was higher in pad4–1 sid2–1 at all examined time points except 10 dpi (Fig. 5e). Three genes from S_A were examined with expression higher in pad4–1 sid2–1 than in Col-0 (Fig. 5f). Apart from some genes at certain time points (11 of 133) that were not consistent with the RNA-seq data, which could be due to the systemic error occurred between different experiments, most of the time points (122 of 133) in the examined genes were consistent with the RNA-seq data (Additional file 3: Table S2). Therefore, our RNA-seq data were reliable and should be significant for future investigation.

Fig. 5

Validation of some DEGs by RT-qPCR. Reverse transcription quantitative PCR (RT-qPCR) analysis on the relative mRNA levels of the indicated four representative resistance type A genes (a), four representative resistance type B genes (b), four representative resistance type C genes (c), one representative resistance type C gene and one gene not in the classified group (d), two representative susceptibility type C genes (e), three representative susceptibility type A genes (f). Note in (a) that the mRNA levels were highly induced at and after 3 dpi in Col-0, but remained at low levels in pad4–1 sid2–1, and in (f) that the mRNA levels of CLE1 and At3g59930 were remarkably up-regulated at and after 6 dpi in pad4–1 sid2–1. Relative mRNA levels were normalized to 0 dpi of Col-0 that was arbitrarily set as 1. Values are means of three technique replications. Error bars indicate SD. Asterisks above the bars indicate significant differences between Col-0 and pad4–1 sid2–1 at P < 0.01 as determined by student’s t-test. All the experiments were repeated two times with similar results

Discussion

Multiple pattern-triggered immune pathways may constitute the basal defense preventing the growth of tobacco powdery mildew in Col-0

In this study, over a time course of infection by a poorly-adapted tobacco powdery mildew strain in Arabidopsis, we detected high induction of genes in receptor protein kinase signaling pathway, MAPK cascade and protein kinase activity, such as LYK5, WAK1, RLK5, RLP23 and 16 CRK genes, indicating that multiple PTI pathways may be involved in preventing the growth of tobacco powdery mildew in Col-0. First, chitin is a major component of fungal cell wall, presenting in conidia and at the tips of growing hyphae of powdery mildew (Ramonell et al. 2005). The chitin receptor gene LYK5 encodes a LysM-containing receptor-like kinase. Upon cognition of chitin, LYK5 forms a complex with CERK1 to activate down-stream innate immunity (Cao et al. 2014). Therefore, chitin-triggered immunity could be one contribution to the post-penetration defense against tobacco powdery mildew in Col-0 (Fig. 6). Second, Wall-associated Kinase 1 (WAK1) acts as a receptor of oligogalacturonides (OGs) derived from the plant cell wall after wound (Brutus et al. 2010; Gramegna et al. 2016). Here, we detected that WAK1 was expressed higher in Col-0 than in pad4–1 sid2–1 over the time course (Additional file 6: Table S5). This is consistent with its roles in response to wound because powdery mildew infection could cause wound on the host. Third, SOBIR1 forms a complex with the PRR protein RLP23 (AT2G32680), and recruits BAK1 into a tripartite complex upon recognizing necrosis- and ethylene-inducing peptide 1-like proteins (NLPs) that are produced by fungal, oomycete and bacterial pathogens to mediate immune activation (Albert et al. 2015). SOBIR1 acts in both penetration and post-penetration resistance to the non-host fungal pathogen Magnaporthe oryzae in Arabidopsis (Takahashi et al. 2016). Thus, RLP23-SOBIR1-BAK1-mediated PTI signaling pathway may play roles in post-penetration defense against tobacco powdery mildew (Fig. 6). Fourth, RLK5/HAE (AT4G28490) is a receptor for the INFLORESCENCE DEFICIENT IN ABSCISSION (IDA)-derived peptide (Wang et al. 2017). Although its roles in response to fungal pathogen are not reported, we detected that RLK5/HAE was expressed higher in Col-0 than in pad4–1 sid2–1 over the time course of infection (Additional file 6: Table S5), implying that RLK5-mediated PTI pathway may be involved in defense against tobacco powdery mildew (Fig. 6).

Fig. 6

Hypothetic working model for activation of post-penetration defense or sustained growth of tobacco powdery mildew in Arabidopsis

In addition, IOS1 (for IMPAIRED OOMYCETE SUSCEPTIBILITY 1, AT1G51800) plays critical roles in chitin-induced PTI via association with the PRR CERK1 and the co-PRR BAK1 (Yeh et al. 2016). SARD1 (At1g73805) and CBP60g (AT5G26920) are key regulators controlling SA synthesis via regulation of the expression of ICS1. In turn, the expression of SARD1 and CBP60g is regulated by PCRK1 (AT3G09830) and PCRK2 who also interact with FLS2 (Kong et al. 2016). SARD1 expression is also regulated by WRKY70 (Zhou et al. 2018). Here, we found that IOS1, SARD1 and WRKY70 were expressed higher in Col-0 than in pad4–1 sid2–1 (Additional file 6: Table S5 and Table S6). Therefore, IOS1-CERK1 may crosstalk with SARD1-CBP60g through PRRs and contribute to basal defense against tobacco powdery mildew (Fig. 6).

Finally, CRKs are one large subfamily of RLKs with 44 members in Arabidopsis and play important roles in response to abiotic and biotic stresses (Wrzaczek et al. 2010). We detected 16 CRKs that were expressed higher in Col-0 than in pad4–1 sid2–1 over the time course of infection, including CRK1, CRK4, CRK5, CRK6, CRK7, CRK8, CRK14, CRK15, CRK23, CRK24, CRK26, CRK36, CRK37, CRK38, CRK39 and CRK45 (Additional file 6: Table S5 and Additional file 7: Table S8). Previously, CRK4, CRK5, CRK6, CRK7, CRK8, CRK14, CRK15, CRK23, CRK24, CRK36, CRK37, CRK38 and CRK39 were found to be associated with cell death in bak1/serk4 mutant through regulation of ER quality control (ERQC) (de Oliveira et al. 2016). Overexpression of certain CRKs, such as CRK4, CRK5, CRK13, CRK19 and CRK20, can induce cell death in Arabidopsis transgenic plants (Chen et al. 2003; Chen et al. 2004; Acharya et al. 2007; de Oliveira et al. 2016). Moreover, overexpression of other CRKs, such as CRK4, CRK6 and CRK36 also leads to enhanced PTI responses (Yeh et al. 2015). Some CRKs are associated with SA-signaling pathway. SA can induce the expression of CRK4, CRK5, CRK6, CRK10, CRK11, CRK19, CRK20 and CRK45, while overexpression of CRK5, CRK13 and CRK20 results in higher amount of SA and up-regulation of PR1, PR5 and ICS1, marker genes of the SA-signaling pathway (Chen et al. 2003; Acharya et al. 2007; Ederli et al. 2011; Zhang et al. 2013). Therefore, CRKs may also play roles in basal defense and cooperate with SA-signaling pathway to activate post-penetration defense against tobacco powdery mildew, and some of them may contribute to cell death at late infection stage in Col-0 (Fig. 6).

Taken together, our data indicate that multiple PTI signaling pathways may constitute post-penetration defense against tobacco powdery mildew in Arabidopsis (Fig. 6).

Many SA-signaling-associated genes are up-regulated in Col-0 upon Gc SICAU1 infection

Previously, several labs have reported that SA-signaling pathway is required for host defense against powdery mildew (Nishimura et al. 2003; Xiao et al. 2005; Fabro et al. 2008). SA-signaling pathway also contributes to post-penetration defense against tobacco powdery mildew in Arabidopsis (Zhang et al. 2015). Here, we detected many SA-signaling-associated genes that were expressed higher in Col-0 than in pad4–1 sid2–1. First, several well-known SA-signaling pathway marker genes were expressed much higher in Col-0 than in pad4–1 sid2–1 upon infection of Gc SICAU1, including the Pathogenesis-Related gene 1 (PR1), PR2 (BGL2), PR5, DLO1, DOX1, NIMIN1, PHYTOALEXIN DEFICIENT 3 (PAD3), ACD6 and LURP1, and the SA synthetic pathway genes SID2 and EDS5 (Fig. 5a and Additional file 6: Table S6). Second, SA-signaling pathway associates with many WRKY transcription factors. We detected 7 WRKY transcription factors, i.e. WRKY38 (AT5G22570), WRKY51 (AT5G64810), WRKY53 (AT4G23810), WRKY54 (AT2G40750), WRKY58 (AT3G01080), WRKY62 (AT5G01900) and WRKY70 (AT3G56400), whose expression were higher in Col-0 than in pad4–1 sid2–1 over the time course (Additional file 6: Table S6). While WRKY70 acts as an activator of SA-signaling pathway but a repressor of JA-signaling pathway, WRKY38 and WRKY62 negatively regulate SA-signaling pathway, WRKY53 and WRKY54 are SA-responsive transcription factor acting as a positive and negative regulator of leaf senescence, respectively (Li et al. 2004; Li et al. 2006; Kim et al. 2008; Besseau et al. 2012). The expression of WRKY70 is positively regulated by the R2R3 Myb family TF MYB44 (AT5G67300) (Shim et al. 2013). Consistently, we found that MYB44 was down-regulated in both Col-0 and pad4–1 sid2–1 over the time course (Additional file 4: Table S3). Third, the Glutaredoxin-encoding gene GRX480 associates with TGA transcription factors and is involved in SA/JA cross-talk by suppressing JA-responsive PDF1.2 transcription (Ndamukong et al. 2007). GRXS13 is also required to maintain proper levels of superoxide radicals to protect cells from oxidative damage (Laporte et al. 2012). Here, we detected that the expression of GRX480 and GRX13 were higher in Col-0 than in pad4–1 sid2–1 over the time course (Additional file 6: Table S6). Fourth, Downy Mildew Resistant 6 (DMR6) and its closest homologs DMR6-Like Oxygenases (DLOs), DLO1 and DLO2, are responsive to SA. They negatively promote susceptibility to the downy mildew Hyaloperonospora arabidopsidis, the bacterium P. syringae and the oomycete Phytophthora capsici (Zeilmaker et al. 2015). Here, both DMR6 and DLO1 were up-regulated in Col-0, but down-regulated in pad4–1 sid2–1 upon infection of Gc SICAU1 (Additional file 6: Table S6 and Additional file 4: Table S3), implying that they play roles in response to tobacco powdery mildew.

Detoxification pathway may play critical roles for the sustained growth of tobacco powdery mildew in pad4–1 sid2–1

To scavenge the endogenous and/or exogenous toxic substances, most organisms exploit a three-phase detoxification system (Lallement et al. 2014). In phase one, enzymes such as cytochrome P450 monooxygenases (CYPs) catalyze reactions to introduce a functional moiety to hydrophobic substrates. In phase two, the newly modified substrates are conjugated with glutathione by specific classes of GSTs. Then the glutathionylated products are either sequestered in the vacuole or exported to the apoplast. Therefore, CYPs, hyroxylases and peroxidases act as the major players in phase 1, GSTs and glycosyl transferease act in phase 2, whereas, exocytosis and apoplastic signaling are required in phase 3. In our transcriptome data, we detected 3 CYP, 5 GST, 7 UGT genes, and 7 genes in apoplastic signaling that were expressed higher in pad4–1 sid2–1 than in Col-0 (Additional file 6: Table S7). While CYPs, GSTs, and UGTs may act in phase 1–2, apoplastic signaling genes may act in phase 3. Thus, our data imply that the sustained growth of Gc SICAU1 in pad4–1 sid2–1 may rely on the activation of the detoxification pathway (Fig. 6).

Secondary metabolites are involved in diverse stress-induced responses. These compounds rarely accumulate in their free form but are often conjugated to Glc through the action of glycosyltransferases that are encoded by UGT gene family. UGT members are involved in stress-induced responses (Langlois-Meurinne et al. 2005). Here, we detected seven UGT genes that were induced over the time course of infection with expression higher in pad4–1 sid2–1 than in Col-0 (Additional file 6: Table S7) and the other 27 UGTs that were differentially expressed in Col-0 and pad4–1 sid2–1 upon infection of Gc SICAU1, most of them expressed higher in pad4–1 sid2–1 than in Col-0 (Additional file 8: Table S9). Among them, UGT78D1 (AT1G30530) encodes a UDP-rhamnose:flavonol-3-Orhamnosyltransferase involved in flavonol glycoside biosynthesis (Jones et al. 2003). Interestingly, most of the above mentioned UGTs belong to Arabidopsis Glycolsytransferase Family 1 that were reported to be significantly up-regulated by treatment of 2,4,6-trinitrotoluene (TNT) (Gandia-Herrero et al. 2008). Therefore, UGT-mediated detoxification pathway genes may play roles for the sustained growth of Gc SICAU in the pad4–1 sid2–1 mutant (Fig. 6).

CLE (CLAVATA3/Endosperm surrounding region-related) family genes encode secreted peptides acting in regulation of plant growth, development and response to environmental stimuli (Wang et al. 2015). Previously, CLE1 was found to act in signaling regulating root development in nitrogen-dependent manner (Araya et al. 2014). In our data, we identified that CLE1 (AT1G73165) was highly and particularly induced in pad4–1 sid2–1 from 6 to 12 dpi (Fig. 5f and Additional file 6: Table S7). Thus, CLE1 may play a role in apoplast signaling to regulate the sustained growth of Gc SICAU in pad4–1 sid2–1 mutant and this could be an interesting focus in future research.

It seems that the genes involved in the detoxification pathway were also activated in Col-0 at earlier time points (1–3 dpi) but dropped back to background level presumably due to the activation of defense-related signaling pathways. Particularly, SA signaling pathway might suppress the activation of detoxification. On the contrary, genes involved in detoxification may be continually increased in pad4–1 sid2–1 due to the lack of SA signaling (Fig. 6).

Conclusions

Powdery mildew fungi are complex pathogens that infect many agriculturally important plants such as wheat, barley and rye. In a previous study, we identified the tobacco PM strain Gc SICAU1 that can grow well in pad4–1 sid2–1 but not in Col-0 (Zhang et al. 2015). Here, through transcriptome analysis on a time course of infection, we found that 3 dpi was important for activation of post-penetration defense in Col-0 and 6 dpi was crucial for activation of the detoxification pathway for the sustained growth in pad4–1 sid2–1. We also identified a number of genes that may contribute to post-penetration defense and the sustained growth. From these data and combining with the function of related genes reported in the literatures, we proposed a working hypothesis for the post-penetration defense in Col-0 and the sustained growth in pad4–1 sid2–1 (Fig. 6). First, upon infection of the tobacco powdery mildew in Col-0, a series of defense pathways may be activated. Quite a few receptors, such as LYK5, WAK1, RLP13 and RLK5, that recognize chitin, OGs, NLPs and IDA, respectively, were induced over the time course of infection. Subsequently, genes involved in MAPK cascade were also up-regulated, which in turn presumably led to activating basal defense. Meanwhile, a batch of CRKs was up-regulated, which probably led to activation of SA-signaling pathway. Eventually, these defense responses were culminated into cell death, which arrested the growth of the tobacco powdery mildew. On the contrary, because basal defense and SA-signaling were defective in pad4–1 sid2–1, a variety of genes involved in the detoxification pathway were remarkably up-regulated particularly at and after 6 dpi of Gc SICAU1. The detoxification pathway may be antagonistic to SA-signaling pathway and presumably recruited members of CYP, GST and UGT gene families that are involved in toxin catabolic process, proline/amino acid transport, flavonoid biosynthetic process and apoplastic signaling. In turn, the detoxification pathway led to the sustained growth of Gc SICAU1 in pad4–1 sid2–1. In theory, genes involved in the sustained growth of Gc SICAU1 could be the candidates for engineering disease resistance by gene-knock-out approaches and thus are significant in future investigations.

Methods

Plant materials and powdery mildew inoculation

Arabidopsis thaliana wild type (WT) accession Col-0 and the double mutant pad4–1 sid2–1 were grown under 8 h/16 h light/dark regime at 22 °C. Tobacco powdery mildew isolate Gc SICAU1 was maintained on tobacco leaves at 23 °C (16 h light, 8 h dark) in a growth room as previously reported (Zhang et al. 2015). Rosette leaves from 6-week-old plants were inoculated with powdery mildew strain Gc SICAU1. Inoculum was prepared by growing of Gc SICAU1 on tobacco leaves for 10 days, and then spores were collected to inoculate on WT and pad4–1 sid2–1. Samples were collected at 0, 1, 3, 6, 8, 10 and 12 dpi for transcriptome analysis.

RNA-seq analysis

RNA was extracted from two biologically duplicated leaf samples using TRIzol® Reagent (Life Technologies, USA) following the manufacturer’s protocol. RNA quantity and quality were determined using a 2100 BioAnalyzer (Agilent Technologies Canada Inc., Mississauga, ON, Canada). RNA libraries were constructed with 2 μg total RNA and subjected to high-throughput sequencing by a HiSeq-4000 sequencer (Illumina). Approximately 6 Gb of reads for each sample were obtained. After filtering the raw reads by removing adapter and reads with 5% “N”, clean data were mapped to the TAIR10 reference genome by TopHat2 with default parameters (Kim et al. 2013).

The gene quantity was calculated by HTseq (Anders et al. 2015). The FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) method (Lee et al. 2011) was used to calculate the normalized expression data of each gene (Mortazavi et al. 2008). DEGs were selected by DEseq2 (Love et al. 2014), with the criteria of absolute log₂ fold change ≥2 and false discovery rate (FDR) ≤ 0.01 between any two samples. GO (Ashburner et al. 2000) database was used for gene ontology analysis. GO enrichment analysis was performed using GOseq (Young et al. 2010).

Clustering DEGs

Hierarchical clustering analysis of DEGs was run with hclust in R version 3.2.0. The FPKM of genes were log-converted and centered, and the complete linkage method was used. The median of log₂(FPKM+ 1) for each transcript was subtracted from the log₂(FPKM+ 1) value of each sample. Then the distance matrix was run with hclust in R, and a dendrogram was plotted that displays a hierarchical relationship among the vehicles. The hclust function in R uses the complete linkage method for hierarchical clustering by default. At every stage of the clustering process, the two nearest clusters were merged into a new cluster. The process was repeated until the whole data set was agglomerated into one single cluster. The final clusters were generated by cutting trees at 20% of the height.

Quantitative RT-PCR analysis

To validate the transcriptome data, reverse-transcription quantitative PCR (RT-qPCR) was performed for 19 genes representing different expression patterns. cDNA was synthesized from 0.5 μg total RNA using ReverTraAce qPCR RT Master Mix with gDNA Remover (TOYOBO, Japan). A two-step real-time PCR reaction was performed using an QuantiFast™ SYBR® QPCR kit (QIAGEN, Germany) with 100 ng template DNA and 10 nM each primer in a final volume of 10 μL according to the following protocol: polymerase activation at 95 °C for 5 min followed by denaturation at 95 °C for 10 s and concurrent annealing and extension at 65 °C for 30 s. The Arabidopsis actin (ACT2) gene was used as the internal control. Primer sequences were listed in Additional file 10: Table S10.