A local regulatory network around three NAC transcription factors in stress responses and senescence in Arabidopsis leaves.

A model is presented describing the gene regulatory network surrounding three similar NAC transcription factors that have roles in Arabidopsis leaf senescence and stress responses. ANAC019, ANAC055 and ANAC072 belong to the same clade of NAC domain genes and have overlapping expression patterns. A combination of promoter DNA/protein interactions identified using yeast 1-hybrid analysis and modelling using gene expression time course data has been applied to predict the regulatory network upstream of these genes. Similarities and divergence in regulation during a variety of stress responses are predicted by different combinations of upstream transcription factors binding and also by the modelling. Mutant analysis with potential upstream genes was used to test and confirm some of the predicted interactions. Gene expression analysis in mutants of ANAC019 and ANAC055 at different times during leaf senescence has revealed a distinctly different role for each of these genes. Yeast 1-hybrid analysis is shown to be a valuable tool that can distinguish clades of binding proteins and be used to test and quantify protein binding to predicted promoter motifs.

Introduction

As plants are sessile organisms, constant exposure to biotic and abiotic stresses has driven evolution of ever more complex mechanisms to defend against attack from other organisms and to deal with environmental perturbations. Elicitation of these stress responses requires signal perception, signal transduction and large-scale gene expression changes. Although some gene expression may be specific to each stress there is often a large overlap between stresses, suggesting that common features are likely to be involved (e.g. Kilian et al., 2007).

Genetic studies have highlighted the important role of the plant-specific NAC (NAM/ATAF1, 2/CUC2) family of transcription factors (TFs) in the regulation of stress responses (reviewed in Puranik et al., 2012). NAC proteins share a conserved N-terminal NAC domain and form one of the largest plant TF families with over 100 members encoded in the Arabidopsis genome (Ooka et al., 2003). NAC TFs have been implicated functionally in a variety of stress-related programs such as the response to drought, high-salinity, bacterial pathogens, fungal pathogens and senescence (Fujita et al., 2004; Tran et al., 2004; Guo and Gan, 2006; Hu et al., 2006; Balazadeh et al., 2010). Many members of the NAC family have overlapping expression patterns and are induced by several stresses, a situation that suggests common roles and regulation in multiple stress responses (Balazadeh et al., 2010; Nakashima et al., 2012).

Three members of the NAC family that have been implicated in overlapping stress responses are ANAC019, ANAC055 and ANAC072. These TFs form part of a closely related clade of stress-related NAC proteins of the SNAC-A class as defined recently in Nakashima et al. (2012) or group III-3 as defined in Jensen et al. (2010). There is evidence that these three genes have overlapping functions, but that they may also act individually in different stress responses. Tran et al. (2004) have shown that overexpression of any of the three genes results in increased drought tolerance, that all were induced in expression following salt stress, abscisic acid (ABA) or jasmonic acid (JA) treatment and that all three proteins could bind the promoter of the drought-inducible ERD1 gene. However, ANAC072 (RD26) has been shown to have a role in cold and desiccation stress in addition to a central role in ABA signalling (Fujita et al., 2004), while ANAC019 and ANAC055 are implicated in JA and/or ethylene signalling following pathogen infection (Bu et al., 2008).

The elucidation of the individual roles and regulatory pathways for these NACs is complicated by their overlapping activities, but here we present a dual approach to address this question and to identify common and distinguishing factors that lie upstream and downstream of ANAC019, ANAC055 and ANAC072 in stress-specific gene regulatory networks (GRNs). A key mechanism in the regulation of differential gene expression is through the sequence-specific binding of TFs to DNA motifs present in the promoters of targets. The physical interactions between these trans- and cis-acting factors and the nature of their expression-controlling activity form the basis of GRNs that tightly coordinate the spatiotemporal expression of gene products that are responsible for stress tolerance at the biochemical and physiological level. High-throughput yeast 1-hybrid (Y1H) assays are becoming an increasingly popular ‘gene-centred’ approach for delineating GRNs in Arabidopsis (Mitsuda et al., 2010; Brady et al., 2011; Castrillo et al., 2011; Ou et al., 2011). We used Y1H to identify TFs that interact with the promoters of ANAC019, ANAC055 and ANAC072 and, by modelling expression patterns of these potential regulators, excellent candidates for direct upstream regulation in different stress responses were predicted. We complemented this approach with microarray analysis of these potential upstream regulators to test the model predictions and also investigated the consequences of knockout mutants of the NAC genes to reveal genes that lie downstream in the GRN in senescence.

Results

Y1H assays identify TFs that interact with ANAC019, ANAC055 and ANAC072 promoters

The Y1H technique was used to identify TFs capable of interacting physically with the promoters of ANAC019, ANAC055 and ANAC072. Distal effects on transcription in yeast were minimised using a ‘promoter-hiking’ (Pruneda-Paz et al., 2009) approach with a series of overlapping promoter sequences of approximately 400 bp and that spanned 800–1200 bp upstream of the predicted transcription start site (TSS). These promoter fragments were used as bait against an Arabidopsis TF library (Castrillo et al., 2011) that was comprised of approximately 1500 Arabidopsis TFs. All interactions were verified pairwise with each fragment–TF pair.

The Y1H assays identified many TFs that were capable of binding to the promoters of ANAC019, ANAC055 and ANAC072 respectively: a total of 71 interactions between TFs and all bait fragments (Table 1). The interacting TFs comprised members of several TF families including bZIP, bHLH, MYB, AP2/ERF and homeodomain and are all expressed above background in green and/or stressed leaf tissues [based on internal gene expression data and/or eFP Browser analysis (Winter et al., 2007)]. In silico analysis of each fragment sequence revealed the presence of putative binding sites for all the interacting TFs identified (Figures S1, S2 and S3).

Table 1

Positive interactions identified by yeast 1-hybrid (Y1H) screens between named transcription factors (TFs) and the promoter fragments of ANAC019, ANAC055 and ANAC072. Fragments are numbered F1–F5 with F1 being closest to the transcription start site (TSS)

		ANAC019					ANAC055				ANAC072
AGI	Name	F1	F2	F3	F4	F5	F1	F2	F3	F4	F1	F2	F3	F4
AT4G34000	ABF3				•		•		•		•
AT3G19290	ABF4	•			•		•		•		•		•	•
AT2G40220	ABI4	•					•
AT1G19350	BES1													•
AT1G59640	BPEP				•
AT4G25490	CBF1										•
AT4G25470	CBF2										•
AT4G25480	CBF3										•
AT5G51990	CBF4										•
AT4G36730	GBF1				•		•		•
AT4G01120	GBF2				•		•		•
AT2G44910	HB04	•
AT5G65310	HB05			•
AT3G61890	HB12			•
AT5G66700	HB53	•
AT5G11260	HY5								•
AT3G17609	HYH								•
AT5G54680	ILR3				•				•
AT3G06490	MYB108					•			•	•				•
AT1G48000	MYB112					•	•		•	•				•
AT1G25340	MYB116					•	•		•	•				•
AT2G47190	MYB2					•	•		•	•				•
AT3G27810	MYB21					•	•		•	•				•
AT1G06160	ORA59	•						•			•
AT5G61270	PIF7	•			•			•	•		•		•
AT2G22750	bHLH				•				•
AT4G14410	bHLH104								•
AT4G34590	bZIP11								•
AT2G35530	bZIP16								•
AT1G75390	bZIP44						•		•
AT1G06850	bZIP52								•
AT3G44460	bZIP67								•

Analysis of the biological processes associated with the interacting TFs using gene ontology (GO) classifications revealed significant enrichment for TFs associated with ‘response to hormone stimulus’ and ‘response to ABA’. Expression of all three NAC TFs is induced by ABA (Tran et al., 2004) and thus targeting by TFs that are known to mediate ABA signalling would be expected. Indeed, the key ABA-related TFs ABF3 and ABF4, which recognise the ABA-responsive element (ABRE) (Yoshida et al., 2010), were found to bind to all three promoters (Table 1). ABRE occurs multiple times within each NAC promoter fragment (Figures S1, S2 and S3) and recent motif analysis revealed two ABRE-like motifs in the promoter of ANAC019 that were sufficient to drive expression in response to stress (Zou et al., 2011). The promoter fragments of ANAC019 and ANAC055 also interact with another key mediator of ABA expression, ABI4 (Table 1), an AP2 family TF whose predicted binding motif (Niu et al., 2002) is present in the promoters of these genes (Figures S1 and S2).

Members of a MYB TF family group, MYB2, MYB21, MYB108, MYB112 and MYB116, bound to the promoters of all three NAC genes (Table 1). An example of these interactions is shown in Figure 1(a) and demonstrates the interaction of overlapping fragments of the ANAC055 promoter (ANAC055_F3 and ANAC055_F4) with these five MYB TFs. Consistent with this observation is the presence of a MYB binding motif (TRANSFAC ID: M00218, Solano et al., 1995) in the overlapping region of these two fragments (Figure 1b). The five MYB TFs form part of a phylogenetic clade within the MYB family, based on the similarity of the DNA- binding domain (Stracke et al., 2001). Interestingly MYB21, the least related within this group, shows a weaker interaction (i.e. less growth) with the NAC promoters (Figure 1a). Members of this MYB subgroup are implicated in regulating the same stress responses as the NAC TFs. MYB2 is involved in the ABA-mediated regulation of salt and drought-responsive genes (Abe et al., 2003), whilst MYB108 is induced in response to ABA, JA and ethylene and is involved in regulating the response to infection by B. cinerea and other pathogens (Mengiste et al., 2003; Mandaokar and Browse, 2009). The observation that only these related members of the MYB family demonstrate positive interactions with the three NAC promoters suggests that this binding occurs in a selective, sequence-specific manner.

Figure 1

Related MYB transcription factors (TFs) bind to the promoter of ANAC055.

(a) Interactions between MYB TFs and ANAC055 promoter fragments 3 and 4 linked to the HIS3 reporter gene were transformed into yeast with plasmids carrying an Activation Domain (AD)-MYBTF translational fusion or the control plasmid expressing AD alone. Growth was assessed in the presence of increasing 3-amino-1,2,4-triazol (3AT) concentrations.

(b) The ANAC055 promoter indicating the putative MYB binding site highlighted in bold (TRANSFAC ID: M00218, Solano et al., 1995) in the overlapping region of fragments 3 and 4.

In contrast to the common interactions of ABA-related and MYB family TFs with all three NAC promoters, binding of four highly related proteins, CBF1 (C-repeat/dehydration responsive element (DRE) Binding Factor 1) CBF2, CBF3 and CBF4, was detected to the ANAC072 promoter only (Figure 2a). The CBFs bind via the DRE, which has the consensus (A/G)CCGAC (Sakuma et al., 2002) and is present in the interacting promoter fragment (Figure 2b). To confirm that CBF1–CBF4 are interacting via the DRE motif, mutational analysis was performed and the protein/DNA interactions were quantified. Firstly the DRE motif was completely replaced with the sequence TTTTT (Figure 3a – M1). This replacement resulted in the abolition of the interaction between the ANAC072 promoter fragment 1 and all four of the CBFs (Figure 3c). Hao et al. (2002) showed that the core DRE was CCGAC by testing the effect on CBF1 binding of mutating the nucleotides at either end of this motif. As this study identified only the boundaries of the DRE motif, mutation of the central nucleotide (CCGAC) was tested (Figure 3a – M2) and was shown to eliminate the binding of the CBFs (Figure 3b,c). Furthermore, the DRE consensus includes a sixth residue at the 5′ end of this core motif, an A or a G. To investigate the importance of these alternative residues in CBF binding the ANAC072 motif was mutated to ACCGAC (Figure 3a – M3). This mutated motif still retained some ability to bind all four CBFs but all are compromised, with the interaction with CBF1 and CBF2 being more severely affected than that with CBF3 and CBF4 (Figure 3c).

Figure 2

The C-repeat/dehydration-responsive element (DRE) binding factor genes (CBF1–CBF4) bind to the promoter of ANAC072.

(a) Interactions between CBF1-CBF4 and fragment 1 of the ANAC072 promoter (ANAC072_F1) were confirmed by transforming yeast that carried the HIS3 reporter gene under the control of ANAC072_F1 with plasmids that expressed the AD-CBF fusion proteins or a control plasmid that expressed AD alone and by assessing growth on media that lacked histidine.

(b) The DRE motif (consensus (A/G)CCGAC, bold text) in the ANAC072 promoter on the reverse strand reading 5′ to 3′.

Figure 3

The CBF family binds to the ANAC072 promoter via the dehydration-responsive element (DRE) motif.

(a) The DRE motif is shown underlined within the ANAC072 promoter sequence. Note that the sequence is the reverse strand reading 5′ to 3′ relative to the direction of transcription, to allow comparison of sequence to the DRE consensus. Mutated bases in the three mutants are shown in bold, lower case.

(b) Yeast that carried the HIS3 reporter gene under the control of each mutant promoter were transformed with either pDEST22 alone or each of CBF1-CBF4 and the resulting strains were examined for growth on media that lacked histidine. Example plates are shown for the mutant promoters plus either pDEST22 or CBF1.

(c) Quantification of the growth on the plates described in (b). Relative density of the growth at each cell dilution was measured using ImageJ.

The absence of interactions between the CBFs and the promoters of ANAC019 and ANAC055 suggests that a different control mechanism regulates ANAC072. Representation of the interaction data in Table 1 in the form of a network (Figure 4a) hints at other distinct mechanisms, with the bZIPs binding only to ANAC055 and the HB TFs binding solely to ANAC019. Taken together, this Y1H network shows clear commonalities and distinctions between the TFs that are capable of binding to the NAC genes in question. However, whilst Y1H analysis identifies potential regulators, showing which proteins are capable of binding, this technique does not predict the in vivo conditions under which these TFs may bind.

Figure 4

Hierarchical causal structure identification (hCSI) modelling predicts stress-specific sub-networks using different stress timecourse datasets based upon the core yeast 1-hybrid (Y1H) network (a). Networks for (b) B. cinerea infection and (c) developmental senescence are shown. True Y1H connections predicted by CSI are shown as solid black lines in which the thickness of the line represents the likelihood of that interaction, and the arrows show direction.

Network inference identifies context-dependent sub-networks

The availability of time series gene expression datasets enables the use of modelling approaches to predict interactions likely to occur under various stress treatments, although it should be noted that the action of TFs that are regulated through post-translational modifications rather than transcriptionally will not be identified using such methods. The hierarchical causal structure identification algorithm (hCSI; Penfold et al., 2012) was used to infer a separate network structure for each stress/condition using time series gene expression datasets from Arabidopsis leaves – developmental senescence (Breeze et al., 2011), Botrytis cinerea infection (Windram et al., 2012), osmotic, cold and salt stress (Kilian et al., 2007) – whilst jointly constraining the topology of the networks to favour similar structures. This approach allows for the possibility that some regulators may be universal amongst the stress response, whilst others may be condition specific.

The inferred networks for the Arabidopsis response to B. cinerea and during developmental senescence are shown in Figure 4(b,c) respectively, in which the thickness of the lines represents the marginal probabilities of a connection in that stress. Further networks that predict connections during cold, osmotic and salt stresses are shown in Figure S4. Clear differences are obvious, for example in the differential stress-dependent binding of the four CBF family members to the promoter of ANAC072. CBF1–CBF3 are predicted to bind to and regulate ANAC072 under conditions of cold stress, with the likelihood of CBF4 involvement being less significant (Figure S4). CBF1 and CBF2 appear to be important in the senescence response (Figure 4c) whereas CBF4 is predicted to regulate ANAC072 expression during B. cinerea infection (Figure 4b). In contrast, CBF3 is predicted to have the most influence in cold, osmotic and salt stress (Figure S4).

The hCSI analysis indicates that, of the five MYB proteins shown to bind the three NAC promoters in Y1H, MYB116 appears to have no significant role in developmental senescence or the stress responses modelled here (Figures 4 and S4). Interestingly the predicted regulatory role of MYB2 is different between senescence and B. cinerea. Although predicted to influence all three NACs, the models suggest that MYB2 is highly likely to regulate ANAC019 and ANAC055 in B. cinerea infection (Figure 4b), whilst in developmental senescence its main role is predicted to be in the regulation of ANAC072 (Figure 4c). MYB108 is predicted to influence the expression of all three NACs in all stresses modelled, particularly ANAC055 in B. cinerea infection (Figure 4b) and salt stress (Figure S4a). These predictions were tested experimentally as described in the next section.

Expression of ANAC019, ANAC055 and ANAC072 is perturbed in myb2 and myb108 mutants in different conditions

As shown above, a related group of five MYB TFs can bind to the promoters of ANAC019, ANAC055 and ANAC072 in Y1H, with a predicted stress-dependent differential contribution to the control of the expression of the three NACs. Following B. cinerea infection, MYB108 was predicted to regulate the expression of all three NAC genes, but with a more significant influence on ANAC055 (Figure 4b); and MYB2 was predicted to influence the expression of ANAC019 and ANAC055 and, to a lesser extent, ANAC072 (Figure 4b). The hCSI algorithm does not give the sign (activation or inhibition) for each interaction. However, through comparison of MYB and NAC expression profiles in B. cinerea infection and developmental senescence (Figure S5; Windram et al., 2012 and Breeze et al., 2011 respectively) it is predicted that the action of the MYBs would be to activate NAC expression. To test these predictions the expression of the three NAC genes was determined in leaves of myb108 and myb2 T-DNA insertion mutants at two time points following infection with B. cinerea (Figure 5a). None of the NACs shows a significant expression change in myb108 compared with wild type (WT) at 26 h post infection (hpi) however, at 30 hpi all genes show reduced levels of upregulation as was predicted in our model, with ANAC055 showing the greatest reduction in expression. In the myb2 mutant, as with myb108, there is no differential expression of the three NAC genes at the earlier time point (24 hpi), but by 30 hpi the expression of ANAC019 is statistically significant lower whereas ANAC072 expression is not altered significantly. However, the predicted effect on ANAC055 expression is not observed.

Figure 5

Expression levels of ANAC019, ANAC055 and ANAC072 in myb2 and myb108 mutant backgrounds under different stress conditions. Heatmap representation (P-values) of expression changes of ANAC019, ANAC055 and ANAC072 shows downregulation in response to infection by B. cinerea (hours post infection; hpi) (a, b) and dark-induced senescence (DIS) at two time points (T1, T2) (c, d). Numbers shown on heat map are ratios of gene expression in mutant compared with wild-type (WT). (e) Colour scale for P-value significance.

To test the predictions of the senescence model, the expression of the three NAC genes was determined in leaves of the myb108 and myb2 mutants subjected to dark-induced senescence (DIS) at two time points termed DIS1 and DIS2 (8 and 9 days after darkness respectively). In the senescence hCSI model, MYB108 is predicted to have an effect on the expression of all three NACs. However, in the myb108 mutant, at both time points only ANAC019 and ANAC055 show a significant reduced level of expression (Figure 5c). The expression of ANAC072 is not significantly altered in the mutant. In contrast, MYB2 is predicted to significantly affect the expression of ANAC072, whilst having less influence on the expression of ANAC019 and ANAC055. In the myb2 mutant expression of both ANAC019 and ANAC055 is significantly reduced at DIS2 but ANAC072 is unaffected (Figure 5d). The discrepancy between the predicted and experimental results with respect to expression of ANAC072 indicates that additional components are involved in the regulation of ANAC072 during senescence. The data points in the senescence time course used for the CSI modelling are taken at 2-day intervals. In contrast, at 2 hourly intervals the data points in the B. cinerea time course are more highly resolved thus providing a more detailed picture of the immediate consequences of temporal changes in gene expression. Additionally, the CSI modelling used data from a developmental senescence experiment. In contrast, analysis of the myb2 and myb108 mutants was performed following DIS.

The use of expression data to predict the most likely interactions from the protein DNA-binding studies generated testable stress-specific models. In this example, we tested the predicted interaction of two of the many upstream interactors in two different stresses. Most of the predictions tested were confirmed, a finding that indicated that many of the additional interactions shown in Figures 4 and S4 are also likely to be true.

Microarray analysis identifies genes and pathways regulated by ANAC019 and ANAC055 during senescence

To identify downstream components in the NAC GRN, microarray analysis of T-DNA mutants anac019 and anac055 was performed at five time points (TP1 to TP5) during developmental leaf senescence. Firstly genes that were differentially expressed over the whole time course were identified using a Gaussian process two-sample test (GP2S) (Stegle et al., 2010). Compared with WT, 2785 and 7457 genes demonstrated altered expression during the time course in the anac019 and anac055 mutants respectively. To identify genes that were differentially expressed from WT at each time point, rather than overall, a t-test analysis was performed with the data at each time point (Tables S1 and S2). This double selection considerably refined the number of differentially expressed genes for GO term analysis.

Table 2

Enriched gene ontology (GO) terms in groups of genes showing higher or lower expression in the NAC gene knockout mutants compared with wild-type (WT) at different times during senescence

	anac019 mutant		anac055 mutant
Days after sowing	Up-regulated	Down-regulated	Up-regulated	Down-regulated
23	Response to stress JA biosynthesis Response to JA		Response to SA Response to timulus
29			Response to stimulus Cell wall	Response to chitin JA stimulus
31	Indole derivative biosynthesis Defence response Photosynthesis	Flavonoid biosynthesis SA signalling	Response to stimulus Transporter activity	Response to stress JA stimulus
33	Chromatin assembly Response to biotic stimulus Photosystem	Flavonoid biosynthesis SA signalling	Response to chitin Defence response	Chloroplast Response to stress
35	JA biosynthesis Response to stimulus Response to JA		Response to chitin Response to stress	Response to stimulus Response to stress

JA, jasmonic acid; SA, salicylic acid.

Gene ontology (GO) term analysis for groups of up- or downregulated genes revealed that there was little overlap in pathways affected by mutations in the two NAC genes, indicating that they have clearly different roles in the senescence response (Table 2 and S1). Expression of both NAC genes was significantly downregulated at every time point in the relevant T-DNA mutant, but was not affected by the other mutation and indicated that neither gene is dependent on the other for expression (Figure 6a). Several photosynthesis genes, including LHCA4 (Figure S6a), LHCB6, PSBX and PSAK, showed increased expression in anac019 compared with WT – indicating a delay in aspects of the senescence process. Conversely, certain chloroplast-related genes, such as SBPase, are downregulated in anac055 – possibly indicating a more advanced decay and accelerated senescence in this mutant (Figure S6a).

Figure 6

Expression patterns of selected differentially expressed genes in the NAC gene mutants.

(a) Expression of the target NAC genes ANAC019 and ANAC055 in each knockout (KO) mutant: (i) in anac019 KO; and (ii) in anac055 KO.

(b) Expression of chitin response (e.g.WRKY53) and jasmonic acid (JA) signalling genes (e.g. JAZ10) showing altered expression in the anac055 KO mutant.

(c) Expression of genes involved in flavonoid biosynthesis (e.g. DFR and TT8) showing delay in expression in the anac019 mutant. The solid line shows expression of the gene indicated in the mutant, a dashed line shows the wild-type (WT) expression.

The GO terms ‘response to stress’, ‘response to stimulus’ and ‘defence response’ are over-represented in genes upregulated in the anac019 mutant. The most consistent differentially expressed stress response gene is SAT32, a salt responsive gene, that is upregulated at all five time points (Figure S6b). Overexpression of this gene in Arabidopsis increases tolerance to salt (Park et al., 2009). Other upregulated stress-related genes include PHYTOALEXIN DEFICIENT 3 (PAD3), PR4, PDF1.2, ACCELERATED CELL DEATH 6 (ACD6) and RECEPTOR LIKE PROTEIN 46 (AtRLP46), all of which are involved in the defence response and could improve stress tolerance (Table S1 and Figure S6b). In contrast, the anac055 mutant demonstrates reduced or delayed expression of several abiotic stress response genes, including LEA14, SAT32, RHA2A, EDS5 and ERD5 (Table S2 and Figure S6c).

Gene ontology (GO) term analysis indicates that these two closely related NACs may have reciprocal roles in the regulation of JA and salicylic acid (SA) signalling. ANAC055 may be required for normal JA signalling, while ANAC019 could enhance SA and repress JA signalling. Enriched GO terms for JA biosynthesis and signalling at TP1 and 5 in the anac019 mutant are illustrated by increased expression of JA biosynthesis genes, including LOX2 and ALLENE OXIDE SYNTHASE (AOS) (Figure S6b), and JA response genes including PR4 and PDF1.2 (Table S1). In contrast, reduced expression of JA signalling genes such as JAZ7 and JAZ10 is observed in the anac055 mutant (Figures 6b and S6c and Table S2). SA signalling appears to be the dominant pathway in the anac055 mutant with GO term enrichment for this response associated with upregulated genes at TP1. This is illustrated by early enhanced expression of genes such as CELL WALL-ASSOCIATED KINASE (WAK1), usually expressed in response to SA (Figure S6c) (He et al., 1999), and EDS1, which is involved in SA-mediated signalling in plant defence (Feys et al., 2001). This change may result in the apparent inhibition of the JA pathway in the anac055 mutant. Downregulation of the SA pathway in the anac019 mutant is illustrated by decreased expression of AHBP-1B (Figure S6b), a transcriptional repressor implicated in SA signalling (Fan and Dong, 2002).

ANAC019 may act to enhance expression of the flavonoid biosynthesis pathway but repress the activity of the camalexin pathway (possibly as a consequence of the repression of JA signalling). The flavonoid biosynthesis pathway is significantly downregulated in the anac019 mutant, including genes such as DFR, LDOX, F3H, TT4 and TT5, which show reduced expression at TP3 and TP4 in the mutant when compared with WT (Table S1 and Figure 6c). In addition, two potential regulatory genes MYB90 and TT8, both of which are TFs implicated in regulation of flavonoid biosynthesis (Borevitz et al., 2000; Baudry et al., 2004), are also downregulated, indicating that these TFs could be a primary target for ANAC019. Example expression patterns (for DFR and TT8) in Figure 6(c) show that the rapid induction of expression of these genes in the WT is blocked in the mutant, but that expression recovers to WT levels by TP5, indicating that the lack of this TF can be compensated for later in senescence.

Genes involved in the synthesis of the antifungal phytoalexin camalexin are associated with the GO term ‘indole biosynthesis’, over-represented in upregulated genes in the anac019 mutant. These genes include ANTHRANILATE SYNTHASE ALPHA SUBUNIT 1 (ASA1), PHOSPHORIBOSYL ANTHRANILATE TRANSFERASE 1, and INDOLE-3-GLYCEROL PHOSPHATE SYNTHASE (IGPS), which all function in the biosynthesis pathway from chorismate to tryptophan and in addition to PAD3 and CYTOCHROME P450 MONOOXYGENASE 79B2 (CYP79B2), which are required for the production of camalexin from tryptophan (Schuhegger et al., 2006). The presence of ANAC019 causes a delay in the early expression of camalexin synthesis genes (see PAD3 example in Figure S6b); by TP5 the levels of expression are the same in both mutant and WT.

In the anac055 mutant there is a striking group of genes with the GO annotation ‘response to chitin’ that show lower expression than WT at TP2, but that exhibit higher expression at TP4 and TP5 (illustrated by WRKY53 expression in Figure 6b). This group contains several TFs, including WRKY33,WRKY53,WRKY11 and ERF5, all of which are enhanced in expression in response to chitin and in defence responses (Libault et al., 2007). These genes show a peak in expression at TP2 that is considerably delayed in the mutant implying a requirement for ANAC055 for this rapid increase in expression. At subsequent time points, expression in the WT decreases while that in the mutant increases. ANAC055 may play a role in the WT induction of these genes, with other genes able to bring about the induction, although less efficiently, in the absence of ANAC055.

In summary, the microarray analysis of the two NAC mutants indicates that they have very different roles to play in regulating gene expression during developmental leaf senescence. Different downstream pathways are affected, for example the flavonoid pathway depends on ANAC019, while the pathogen response pathway induced by chitin appears to be dependent on ANAC055. Also, the genes appear to have opposite roles in regulation of the JA and SA antagonistic interaction.

Discussion

The results in this paper predict a regulatory network around three stress-related NAC genes, ANAC019, ANAC055 and ANAC072 (Figure 7). Potential upstream regulatory genes were identified by a combination of Y1H, stress-specific modelling and mutant analysis; downstream genes were predicted from microarray data.

Figure 7

Schematic showing predicted gene regulatory network around ANAC019, ANAC055 and ANAC072. Upstream genes were predicted by yeast 1-hybrid (Y1H) assays in combination with context-specific network modelling and mutant analysis. Dashed lines indicate interactions validated in senescence. Downstream components were predicted from microarray analysis of mutants anac019 and anca055 during developmental senescence.

This study exploited the use of Y1H to identify interactions between promoter DNA and specific TF proteins. This method is an excellent technique to show protein/DNA binding but is prone to false-negative results as many TFs may not bind in isolation or may require binding at a fixed distance from the start of transcription. For instance, this study did not detect the interaction between the MYC2 protein and the three NAC promoters recently described by Zheng et al. (2012). MYC2 may not be made or processed properly in yeast, or may require other factors to enable binding. It has been demonstrated in yeast that transcriptional activation diminishes with increasing distance of the binding element from the TSS (Dobi and Winston, 2007), thus many true interactions may be missed. However, we identified many different interacting TFs, including members of a phylogenetic clade of the MYB TF family, which bind to all three promoters including a conserved MYB recognition site in two overlapping fragments of ANAC055 (Figure 1a). This result indicates that the chosen fragment length of approximately 400 bp does not necessarily preclude the identification of TFs binding either end of the promoter fragment.

The binding of members of this clade of MYBs shows some degree of specificity; in other experiments other MYBs within our TF library have bound different promoters in the Y1H assay. Current predicted binding motifs are oversimplified but knowledge of such sequences, in combination with Y1H, enables binding specificity to be investigated, as our mutation analysis of the DRE motif has demonstrated. This combination should allow the intricacies of sequence-specific binding to be investigated thus revealing specific TF protein binding motifs beyond the simple gene family motifs we currently employ.

Many TFs occur in large families sharing a similar DNA-binding domain, including NACs, bZIPs and homeodomain TFs (Riechmann et al., 2000). Promoter evolution has been suggested to drive functional differences between members of several stress-related TF families. For example, the CBF TF family, comprised of CBF1, 2, 3 and 4, are important for regulating responses to drought and cold stress. CBF1, 2 and 3 are induced by low temperatures but not dehydration or ABA (Gilmour et al., 1998; Liu et al., 1998; Medina et al., 1999) while CBF4 is induced by dehydration and ABA but not cold (Haake et al., 2002). All members have high similarity at the protein level, yet the CBF4 promoter differs considerably from those of the other CBF genes (Haake et al., 2002). The differential expression of members of TF families such as the CBFs is crucial and it is demonstrated in this paper that although all four members of this family have the ability to bind to the promoter of ANAC072 in Y1H, modelling indicates that their contribution to the regulation differs amongst this family in a context-dependent manner.

Phylogenetic analysis of the promoters of ANAC019 and ANAC055 indicates that they are extremely similar at the promoter level (Ooka et al., 2003; Tran et al., 2004) and this study demonstrates a large overlap in the binding TFs. However, there are also clear differences (Figures 4a and 7). The ANAC019 promoter is bound by a group of homeodomain TFs that show no interaction with the ANAC055 promoter, which is instead bound by a selection of bZIP proteins. Such observations suggest that promoter evolution has refined the regulation of these two NAC genes thus adding to the complexity of GRN in which they act.

A further level of complexity is seen when we consider the context in which the interactions are observed in vivo. The use of modelling algorithms allows prediction of true interactions by considering them in the context of stress-specific expression data. Such analysis indicates that although several members of the same TF family have the ability to bind to the promoters in question, they may not actually bind under all conditions in vivo, as is predicted here with the MYB and CBF TFs. This analysis also demonstrated the importance of highly resolved time series expression data, with the highly resolved B. cinerea dataset providing more accurate predictions than the senescence dataset. In some cases there may be functional redundancy between TFs, which would prevent testing of the model using knockout mutants but should allow the prediction of likely functional homologues. Additionally, it is important to consider that TFs that are required for activation of a gene do not necessarily need to be differentially expressed and would not be predicted using the hCSI algorithm as it relies on differential expression patterns.

Expression analysis of the mutants anac019 and anac055 during developmental senescence indicated involvement of these genes in different signalling pathways. Gene expression in the anac019 mutant indicates that ANAC019 may be an activator of senescence with a role in activating flavonoid and anthocyanin biosynthesis. Conversely, early downregulation of chloroplast-related genes in the anac055 mutant hints at accelerated senescence and this TF appears to be involved in the response to chitin. These genes also appear to have opposing roles in regulating the antagonistic JA and SA pathways (Figure 7). Furthermore the observation that certain genes, including WRKY33 and WRKY53, showed an apparent delay in expression in the anac055 mutant illustrates the importance of measuring the dynamic effects of a mutation.

Previous studies have described similar roles for ANAC019, ANAC055 and ANAC072 when they were constitutively and ectopically expressed, In this paper we describe the use of a combination of experimental and theoretical tools to create a network model around the three genes to identify upstream regulatory genes and downstream pathways. This analysis has illustrated common features in upstream regulators, but also a distinct set of specific interactions that may modulate the expression of each gene depending on the stress experienced. Also, analysis of pathways predicted to be downstream of ANAC019 and ANAC055 has shown that the two genes have very different roles, at least in the process of developmental senescence.

Experimental Procedures

Y1H library screen

The TF library (REGIA + REGULATORS; RR Library) (Castrillo et al. 2011) is a kind gift from the authors and comprises approximately 1500 TFs fused to an N-terminal GAL4 activation domain in pDEST22 (Invitrogen, http://www.invitrogen.com). Yeast strain AH109 (MATa – Clontech, http://www.clontech.com) was transformed with the individual TF clones as detailed by the manufacturer and 24 clones pooled per well in a 96-well plate, in two arrangements.

Gateway Conversion (Invitrogen) was performed on the pHISLEU2 vector described in Çevik et al. (2012) to generate pHISLEU2GW. Overlapping promoter fragments of approximately 400 bp were amplified in a two-step polymerase chain reaction (PCR) from Arabidopsis (Col-0) genomic DNA using KOD DNA polymerase (Merck, http://www.merckmillipore.com) for ANAC019, ANAC055 and ANAC072. Fragments were amplified with sequence-specific oligonucleotides containing half attB Gateway recombination sites (Table S3). Second round PCR was performed with generic attB oligonucleotides (Table S3). Promoter fragments were cloned into the pDonrZeo vector (Invitrogen) using BP clonase II (Invitrogen) and then recombined into pHISLEU2GW using LR clonase II (Invitrogen). Yeast strain Y187 (MATα) (Clontech) was transformed with the pHISLEU2GW-promoter clones to generate bait strains.

The pooled library and bait strains were grown in SD-Trp or SD-Leu media respectively. 3 μl of each promoter strain was spotted onto YPDA (yeast, peptone, dextrose, adenine) plates and overlaid with 3 μl of TF library pools. After incubation for 24 h at 30°C, diploid cells were replica plated onto selective plates [SD-Leu-Trp and SD-Leu-Trp-His ± 1–100 mm 3-amino-1,2,4-triazol (3AT)]. Following overnight incubation, plates were replica-cleaned, then incubated for 4 days. Growth was scored and positive colonies patched onto selective plates and grown overnight at 30°C. Colony PCR was performed by adding a colony to 20 mm NaOH, boiling for 10 min, then these were used as a template in a PCR reaction using pD22 oligonucleotides (Table S3). Products were sequenced to identify the TF showing a positive interaction.

To verify the Y1H results, Y187 was transformed with all promoter constructs and then with pDEST22 or the appropriate TF clone. Cultures were grown in SD-Leu-Trp, diluted to 108 cells/ml, 3 μl spots of serial 10-fold dilutions plated onto selective plates (SD-Leu-Trp and SD-Leu-Trp-His ± 3AT) and grown at 30°C for 3 days before scoring.

Prediction of transcription factor binding sites

Position specific scoring matrices (PSSMs) that model DNA-binding specificities for TFs isolated from the Y1H screen were retrieved from the TRANSFAC (Matys et al., 2006) and PLACE (Higo et al., 1999) databases. PSSMs for a similar TF were used when absent from the databases. The matrix similarity score (Kel et al., 2003) was computed at each position and converted to P-values based on a score distribution of that PSSM on random sequence. Motif instances that achieved a score <0.001 were judged to be candidate binding sites.

Promoter mutations and quantification of Y1H interactions

Promoter mutations were generated by inverse PCR on entry clones containing the relevant promoter sequences using oligonucleotides shown in Table S4. Entry clones were recombined with the pHISLEU2GW plasmid using LR clonase II. Serial five- fold dilutions of Y187 strains containing the promoter mutant clones and relevant TF were plated as described above. Three independent isolates of each promoter-TF pair were plated in triplicate onto selective plates (SD-Leu-Trp and SD-Leu-Trp-His ± 3AT) and grown at 30°C for 3 days before scoring. Using a consistent sized circle, the integrated density function in ImageJ was used to measure the growth of each yeast spot, normalized by subtracting the integrated density of an adjacent equal sized area of empty agar.

GO analysis

Gene ontology (GO) annotation analysis was performed using BiNGO 2.3 (Maere et al., 2005). Over-represented categories were identified using a hypergeometric test with a significance threshold of 0.05 after Benjamini–Hochberg false discovery rate (FDR) correction (Benjamini and Hochberg, 1995) with the whole annotated genome as the reference set except for the analysis of interacting TFs in the Y1H experiment in which all TFs were used as the reference set.

Causal structure identification

The Gaussian process two-sample (GP2S) approach was used to determine differential expression of each gene in the cold, osmotic and salt stress datasets from Kilian et al. (2007). GP2S was implemented as described in Windram et al. (2012), except that a log-likelihood ratio of >8 was chosen as the threshold for indicating differential expression. For the B. cinerea and senescence time series differential expression was from our previous studies (Breeze et al., 2011 and Windram et al., 2012 respectively). The hCSI approach (Klemm, 2008; Penfold and Wild, 2011; Penfold et al., 2012) was used to infer a separate network topology for the three NAC genes using data from each of five datasets described above, using the Y1H network as a constraining hypernetwork. Initial hyperparameters and prior distributions over the hyperparameters for the Gaussian process priors were set as in Penfold et al. (2012). The maximum number of TFs that could bind simultaneously within the algorithm was limited to five if the total number of putative regulators was <15 and 4 otherwise, due to the combinatorial scaling. Five Markov chain Monte Carlo chains were run in parallel, each generating 50 000 samples network structures with the first 10 000 sampled discarded to allow equilibration of the algorithm. The remaining 200 000 samples were thinned by a factor of 5 and used to calculate the marginal probability for each pairwise connection in the Y1H network.

Plant material and stress treatments

The myb2, myb108 and anac055 lines were T-DNA insertion lines Salk_045455, Salk_024059 and Salk_011069 respectively (obtained from the Nottingham Arabidopsis Seed Centre). The anac019 dSpm insertion mutant was identified with gene-specific primers in a pool of SLAT line DNA (Tissier et al., 1999). Arabidopsis plants were grown mostly as described by Windram et al. (2012). For the developmental senescence timecourse, anac019 and anac055 mutants and their WT controls, Col-5 and Col-0, were grown as described by Breeze et al. (2011); leaf 7 was tagged with cotton 18 days after sowing (DAS) and harvested from five randomly selected plants, 8 h into the light period, at 23, 29, 31, 33 or 35 DAS (full senescence).

Botrytis cinerea pepper strain spores (Denby et al., 2004) were prepared and Arabidopsis leaves treated as described in Windram et al. (2012). Col-0, myb2 and myb108 leaves were inoculated with several 10 μl droplets of B. cinerea spores. Replicate samples for the comparison between myb108 or myb2 and Col-0 were harvested at 26 and 30 hpi or 24 and 30 hpi respectively.

For the dark induces senescence screen, nine 3-week old Col-0, myb2 and myb108 rosettes, were cut and transferred to water-saturated filter paper and stored at 20°C in complete darkness. Plates were photographed daily and RGB colour values calculated for leaf 5 of each rosette using the Color Histogram function in ImageJ. RGB intensities were normalized using a white-background reference point and average red–green ratios provided a quantitative measure of leaf yellowing. A red–green ratio of around 0.8 indicates the initiation of senescence. When the average ratio of Col-0 samples was >0.8, leaf 5 for Col-0, myb2 and myb108 lines was harvested (four biological replicates). The same sampling procedure was then performed on consecutive days to sample as senescence progressed.

Microarray analysis

Total RNA was extracted from four leaves per line, labelled and hybridized to CATMA v4 arrays (Allemeersch et al., 2005; http://www.catma.org) as described (Breeze et al., 2011). For analysis of Col-0, myb2 and myb108 samples four replicates were pooled and labelled twice with each dye giving four technical replicates. Comparisons were made pairwise between WT and mutant under each condition. For analysis of the Col-0 and anac055, and Col-5 and anac019, biological replicates were labelled separately, twice with each dye, and comparisons made within and between time points in a ‘loop design’ (Kerr and Churchill, 2001). Analysis of expression differences between Col-0 and myb2 and Col-0 and myb108 under each condition was performed using the R Bioconductor package limmaGUI (Wettenhall and Smyth, 2004) applying PrintTip lowess transformation and quantile-normalization. The data were fitted to a linear model using a least squares method, P-values adjusted to control the false discovery rate (Benjamini and Hochberg, 1995). Analysis of the anac019 and anac055 time course experiment was performed using a local adaptation of the maanova package as described (Breeze et al., 2011). The GP2S approach was used (as described in Windram et al., 2012) to identify differentially expressed genes (log-likelihood ratio of >5) in the anac019 and anac055 mutants compared with WT. A t-test analysis was then performed to identify genes that were differentially expressed at each time point. Genes expressed at a higher or lower level than WT (ratio >1.7 or <0.6 respectively, P-value <0.1) were identified (Tables S1 and S2).

Data repository

The microarray data used in this paper have been deposited in NCBI's Gene Expression Omnibus (Edgar et al., 2002) and have been given a GEO Series accession number GSE46318. These data will be released on publication.