Arabidopsis Defense against Botrytis cinerea: Chronology and Regulation Deciphered by High-Resolution Temporal Transcriptomic Analysis

High-Resolution Time Series Expression Profiling Identifies 9838 Differentially Expressed Arabidopsis Genes Following B. cinerea Infection

Full genome expression profiles were obtained from Arabidopsis leaves following infection with B. cinerea, a fungal necrotroph. Leaf 7 was detached from 192 4-week-old Arabidopsis plants and either inoculated with a suspension of B. cinerea spores or mock inoculated. Over 48 HAI, expanding lesions developed on the pathogen-inoculated leaves (Figure 1). Every 2 HAI (up to 48 HAI) four whole leaves were harvested from each treatment. Expression analysis was performed using CATMA (a complete Arabidopsis transcriptome microarray) arrays (Sclep et al., 2007), cDNA from single leaves, and a statistically designed loop design of hybridizations (see Supplemental Figure 1 online), leading to a high-resolution, highly replicated time series of expression profiles (24 time points separated by 2 h; four biological and an average of three technical replicates for each time point in each condition).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Time Series of B. cinerea Infection on Arabidopsis Leaf 7.

(A) Ten-microliter droplets of a suspension of B. cinerea spores (1 × 10⁵ spores mL⁻¹) were placed on detached leaf 7 from 4-week-old Arabidopsis plants. Images show the same leaf every 2 h after inoculation until 48 h.

(B) A mock-inoculated leaf at 2 HAI (left) and 48 HAI (right).

Bars = 10 mm.

Expression values for each CATMA probe at each time point in each biological replicate were extracted using a mixed-model analysis in a locally adapted version of the R package MAANOVA (for microarray analysis of variance; Wu et al., 2003). This time series data set is longitudinal in that the data reflect the defense process over time, but also cross-sectional in that each sample was one leaf from a different plant (i.e., destructive sampling) so there is no particular connection between the individual biological replicates. Due to this hybrid nature of the data, we investigated three statistical tests for their ability to determine genes differentially expressed between mock-inoculated and B. cinerea–infected samples over time. A detailed description of this process is given in Methods. In brief, we combined the outputs of a standard F test within MAANOVA with that of a Gaussian process two-sample test (GP2S) (Stegle et al., 2010), based on the low false positive rates of these methods (see Supplemental Figure 2 online). We combined the top 10,600 gene probes ranked by GP2S with 236 additional gene probes identified by the F test. Probes that did not map to genes in the TAIR9 annotation and duplicate probes (two or more probes mapping to the same gene) were removed. Thus, the time series expression profiling identified 9838 Arabidopsis genes as differentially expressed between B. cinerea–infected and mock-inoculated leaves over time.

The expression profiles for each individual probe on the CATMA array can be viewed using a Web tool (under the “data” section at http://go.warwick.ac.uk/presta). This plots the expression profiles at all 24 time points for both the infected and mock-inoculated leaves. Variation in expression is shown as a bar representing one standard error. The full data set is available in the Gene Expression Omnibus (GEO) under accession number GSE29642.

As an initial validation of the data set, the profiles of genes previously identified as differentially expressed during B. cinerea infection (Mengiste et al., 2003; AbuQamar et al., 2006; Pré et al., 2008; Dhawan et al., 2009; Wang et al., 2009; Chen et al., 2010; Luo et al., 2010), and in several cases known to influence the progression of disease, were examined (see Supplemental Table 1 online). In the majority of cases (28/41), the genes were identified as differentially expressed in our time series and expression profiles matched that in the literature. In a further four cases, the genes were ranked below the GP2S cutoff, but manual inspection showed that they were differentially expressed and again the profiles matched those in the literature (see Supplemental Table 1 and Supplemental Figure 3A online). Three genes reported to be upregulated in the literature did not show differential expression in our time series (see Supplemental Figure 3B online). Intriguingly, six genes showed differential expression in the opposite direction in our study compared with the literature (see Supplemental Figure 3C online). These included two genes, ANAC002/ATAF1 and HUB1, whose expression level influences defense against B. cinerea (Dhawan et al., 2009; Wang et al., 2009). Using RT-PCR, we tested expression of ATAF1 and LOX2 in RNA from two of the four biological replicate samples at eight time points (see Supplemental Figure 3D online). The profiles matched those from the whole time series, indicating that the difference in expression was not due to the probes on the CATMA arrays; it seems that expression of even key genes varies depending on the environmental conditions, infection strategy, or isolate of B. cinerea being used.

Obviously, changes in transcription are not the only regulatory mechanism employed by plants to regulate their immune response. Our subsequent analysis of these differentially expressed genes (DEGs) and interpretation of such analyses is based solely on transcriptional events, although other regulatory mechanisms are possible.

The Time Series Expression Profiling Spans Multiple Stages of Infection

B. cinerea infection was initiated by pipetting droplets of spore suspension onto the top surface of detached leaves. The first visual symptoms of infection at 20 HAI are a darkening of the leaf surface under the inoculum droplets (Figure 1) and correspond to primary lesion formation following penetration of the host. Expansion of the lesion beyond the inoculum droplets is evident at 36 HAI and continues throughout the 48-h sampling period. We determined the expression of the B. cinerea β-tubulin gene relative to a nonchanging Arabidopsis gene (PUX1, At3g27310) as a measure of fungal growth (Figure 2). An initial rapid increase in tubulin expression/fungal biomass can be attributed to germination of conidiophores and hyphal growth in the inoculum media. A lag phase in growth is apparent between 20 and 28 HAI, during which time initial lesion formation occurs. Trypan blue staining of infected leaf tissue in the middle of the lag phase showed fungal hyphae as well as a claw arrangement of much thicker tubular structures (Figure 2B). These claw-like structures have been associated with penetration of the host and appear to develop from hyphae rather than undifferentiated germ tubes (Kunz et al., 2006). The lesion expansion stage appears to begin by 32 HAI with fungal biomass once again increasing and lesion expansion visible on the leaves from 36 HAI.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Growth of B. cinerea in Arabidopsis Leaves.

(A) Growth measured by expression of the B. cinerea tubulin gene. Expression levels were determined using real-time PCR and are shown as the log₂ ratio of expression of B. cinerea Tubulin relative to Arabidopsis PUX1 (At3g27310). Error bars indicate se of three biological replicates. A one-way analysis of variance was performed to determine which pairs of adjacent time points differed significantly from each other. Significantly different groups (P ≤ 0.05) are labeled a to e.

(B) Trypan blue staining of Arabidopsis leaves infected with B. cinerea 24 HAI. Red arrows: I, filamentous hyphae; II, large claw-like structures. Bar = 25 μm

Clustering of DEG Expression Profiles Reveals Coexpression of Functionally Related Genes

To look at the overall patterns in gene expression during the infection process, the 9838 DEGs were clustered using the SplineCluster algorithm (Heard et al., 2005) on the basis of their expression in infected leaves. Some of these DEGs show diurnal variation in expression in the mock-inoculated leaves (see the section “B. cinerea infection dampens clock gene oscillations” below); however, changes in response to infection override diurnal patterns, hence clustering on the basis of expression profile during infection is valid. Using a prior precision value of 0.001, 44 clusters were obtained that are shown in Figure 3 (two of which are singleton clusters). From a heat map of these 44 clusters, it is clear that a major shift in gene expression (up and down) of infected leaves occurs by ∼26 HAI (see Supplemental Figure 4 online). This major transcriptome change occurs before visible lesion formation during the lag phase of B. cinerea growth (Figure 2). Despite this major change, clusters of genes whose expression changes earlier or later than this point in infection and clusters showing transient changes in expression are also visible. The list of genes in each cluster is given in Supplemental Data Set 1 online.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Infected Leaf Expression Profiles of the Gene Members of Each of the 44 Clusters Generated by SplineCluster.

The individual gene profiles are shown as gray lines with the mean profile in dark blue. The dashed blue lines indicate the mean ± 1 sd of the cluster. The y axes indicate log₂ expression normalized on a per gene basis. hpi, hours post inoculation.

The 44 clusters represent groups of genes that are coexpressed over the time course of infection. We analyzed these groups for overrepresentation of Gene Ontology (GO) terms (Ashburner et al., 2000) using BiNGO (Maere et al., 2005) to ask whether genes in the same cluster are involved in the same biological process, suggesting coordinated regulation of the process. Many different terms were overrepresented in these clusters (see Supplemental Data Set 2 online), suggesting coregulation of genes and highlighting the breadth of the response to this pathogen.

Chronology of the Defense Response

We specifically wanted to identify biological processes that were taking place during the early stages of B. cinerea infection, since these are more likely to influence the outcome of the plant–pathogen interaction and the chronology of the defense response. Two methods were used to investigate the timing of gene expression changes. First, a Gaussian process regression analysis was used to identify the time point at which there is a change in the rate of each gene’s expression. A gradient significantly greater or less than zero indicates expression of the gene is increasing or decreasing, respectively; a gradient of zero indicates a steady level of expression. This analysis was described in detail by Breeze et al. (2011). From the gradient information for each gene, the first time at which at least half of the genes in a cluster have a significantly increasing or decreasing gradient was calculated. This gave a single time point for each cluster indicating the time at which expression of the genes in that cluster began to change following infection (see Supplemental Data Set 2 online) and enables us to chronologically order the biological processes identified through GO analysis of the SplineClusters (Figure 4).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Selected GO Terms Overrepresented in Clusters of Genes Differentially Expressed after B. cinerea Infection of Arabidopsis Leaves.

GO terms are aligned with the time of gradient change and/or time of first differential expression of the cluster (in italics), with red boxes containing GO terms from upregulated genes and blue boxes containing GO terms from downregulated genes.

SplineCluster groups genes on the basis of similarity of their expression profiles over time. We were also interested in whether genes could be grouped in a meaningful way using the time at which a gene is first differentially expressed after infection. A time-local version of the GP2S test (Stegle et al., 2010) was used to determine the time at which each of the 9838 DEGs was first differentially expressed in the B. cinerea–infected leaves compared with mock inoculated leaves. Seventy-four genes were not identified as differentially expressed using this model, but a time of first differential expression (TOFDE) was determined for the remaining 9764 genes. In the time-local GP2S, the 48-h time series was split into 100 increments; hence, TOFDE was calculated to the nearest half hour. TOFDE was used to group the DEGs in bins of 30 min, 1 h, or 2 h. These represent groups of genes that respond to B. cinerea infection of Arabidopsis leaves at the same time and were analyzed for overrepresented GO terms again using BiNGO (Maere et al., 2005). GO terms overrepresented in specific time bins are listed in Supplemental Data Set 3 online with selected terms again highlighted in Figure 4. TOFDE does not separate genes differentially expressed around the middle of the time series as well as the gradient tool analysis; however, it did highlight some additional processes occurring during infection.

Signaling

Key events in many plant responses are the synthesis and/or response to phytohormones. The involvement of ET, auxin, abscisic acid (ABA), and JA was highlighted by overrepresented GO terms. Arabidopsis defense against necrotrophic pathogens, including B. cinerea, is known to involve or be affected by ET, auxin, ABA, and JA (Thomma et al., 1998, 1999; Audenaert et al., 2002; Pandey et al., 2005), but our analysis enables the order of synthesis and/or action of these hormones to be elucidated. At 14 HAI, two genes encoding 1-aminocyclopropane-1-carboxylate synthases (ACS2 and ACS6), enzymes catalyzing the first and rate-limiting step in ET biosynthesis, are upregulated. ACS2 and ACS6 proteins are known to be responsible for the majority of B. cinerea–induced ET production and are phosphorylated and stabilized by MPK3/6 (Han et al., 2010); however, genes encoding ACS enzymes are also known to be transcriptionally activated (Tsuchisaka and Theologis, 2004). Analysis of an acs2 acs6 double mutant suggests that another ACS protein is also involved in ET production in response to B. cinerea infection (Han et al., 2010). However, although all nine ACS genes were on the CATMA arrays, only ACS2 and ACS6 are differentially regulated in our analysis. Synthesis of ET should lead to downstream events; hence, the overrepresentation of the GO term “response to ET” 2 h later and “ET-mediated signaling” highlighted in the time bins analysis. Genes responsible for these terms include EBF2 with a known role in ET signaling (Saito et al., 2004) and two TFs from the ET response factor family, At ERF15 and ORA59.

“Response to JA stimulus” is a GO term overrepresented in the same early cluster of genes responding to ET (16 HAI). Genes corresponding to this term include ERF1 and MYB108, for which overexpressor and knockout lines, respectively, show altered B. cinerea susceptibility (Berrocal-Lobo et al., 2002; Mengiste et al., 2003), four genes with a known role in defense (PROPEP1, ERF4, PEN1, and MYB51) (Huffaker et al., 2006; Kwon et al., 2008; Pré et al., 2008; Clay et al., 2009), and MYB13, whose expression is induced by B. cinerea in a JA/ET-dependent manner (AbuQamar et al., 2006). The remaining gene, VTC5, has no known defense role, but its coexpression with these other defense regulators makes this a viable hypothesis. Interestingly, ERF4 mediates antagonism between the ET and ABA pathways with overexpression of ERF4 leading to decreased sensitivity to ABA (Yang et al., 2005).

“Auxin biosynthesis” is overrepresented in genes upregulated 22 HAI. This group of genes includes anthranilate synthase (ASA1), a rate-limiting step in the synthesis of the auxin precursor Trp, STY1, a transcriptional activator of auxin biosynthesis (Eklund et al., 2010), and two paralogous genes (RGLG1 and 2) thought to be responsible for directional flow of auxin (Yin et al., 2007). At least in roots, ASA1 and ASB1 are required for ET-mediated increases in auxin (Stepanova et al., 2005). Hence, the earlier synthesis of ET we observe suggests that a similar mechanism is operating during response to B. cinerea infection; ET activates auxin biosynthesis via ASA1.

ABA-associated GO terms suggest a strong repression of ABA signaling during infection by B. cinerea (Figure 4). ABA catabolism is overrepresented in the group of genes first differentially expressed 20 HAI due to the upregulation of CYP707A3 and UGT71B6. CYP707A3 catalyzes both 8′- and 9′-hydroxylation of ABA, with 8′-hydroxylation being the major pathway of ABA catabolism (Saito et al., 2004; Okamoto et al., 2011), while UGT71B6 glycosylates ABA to supposedly inactive conjugates (Priest et al., 2005). Two hours later (22 HAI) upregulation of negative regulators of ABA signaling begins. ABI1 and ABI2, two protein phosphatases involved in the core ABA pathway, are induced along with two repressors of ABA responses: ABR1, thought to be a transcriptional repressor (Pandey et al., 2005), and TMAC2, a nuclear-localized protein (Huang and Wu, 2007). AZF2 also has a TOFDE of 22 HAI and is another repressor of ABA signaling (Drechsel et al., 2010). Three of the genes encoding the soluble PYL/PYR/RCAR ABA receptors (PYL8, PYL9/RCAR1, and PYR1) are grouped in cluster 12 and begin to change in expression 18 HAI. All three genes are downregulated, lending weight to the hypothesis that ABA signaling is repressed during B. cinerea infection.

Plant hormones play a major role in defense, and this analysis has provided a timeline of the synthesis and/or action of these during infection. However, other signaling mechanisms are also highlighted by the GO term analysis of clusters. The term “lipid transport” is overrepresented in a cluster of genes that are upregulated very early after infection (10 HAI, cluster 6). Genes annotated with this term include two nonspecific lipid transfer proteins (nsLTPs), LtpV.2 and LtpV.3 (Boutrot et al., 2008), and the two xylogen proteins in Arabidopsis, XYP1 and XYP2, which also contain a nsLTP domain (Motose et al., 2004). XYP1 in particular responds early to B. cinerea infection, and all four genes are downregulated around 25 HAI. The physiological function of nsLTPs is not well understood, and only a few have been demonstrated to bind lipids. They are thought to have a defense function and as such have been characterized as pathogenesis-related protein family 14 (van Loon and van Strien, 1999). However, XYP1 and XYP2 have a role in vascular differentiation (Motose et al., 2004) and have not previously been implicated in plant defense. The coordinated expression of several nsLTP-containing proteins suggests they have a specific role in defense.

Crosstalk between signaling pathways mediating plant responses to biotic and abiotic stress is well known with a growing number of genes involved in these interactions being identified (Fujita et al., 2006). This crosstalk is also evident from our B. cinerea infection time series expression data. “Response to abiotic stress” and more specific child terms are overrepresented in several clusters of DEGs (Figure 5). These clusters include both up- and downregulated genes and span both the primary lesion formation and lag phases of infection. Gradient analysis gave clusters 5, 6, and 34 change times of 8, 10, and 6 HAI, respectively. However, the TOFDE for each of the genes annotated with an abiotic stress GO term was much later, and manual inspection of the gene profiles indicated that the TOFDE values were correct; hence, this is the timing used in Figure 5. This discrepancy occurs because many of these genes show diurnal patterns of expression hence the mock expression profiles are also changing over time.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

GO Terms Relating to Abiotic Stress Responses Overrepresented in Clusters of Genes Differentially Expressed after B. cinerea Infection of Arabidopsis Leaves.

GO terms are aligned with the time of gradient change and/or time of first differential expression of the cluster, with the cluster expression profile shown. All “response to abiotic stimulus” and nonredundant individual abiotic stress GO terms are shown.

Comparison of Gene Expression Patterns during B. cinerea Infection and Developmental Leaf Senescence Shows Considerable Overlap but Reveals Specific Features

Previous work has analyzed a similar time series of gene expression changes during leaf senescence in Arabidopsis and identified over 6000 genes showing differential expression over the 22 d from leaf expansion to senescence (Breeze et al., 2011). These data provide an exciting opportunity to compare and contrast gene expression changes between senescence and defense against B. cinerea infection. The two plant responses are very different in timing but may involve similar signaling pathways. Lists of gene differentially expressed in each treatment were divided into up- or downregulated genes. For the senescence data from Breeze et al. (2011), clusters 1 to 24 were classed as downregulated and clusters 27 to 48 as upregulated, with genes in clusters 25 and 26 being omitted from the analysis due to them showing both up- and downregulation. For genes differentially expressed in response to B. cinerea infection, clusters 23 to 44 were classed as upregulated and clusters 1 to 22 as downregulated.

Overlapping genes in the four lists (B. cinerea up/down and senescence up/down) were identified. A total of 3759 genes showed differential expression in both data sets; however, the vast majority of genes (8126) were specifically differentially expressed under one condition only (Figure 6). In most cases, genes differentially expressed during both senescence and B. cinerea infection were similarly regulated in each response: 1405 genes upregulated and 1767 downregulated. However, 502 genes were upregulated during senescence but downregulated following B. cinerea infection, while only a small group showed the opposite profile (85 genes). Genes in all eight groups are listed in Supplemental Data Set 4 online.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Number and Function of Genes Differentially Expressed during Both Natural Senescence and B. cinerea Infection.

The number of genes up- and downregulated during senescence and B. cinerea infection and overlaps between these is shown in the Venn diagram. Selected overrepresented GO terms are shown for each subset of genes.

[See online article for color version of this figure.]

GO term analysis was applied to each set of genes to identify pathways and functions common or specific to the two responses. Selected overrepresented terms for each group of genes are shown in Figure 6 with the full set listed in Supplemental Data Set 5 online. As expected, very high over representation of genes involved in photosynthesis, chlorophyll biosynthesis, and starch metabolism was observed in the downregulated genes common to both processes. Genes responsive to ABA, ET, JA, and SA are all overrepresented in genes upregulated in both senescence and B. cinerea infection, highlighting the important role plant hormones play in both of these stress responses. However, there are still condition-specific aspects of hormone involvement indicating complex differential activation in the two processes. Different groups of ABA-responsive genes are upregulated during senescence and B. cinerea infection, upregulated in senescence only, and upregulated in senescence but downregulated during B. cinerea infection. Similarly, ET-responsive genes are overrepresented in both the senescence and B. cinerea upregulated group and B. cinerea infection only, while different SA-responsive genes are upregulated in both conditions as well as downregulated during B. cinerea infection only. Markedly, the involvement of cytokinin and brassinosteroid hormones appears specific to senescence. The specific downregulation of several Arabidopsis response regulator genes and cytokinin response factors during senescence, and not during B. cinerea infection, suggests that the repression of the photosynthetic machinery that occurs during both senescence and infection is not dependent on changes in cytokinin levels.

The aromatic amino acid biosynthesis pathway appears to be important only during B. cinerea infection, and genes encoding enzymes involved in phenylpropanoid, chorismate, Trp, Tyr, and Phe biosynthesis are all upregulated. One of the roles of the Trp pathway during pathogen infection is to provide substrate for camalexin synthesis. The antimicrobial activity of camalexin is clearly not relevant to the senescence process. Upregulation of genes involved in glutathione metabolism also appears specific to B. cinerea infection. Genes encoding glutathione biosynthetic enzymes as well as GSTs and enzymes involved in the reduction of oxidized glutathione and degradation of glutathione conjugates are only upregulated during pathogen infection. This may well reflect the increased burden of oxidative stress that the plant has to cope with during B. cinerea infection, as well as the presence of many toxins. It may also reflect the demand for camalexin, which uses glutathione in its production.

Despite the common involvement of several plant hormones, many genes involved in the regulation of transcription are induced during senescence only or differentially expressed during senescence (up) and B. cinerea infection (down). The NF-Y TF family (corresponding to the CCAAT binding factor complex) is clearly regulated very differently in the two responses. In mammals, a heterotrimeric complex of NF-YA, NF-YB, and NF-YC subunits is required for DNA binding of these TFs to the CCAAT motif. Of the 36 NF-Y genes in Arabidopsis, 26 are differentially expressed during senescence and/or B. cinerea infection. Only one gene, NF-YB4, is similarly expressed (downregulated) during both conditions. The remaining 25 genes are differentially regulated during senescence and B. cinerea infection; six are upregulated and one downregulated in only senescence, nine upregulated and two downregulated in just B. cinerea infection, and seven upregulated in senescence and downregulated in pathogen infection. This family of TFs therefore appears to be key determinants of regulatory specificity in these stress responses. Another level of specificity is also apparent; many NF-YA subunits show increased expression during senescence, while NF-YB and NF-YC subunits are upregulated following B. cinerea infection.

B. cinerea Infection Dampens Clock Gene Oscillations

There has been much circumstantial evidence about the influence of the circadian clock on pathogen defense (Roden and Ingle, 2009). However, a small polypeptide PCC1, whose overexpression leads to resistance against Hyaloperonospora arabidopsidis, is under both defense and circadian regulation (Sauerbrunn and Schlaich, 2004), and, recently, a group of genes involved in basal and R-mediated defense was shown to be under the control of CCA1, a core clock component (Wang et al., 2011). In the latter study, the central clock oscillator appeared unaffected by H. arabidopsidis infection. By contrast, B. cinerea infection of Arabidopsis leaves appears to influence expression of core clock genes, suggesting a stress input into the central oscillator. The timing of core clock gene expression is not perturbed, but the amplitude of cyclical expression is reduced (Figure 7). The reduction in amplitude affects genes expressed at different phases of the clock and appears to start around 24 HAI.

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

B. cinerea Infection Dampens Oscillations of Clock Gene Expression.

The green line indicates the mean expression profile of four biological replicates of mock-inoculated leaves, while the red line indicates the mean expression profile of four biological replicates of infected leaves. Error bars indicate se (n = 4).

In the mock-inoculated gene expression data, we identified 2404 genes that were expressed in a rhythmic fashion with an ∼24-h period (see Supplemental Data Set 6 online). GO analysis of these genes indicated that, as expected, the annotation “circadian rhythm” was significantly overrepresented. “Response to stress,” “response to abiotic stress,” and several terms associated with plant stress responses (for example, “response to cold,” “response to ABA stimulus,” and “response to reactive oxygen species”) were also overrepresented in these rhythmic genes (see Supplemental Data Set 7 online). This is consistent with previous reports (Covington et al., 2008). When these 2404 genes were grouped according to their phase of rhythmic gene expression (using 2-h intervals), the GO term “response to abiotic stimulus” was significantly overrepresented in genes peaking in expression 8 to 12 h after dawn. Genes peaking in expression 10 to 12 h after dawn were overrepresented for several defense-related terms, including “JA biosynthesis,” “response to fungus,” and “response to biotic stress” (see Supplemental Data Set 8 online). Hence, in contrast with the CCA1-controlled defense genes identified by Wang et al. (2011), defense-related genes appear to peak in the early afternoon in our experiment. Over 60% of the 2404 rhythmic genes from the mock-inoculated leaves were differentially expressed in response to B. cinerea infection (1521 genes). In the vast majority of these cases (1407 genes), expression in response to infection overrode rhythmic expression, and these genes were not classified as rhythmic (with a period between 20 and 28 h) in the expression data from infected leaves.

Specific TF Binding Motifs Are Enriched in Groups of Coexpressed Genes

Clustering of DEGs based on their expression profile during B. cinerea infection groups coexpressed genes together and as seen above clearly enables specific processes in defense to be delineated. To ask whether this clustering can also identify coregulated genes, we analyzed the frequency of known TF binding motifs in the promoters of genes in each cluster. Despite the large number of DEGs, many clusters were enriched for specific TF binding motifs (Figure 8; see Supplemental Data Set 9 online). Furthermore, clusters with similar expression profiles are specifically enriched for similar motifs, with a clear difference between the motifs enriched in downregulated (1 to 22) and upregulated (23 to 44) clusters.

• Download figure
• Open in new tab
• Download powerpoint

Figure 8.

Known cis-Regulatory Sequences Associate with Groups of Coexpressed Genes.

Regulatory motifs (represented by sequence logos where character size indicates nucleotide frequency) are differentially enriched in the promoters of genes clustered on the basis of their expression during B. cinerea infection. The blue shaded boxes correspond to raw P value. Expression profiles from selected gene clusters that are enriched for TF binding motifs are shown on the right. Full results used to derive this figure are shown in Supplemental Table 9 online.

Clusters of genes that are suddenly downregulated ∼20 HAI and associated with GO terms related to photosynthesis are enriched for variants of the G-box and I-box motifs, which have been shown to function in light regulated gene expression (Donald and Cashmore, 1990; Menkens et al., 1995). The apparently concerted downregulation of such a large number of genes suggests repression of photosynthetic genes is highly coordinated. The G-box motif has been shown to interact with TFs that act as repressors, suggesting that downregulation of photosynthetic genes could be mediated through these motifs. For example, phytochrome interacting factors bind to the G-box motif present in the promoters of some photosynthetic genes (Huq and Quail, 2002; Huq et al., 2004), and some phytochrome interacting factors, such as PIF3 (Shin et al., 2009) and PIF7 (Kidokoro et al., 2009), have been shown to negatively regulate expression of photosynthetic genes. However, both PIF3 and PIF7 are downregulated during B. cinerea infection, suggesting other TFs are playing this role.

Many of the downregulated clusters are enriched for the site-II motif [TGGGC(C/T)], which serves as a binding site for TCP TFs (Martín-Trillo and Cubas, 2010). In addition to this site-II motif, consensus sequence motifs for both Class I and Class II TCP TFs were found enriched in these specific clusters. TCP proteins are commonly known as regulators of cell proliferation, growth, and development (Li et al., 2005; Hervé et al., 2009; Aggarwal et al., 2011; Kieffer et al., 2011). It is possible that developmental processes may be repressed as part of the defense response to direct resources toward fighting infection. However, three upregulated clusters are also enriched for TCP binding motifs, and TCPs also play a role in the circadian clock (Pruneda-Paz et al., 2009; Giraud et al., 2010) and control of JA biosynthesis (Schommer et al., 2008). Recently, TCP binding sites were shown to be enriched in calcium-responsive gene promoters (Whalley et al., 2011), which may explain the large number of defense-related gene promoters containing TCP binding motifs.

Dof family TFs interact with motifs that are enriched within both down- and upregulated clusters. Dof TFs, of which several are differentially expressed following infection, act as transcriptional activators or repressors in a wide range of biological processes (Yanagisawa, 2004), although as for the TCPs, these mostly include growth and developmental processes. However, expression of a group of three Dof TFs (OBP1-3) is induced by auxin and SA, and the majority of genes differentially expressed in a transgenic line overexpressing OBP2 (Dof1.1) are involved in response to biotic stress (Skirycz et al., 2006). In particular, OBP2 regulates expression of several genes involved in indolic glucosinolate biosynthesis. Although these genes are also involved in the biosynthesis of camalexin and are induced in response to B. cinerea infection, OBP1-3 are not differentially expressed in response to infection. This suggests additional members of the Dof family also play a role in biotic stress.

Members of the WRKY TF family are known to be involved in regulation of plant defense responses as both positive and negative regulators (Eulgem and Somssich, 2007). The W-box motif, which is bound selectively by members of the WRKY TF family, is overrepresented in clusters of genes rapidly induced 18 to 24 HAI (Figure 8). W-boxes are known to be enriched within the promoters of genes induced by biotic stress, and many WRKY TFs have a demonstrated function in defense against B. cinerea and other pathogens (for example, WRKY3, 4, 46, 53, and 70; Lai et al., 2008; Hu et al., 2012). WRKY TFs are overrepresented in TFs upregulated at many time points following B. cinerea infection (Figure 9). Despite this, W-boxes are enriched in specific clusters with a similar expression pattern.

• Download figure
• Open in new tab
• Download powerpoint

Figure 9.

TF Families Significantly Overrepresented for DEGs, Indicating Distinct Periods of Regulation.

The plots show a number of TF families significantly upregulated ([A], red) and downregulated ([B], green) following B. cinerea infection. Color bars indicate P values (after FDR correction; Benjamini and Hochberg, 1995) with a range of significance thresholds (0.01, 0.05, 0.1, 0.25, and 0.5).

Fifty-three AP2/EREBP TFs are differentially expressed following B. cinerea infection; however, the GCC-box, which is bound by AP2/EREBP TFs, is only overrepresented in cluster 22. This could indicate that multiple AP2/EREBP TFs operate at different times during infection, mediating different target gene expression patterns, or that alternative motifs can be bound by these TFs. A MYB TF binding motif is enriched within the promoters of genes in cluster 26 that are upregulated in response to infection. The motif was identified as the binding site for MYB80; however, neither MYB80 nor the six other MYBs most closely related to MYB80 are differentially expressed after B. cinerea infection, suggesting other MYB TFs may be involved. Several members of the MYB family have regulatory roles in response to biotic stress (Mengiste et al., 2003; Clay et al., 2009; Ramírez et al., 2011a), and 36 of the 132 MYB TFs in Arabidopsis (Stracke et al., 2001) are differentially expressed in response to B. cinerea infection.

NAC TFs are heavily linked with the regulation of abiotic stress responses (Nakashima et al., 2012). However, recent studies have indicated that this family of TFs also plays crucial roles in the response to biotic stress. Transgenic lines with reduced or increased expression of ANAC002, ANAC019, ANAC055, ANAC081, and ANAC092 have altered susceptibility to pathogen infection (Delessert et al., 2005; Bu et al., 2008; Carviel et al., 2009; Wang et al., 2009; Wu et al., 2009) and all are differentially expressed following B. cinerea infection. Binding sites for NAC TFs are significantly overrepresented in a cluster of genes induced between 18 and 24 HAI with additional enrichment in clusters with similar expression profiles. The pattern of enrichment is consistent with the upregulation of many members of the NAC family around this time (Figure 9).

The evening element (Harmer et al., 2000) is overrepresented in two clusters of DEGs: one downregulated cluster (5) and one upregulated cluster (26) (see Supplemental Data Set 9 online). This element is required for evening-phased circadian regulation (Harmer and Kay, 2005) but has also been shown to play a role in the regulation of cold response genes (Mikkelsen and Thomashow, 2009). Intriguingly, 63 out of the 74 genes in cluster 5 are rhythmically expressed in mock-inoculated leaves, compared with only nine out of 134 genes in cluster 26. Perhaps cluster 5 represents the clock role of the evening element, while its overrepresentation in cluster 26 is indicative of a role in regulation of genes in response to infection.

Specific Families of TFs Are Differentially Expressed at Varying Times during B. cinerea Infection

From the TF binding motifs overrepresented in specific clusters of coexpressed genes, we obtain a pattern of TF activity. To investigate this further, we tested whether specific TF families were differentially expressed at coordinated times during the onset of infection. Using the GP2S time-local model, the times at which a gene is differentially expressed can be probabilistically identified. Using a specific threshold of P ≥ 0.5, a binary time series model for each mRNA transcript was obtained indicating whether the transcript is differentially expressed (1) or not (0) at a given time. Using these models, and family-specific gene groupings, we identified a number of TF families that were significantly overrepresented for DEGs at each time point during infection. Heat maps indicating the significance of each family’s overrepresentation are shown in Figure 9 (separated into up- and downregulated). Numeric data are in Supplemental Data Set 10 online.

A number of TF families were significantly overrepresented for upregulated genes (adjusted P < 0.05), indicating significant coordinated transcriptional activity. The earliest highly significant overrepresentation is the WRKY family, around 18 HAI. Consistent with coordinated expression of this gene family, numerous regulatory interactions between WRKY TFs have been elucidated (Eulgem and Somssich, 2007; Pandey and Somssich, 2009). Furthermore, the coordinated upregulation of WRKY TFs from 18 HAI matches the overrepresentation of W-box motifs in clusters of genes differentially expressed over this period. Shortly afterwards, the NAC family is significantly upregulated, with the timing again matching the expression profiles of clusters enriched for NAC binding motifs. The AP2/EREBP TF family is significantly upregulated at the same time and contains many genes induced by hormones and/or biotic stress (Gutterson and Reuber, 2004). As indicated above, surprisingly, the binding motif associated with this family is only enriched in one downregulated cluster. Another large TF family, C3H, shows significant overrepresentation for upregulated genes. The proteins encoded by genes in this family possess a RING-type zinc finger domain, but the family structure, membership, and function are relatively uncharacterized. Significant coordinated expression during B. cinerea infection suggests this family may be involved in plant defense responses. Coordinated expression of Trihelix TFs is a later response to infection, from 28 HAI. This family (30 members in Arabidopsis; Kaplan-Levy et al., 2012) includes GT factors that bind to GT-boxes and repress light-inducible genes. However, as with the NAC family, the function of this family is expanding to include regulation of growth and abiotic/biotic stress responses. The WRKY, NAC, AP2-EREBP, and C3H families are also overrepresented in upregulated genes during senescence (Breeze et al., 2011), suggesting interconnected roles of these TFs in response to B. cinerea and senescence.

The only TF family showing significant coordinated downregulation is the NF-YA family at 32 HAI. This family contains genes encoding A subunits of the heterotrimeric NF-Y TF complex (a trimer of A, B, and C subunits). Individual NF-Y subunits have been shown to regulate several developmental processes and tolerance to abiotic stress (for example, embryo development and flowering time; Lee et al., 2003; Wenkel et al., 2006); however, a functional trimer has only been demonstrated for NF-YA4/NF-YB3/NF-YC2, which regulates endoplasmic reticulum stress-induced genes (Liu and Howell, 2010). During senescence, NF-YA genes also show coordinated expression; however, they are overrepresented for upregulated genes (Breeze et al., 2011). The activity of the NF-Y complex could be regulated by the expression of either A, B, or C subunits; however, the coordinate regulation of NF-YA genes following both B. cinerea infection and during senescence suggests expression of A subunits may be an important control mechanism, as is the case in mammals (Manni et al., 2008).

Causal Structure Identification Network Modeling Highlights Potential Impact on Pathogen Growth

An advantage of extensive time series expression data is that it can be used in biological network inference: the prediction of the topology of a gene regulatory network. Understanding how genes interact and function together in networks to regulate plant responses is a crucial step toward being able to accurately predict the effect of genetic perturbations (i.e., phenotypic predictions from genotype). However, given the number of genes differentially expressed during B. cinerea infection, the number of time points and replicates in our data set are still insufficient to be able to generate a genome-wide network model; hence, the selection of genes to include in the model is necessary. As each cluster represents a group of coexpressed genes, we used the cluster mean as a representative of each group of genes and inferred network topology between the 44 clusters.

The expression of B. cinerea tubulin highlights two potential lag phases when pathogen growth appears to be arrested (Figure 2): 12 and 20 HAI. We compared this pathogen growth profile to the progression of transcriptional reprogramming in the host. Plotting the TOFDE obtained using the GP2S test indicates key times of transcriptional change; there is a small peak 11 HAI and a subsequent larger response beginning at 18 HAI (Figure 10). The two prospective lag phases in pathogen growth occur shortly after these peaks of transcriptional change, suggesting a possible causal relationship between transcriptional change in Arabidopsis and arrested growth of B. cinerea. To test this, we included the expression profile of B. cinerea tubulin (Figure 2), as an indicator of pathogen growth, in the network modeling. This B. cinerea growth profile was generated using the 12 observed levels of tubulin compared with PUX1 (Figure 2) and interpolating over intermediate time points using Gaussian process regression (Rasmussen and Williams, 2006).

• Download figure
• Open in new tab
• Download powerpoint

Figure 10.

The Number of Genes First Differentially Expressed at Each Time Point.

TOFDE for all 9838 DEGs following infection with B. cinerea (HAI [hpi]) is shown. The inset shows early time points (5 to 15 HAI) in more detail.

[See online article for color version of this figure.]

The discrete-time causal structure identification (CSI) algorithm of Klemm (Klemm, 2008; Penfold and Wild, 2011) was used to infer a regulatory network from the cluster means and B. cinerea growth. A section of the predicted network is shown in Figure 11, and several regulatory predictions can be made from this network. A single NAC TF, ANAC055, is present in cluster 23, and two clusters downstream of this (24 and 27) are enriched for a NAC binding motif in their gene promoter sequences. This suggests that ANAC055 regulates target genes in these clusters. Indeed, several genes differentially expressed in a knockout line of ANAC055 compared with the wild type are present in cluster 27, including ATG18a (R. Hickman, C.L. Hill, S. Ott, and V. Buchanan-Wollaston, unpublished data). Furthermore, ANAC055 acts downstream of MYC2 in JA-mediated defense against B. cinerea (Bu et al., 2008), and clusters 24, 26, 27, and 28 are all overrepresented for genes differentially expressed in a MYC2 knockout line (P < 0.01), suggesting this TF is a common upstream regulator. The activity of MYC2 is controlled by binding to JAZ proteins (Chini et al., 2009); consistent with this, MYC2 is not differentially expressed during B. cinerea infection and hence is itself not included in any cluster.

• Download figure
• Open in new tab
• Download powerpoint

Figure 11.

Inferred Network Model Using the Discrete-Time Causal Structure Identification Algorithm.

Numbered nodes represent a cluster from the SplineCluster clustering of genes differentially expressed during B. cinerea infection. The expression profile of B. cinerea tubulin was used as a proxy for pathogen growth. Selected TFs present in clusters are indicated under nodes. Colored boxes adjacent to nodes indicate motifs enriched in the promoter sequences of cluster genes, with TFs from the corresponding binding family highlighted in the same color. The color of nodes indicates the proportion of cluster genes that are differentially expressed in the tga3-2 mutant compared with the wild type (either in the 193 or 1233 set of potential target genes).

Two clusters (20 and 26) are enriched for genes containing MYB binding sites in their promoters (Figure 8) and are predicted to be downstream of clusters containing a MYB TF. Cluster 24 contains MYB2 (a known stress-associated gene that plays a role in ABA signaling; Abe et al., 2003), and cluster 32 contains MYB54. MYB54 can induce genes of secondary cell wall biosynthesis (Zhong et al., 2008). Such genes are not overrepresented in cluster 20, but MYB54 may well have additional roles in the plant. Clusters 26 and 28 are enriched for W-box motifs and predicted to be downstream of cluster 24 containing WRKY75. WRKY75 is known for its role in response to phosphate starvation but also influences basal and R-mediated defense against P. syringae (Encinas-Villarejo et al., 2009); hence, its predicted role in regulation of defense against B. cinerea is worth testing. The CSI network therefore has enabled predictions about the regulation of plant defense to be made. Although each node in the network represents a cluster of genes and hence the “causal” gene(s) in the node is not identified, specific hypotheses can be formed by integrating TF binding motif analysis.

In the CSI network model, growth of B. cinerea is upstream of at least nine clusters and is hence predicted to have a major effect on the transcriptome. Interestingly, pathogen growth is upstream of two clusters containing known clock genes, LHY and GI (clusters 6 and 34, respectively) potentially modeling the dampening of clock gene oscillations we observed in Figure 7. Only one cluster was found to be upstream of the B. cinerea growth curve, cluster 5. Cluster 5 contains two known TFs, TGA3 and ABF1. Differential expression of TGA3 starts around 18 HAI, shortly before the second pathogen lag phase, consistent with the hypothesis that downregulation of this TF may lead to the temporary arrest of pathogen growth. Therefore, we tested the effect of TGA3 expression on susceptibility to B. cinerea using two independent knockout lines, tga3-2 and tga3-3. Both mutant lines showed altered immunity, but surprisingly both knockout lines showed increased susceptibility to B. cinerea infection (Figure 12). However, the prediction of TGA3 expression influencing B. cinerea growth from network modeling has led to identification of a new player in the defense response against this pathogen. More generally, the inclusion of phenotypic information into network inference may allow for important insights that could otherwise be missed.

• Download figure
• Open in new tab
• Download powerpoint

Figure 12.

Susceptibility of tga3 Mutants to B. cinerea.

Altered susceptibility of tga3-2 and tga3-3 T-DNA insertion lines compared with Col-0 control. Lesion perimeters are a mean of 20 leaves drop inoculated with 10 μL of 10⁵ spores/mL of B. cinerea. Significantly different lesion perimeters in insertion lines compared with their respective Col-0 control were determined using a two-tailed Student’s t test assuming equal variance. **P value < 0.001 and ***P value < 0.0001. Error bars show se. hpi, hours post inoculation.

Identification of TGA3-Regulated Genes during B. cinerea Infection

Having demonstrated that TGA3 expression influences susceptibility to B. cinerea, we wanted to determine how TGA3 exerts its effect on defense by identifying the regulatory targets of this TF. The transcriptome of tga3-2–infected leaves was compared with tga3-2 mock-inoculated controls and Columbia-0 (Col-0) infected leaves at three time points (16, 24, and 32 HAI). Genes differentially expressed between the wild-type and tga3-2 infected leaves and between tga3-2 mock-inoculated and infected leaves at one or more time points were compared with the 9838 DEGs from the time series (see Supplemental Figure 5 and Supplemental Data Set 11 online). A total of 2479 genes were differentially expressed between tga3-2 and Col-0–infected leaves; 1426 of these are also differentially expressed in the time series infection data and hence are most likely to include target genes of TGA3 reproducibly regulated during B. cinerea infection. These 1426 genes can be partitioned into a group of 193 that are not differentially expressed between tga3-2–infected and mock-inoculated leaves (i.e., regulation during infection appears to be totally dependent on TGA3) and 1233 genes that are still differentially expressed between tga3-2–infected and mock-inoculated leaves. This latter group represents additional potential target genes of TGA3 with altered, but not abolished, expression in the knockout lines. As expected, TGA3 fell into the group of 193. Looking at where tga3-2 DEGs are in our network model (Figure 11), it is apparent that the 193 genes whose expression is totally dependent on TGA3 are more prevalent in groups of genes predicted to be downstream of cluster 5 (containing TGA3) or cluster 5 itself, than in cluster 23 (containing ANAC055) or groups predicted to be downstream of this cluster only.

The 1426 genes differentially expressed in tga3-2 knockouts compared with the wild type during B. cinerea infection will contain both direct and indirect targets of TGA3. Screening upstream promoter sequences of these potential targets identified 395 genes that had one or more exact matches to the consensus binding sequence TGACGT in their promoters (see Supplemental Data Set 12 online). These represent the most likely direct targets of TGA3 and include PR1, a known target of TGA3 (Johnson et al., 2003). A few likely direct targets are worth noting and are shown in Figure 11. Like PR1, expression of WRKY70 is reduced in the tga3-2 mutant compared with the wild type after B. cinerea infection. WRKY70 has a known role in defense and appears to integrate SA and JA responses (Li et al., 2006; Knoth et al., 2007). Two genes involved in the camalexin biosynthetic pathway (TSB2 and CYP79B2) are both potential direct targets of TGA3 and are coexpressed during B. cinerea infection (Figure 11), suggesting a role for TGA3 in regulating this pathway. Two other potential direct targets of TGA3 (BAK1 and BRL3) are highly upregulated during B. cinerea infection and may functionally interact. BAK1 is a Leu-rich repeat receptor kinase originally isolated through its role in brassinosteroid responses (dimerization with the Leu-rich repeat receptor kinase BRI1) (Li et al., 2002). BAK1 also dimerizes with FLS2, a pattern recognition receptor in plants, and absence of BAK1 increases susceptibility to biotrophic and hemibiotrophic pathogens (Roux et al., 2011). BRL3 is similar to BRI1 and has been shown to also bind brassinosteroids (Caño-Delgado et al., 2004). The positive regulator of ABA responses, SNRK2.3, is upregulated in the tga3-2 mutants compared with the wild type and contains a motif matching the TGA consensus within its upstream sequence. As with other ABA signaling components, SNRK2.3 is downregulated during B. cinerea infection, and TGA3-mediated expression appears to be one mechanism for this. Interestingly, this target gene analysis suggests two potential regulatory mechanisms of TGA3 itself. At a protein level, TGA3 interacts with NPR3; NPR3 is downregulated in tga3-2 mutants compared with the wild type and its promoter contains a TGA binding site. Similarly, the promoter of TGA3 contains a TGA binding site, suggesting autoregulation or regulation by other TGA factors. Despite including clear defense-related genes, the TGA3 target genes are enriched for annotation with the GO terms “response to abiotic stimulus,” “response to oxidative stress,” and “response to water deprivation,” suggesting TGA3 may play a wider role in plant responses to stress.