Differentiation of Arabidopsis guard cells: analysis of the networks incorporating the basic helix-loop helix transcription factor, FAMA

: Nearly all extant land plants possess stomata, the epidermal structures that mediate gas exchange between the plant and the environment. The developmental pathways, cell division patterns and molecules employed in the generation of these structures are simple examples of processes used in many developmental contexts. One specific module is a set of “master regulator” basic helix-loop helix (bHLH) transcription factors that regulate individual consecutive steps in stomatal development. Here we profile transcriptional changes in response to inducible expression of FAMA, a bHLH whose actions during the final stage in stomatal development regulate both cell divisions and cell fate. Genes identified by microarray and candidate approaches were then further analyzed to test specific hypothesis about the activity of FAMA, the shape of its regulatory network and to create a new set of stomatal-specific or stomatal-enriched reporters. activator and repressor. Question mark indicates uncertainties concerning the exact nature of the regulatory feedback loops detected in our assays. Black lines represent trends that are backed up by ChIP data and/or reporter line analysis while grey lines represent plausible but yet-to-be-demonstrated regulatory steps.


Introduction
Stomata are pores in the plant epidermis that serve as the main conduits for gas exchange between the plant and the atmosphere. Stomata are found in virtually all land plants and their activity, when considered collectively, plays a pivotal role in global carbon and water cycles (Hetherington and Woodward, 2003). Individual plants use these pores to efficiently balance photosynthesis with water availability. They do so in the short term by regulating the aperture of existing stomata, but in the long term, rely on developmental processes to alter the number of stomata on emerging leaf surfaces.
Each step of stomatal development, from initiation to terminal differentiation, is highly organized and regulated (reviewed in Bergmann and Sack, 2007;Nadeau, 2009;Abrash and Bergmann, 2009;Peterson, et al., 2010). In dicots like Arabidopsis, stomatal formation begins with asymmetric divisions of protodermal cells to generate dispersed meristemoids. Meristemoids may continue dividing asymmetrically several times before differentiating into a guard mother cell (GMC), the immediate precursor of the paired guard cells of the stomatal complex. Grasses, in contrast, create GMCs by asymmetric divisions in distinct cell files. The GMC then recruits subsidiary cells from the neighboring cell files before completing its division into a pair of guard cells (Liu, et al., 2009;Hernandez, et al., 1999). As seen in these two examples, diversity and flexibility are the hallmarks of the early stages in stomatal development, but later these trajectories converge on a common final step wherein the GMC undergoes a single symmetric division, and the resultant daughters undergo significant shape and gene expression changes to become the two functional guard cells of the mature stoma ( Figure 1A). Not only is the transition from GMC to guard cells the most conserved developmental process in stomatal formation, but current data indicate that the underlying genetic control of the transition is also conserved (Liu, et al., 2009).
Genes regulating stomatal development are best known in Arabidopsis. Here, a combination of cellcell signaling mediated through ligands and cell surface receptors is coordinated with (and with SCRM/ICE1 and SCRM2 in the yeast two-hybrid assay and in bimolecular fluorescent complementation (BiFC) assays in plants (Kanaoka, et al., 2008).
Stomata are functionally important, the transition between division competence and terminal differentiation is a widely faced regulatory problem and bHLHs are a large fraction of the plant transcription factor repertoire. Analysis of FAMA in the GMC to guard cell transition, therefore, has the potential to illuminate key questions of broad interest. Previous data have identified a few potential partners and targets of FAMA (Ohashi-Ito and Bergmann, 2006;Kanaoka, et al., 2008) and high-throughput studies in Arabidopsis and Brassica have significantly advanced our knowledge of the mature guard cell transcriptome, translatome and proteomes (Mustroph, et al., 2009;Yang et al., 2008;Zhang et al., 2008). How guard cells are made, and how FAMA directs this transition, however, is not well defined. In this study we identified FAMA targets and regulators of guard cell development by microarray-based transcript profiling of genes differentially modulated during a timecourse of FAMA induction. We further identified a subset of these genes as potential direct targets of FAMA through expression and chromatin immunoprecipitation (ChIP) analysis. Our results from these targets lead to specific hypotheses about the activity of FAMA and the shape of its regulatory network.

Generation of transcriptional profile in response to FAMA induction
FAMA was previously shown to induce morphologies and gene expression patterns typical of guard cells in epidermal cells and, to a lesser extent, in mesophyll cells throughout the plant (Ohashi-Ito and Bergmann, 2006). Because of this ability to drive multiple cell types into guard cell identity, we used this previously described estrogen inducible expression system in a timecourse experiment to identify genes that were differentially regulated by FAMA and/or by the acquisition of guard cell identity (details in methods). Two time points were analyzed, 4 hours (h) and 48 h post-induction; these time points were chosen to enrich for potential direct targets of FAMA (4h) and to capture the set of genes expressed when guard cells differentiate (48h). Five-day-old seedlings were used as the starting material for induction. No morphological changes were observed after induction of FAMA expression for 4h, but examination of these plants several days after induction treatment revealed some ectopic stomata, suggesting that this short induction was sufficient to induce FAMA activity (not shown). Subsequent monitoring of FAMA expression confirmed strong induction of expression after 4h, (Figure 2, cluster I). Three biological replicates for each sample and control at each timepoint were collected. RNA was isolated from whole seedlings at these timepoints, processed and labeled using standard protocols and hybridized to Affymetrix ATH1 arrays (see methods).

Analysis of microarray data
We identified genes modulated by FAMA expression in the microarray data using two-way analysis of variance and a Benjamini and Hochberg multiple testing correction with a P<0.001 and three different fold-cutoff (FC) values (2, 1.5 and 1.2, Supplemental Figure 1). These FCs yielded gene lists containing 942, 1969 and 2986 genes modulated 4 and/or 48h after FAMA induction, respectively. After examining the resultant gene lists for "known" stomata regulators, we chose a FC value of 1.2 (Supplemental Table 1) for subsequent analyses. To identify dominant trends in the data, this filtered dataset was clustered using a K-means clustering algorithm (Matrix distance: Pearson correlation) leading to the identification of 9 different clusters ( Figure 2).
After 4h of FAMA induction, 357 genes were upregulated (with FC= 1.2). A significant fraction (44%, n=164) of these genes was also upregulated at 48h but interestingly, 25% (n=89) were downregulated at the 48h timepoint ( Figure 2 clusters I,II and V; Supplemental Figure 1). Because of the highly restricted expression pattern of FAMA, this latter class is most consistent with the behavior of genes involved in promoting the GMC to guard cell developmental transition, while the genes that remain upregulated might be predicted to also have a role in guard cell identity or function. We also predicted the 4hr up classes would be enriched in direct FAMA targets; in support of this idea, we found a significant enrichment of G-boxes in the 2kb upstream regions of these early upregulated genes relative to 30067 2kb 'promoter regions" that represent the whole genome (23% as calculated by the Athena data mining tool, O'Connor, et al., 2005). This was true when we considered the whole group (110/357, 30%, p<10 -6 ) or just the 4h up and 48 down class (27/89, 30%, p<10 -4 ). These early genes encode proteins with a variety of functions. Some expected classes based on FAMA's phenotypes were transcription factors and cell cycle controllers, receptors and signaling proteins (Table 1 and Supplemental Table 2). However, there were significant numbers of genes associated with cell wall modification and metabolic processes and the most significantly enriched GO categories corresponded to catalytic (GO:0003824) and glycosyltransferase (GO:0016758; GO:0016757) activity (Supplemental Table 1).
At the 48 h timepoint, many more genes were upregulated and each to a much greater extent than at 4h. For example, the fold-change in expression in the top twenty genes from the 4h sample ranged between 4.7-2.0x, whereas the top twenty from 48h ranged between 188.7-19.6x (Figure 2 and   Supplemental Table 1). Most (1203, 88.0%) of these 48h up genes were not significantly different from controls at 4h (Figure 2, clusters I,II, III and VII). We reasoned that this later timepoint was capturing the trans-differentiation of cells throughout the plant into guard cells and therefore these "late upregulated" genes were likely to encode components present in maturing stomata. This hypothesis was confirmed by comparing these data to profiles derived from mature guard cell protoplasts (Yang, et al., 2008). 75% (1025/1367) of the genes upregulated after 48h were expressed in guard cells. This overlap was significantly different (P <0.01) than the overlap between guard cells and a larger control sample (48% 2700/5650) derived from root cells (Birnbaum, et al., 2003). The intersection of the guard cell transcriptome with genes upregulated at 4h (55%, 197/357), or downregulated at 4h (53%, 601/1126) or 48h (58%, 840/1453), was not significantly different from the control (root) population.
FAMA induction also triggered a significant downregulation of a large number of genes at 4h (n=1126, Figure 2, clusters II, III, IV, VII, VIII, IX) and 48h (n=1453, Figure 2, clusters IV, V, VI, VIII). While this trend could be due to a direct role of FAMA as a transcriptional repressor on many targets, we noticed that many of the genes in this group have tissue-specific expression in nonstomatal lineage cells. For example, this group included root specific MYB49 (At5g54230) and HYPERSENSITIVITY TO LOW PI-ELICITED PRIMARY ROOT SHORTENING 1 (HRS1, At1g13300) genes. Furthermore, 9/40 of the most highly downregulated genes after 4h were designated as root-specific by microarray expression analysis on sorted cells (Birnbaum, et al., 2003). The simplest explanation for these results is that as cells are being converted to guard cell identity, they lose their previous identities, and with them, their specific gene expression. Because of this potential developmental confound, we focused more heavily on the upregulated genes in subsequent analyses. (Nakamura, et al., 1995). KAT2 (At4g18290, 3.6 fold up at 48h) a redundantly acting homologue of KAT1 (Pilot, et al., 2001) GORK (At5g37500, upregulated 3.5 fold at 48h) a guard cell-specific major outwardly rectifying K + channel (Hosy, et al., 2003) and ABA signal transducers MPK12 (At2g46070, 8.5 fold up at 48h) and MPK9 (At3g18040, 2.5 fold up at 48 h) (Jammes, et al., 2009).

Production of reporters indicated new genes expressed in GMCs and GCs
Finding known stomata-expressed genes suggested that we would be able to use our dataset to identify new stomatal lineage enriched genes. As a first survey of expression patterns, we generated transcriptional reporters using ~2.5KB of 5' regulatory sequences to drive expression of green fluorescent protein (GFP) and/or β-glucuronidase (GUS) for 28 genes (Figure 3). These genes were chosen based on timing and degree of modulation by FAMA expression, predicted gene function, and comparison to other stomatal expression datasets described in the previous section.
Top priority genes were transcription factors, especially of the bHLH or HLH classes, given the fact that this family plays key roles in several stomatal cell fate decision points and bHLHs often participate in interconnected regulatory loops. Canonical bHLHs included bHLH039 (At3g56980, up 4.8 fold at 48h), bHLH101 (At5g04150, up 4.6 fold at 48h) and At5g08130 (up 1.4 at 4h and 2.7 at 48h). The most highly upregulated gene at 4 hr (At2g40435, 4.7 fold up, 26.6 at 48h) resembles an HLH protein and we also included this in the analysis. We made a reporter for the most highly upregulated transcription factor at 48hr, Ethylene response factor (ERF) 51 (At3g57600, Nakano, et al., 2006) and also generated tools to follow the expression of cell cycle regulators, signal transduction elements and cell wall modifying proteins. We chose the remaining genes based solely on expression profiles generated in our experiment and in comparison to other high-throughput approaches targeted toward guard cells (specific gene names and AGI codes in Figure 3 and Supplemental Table 4).
Of the 28 reporters, we could detect expression of 21 in the leaf epidermis and 19 among those were expressed in guard cells (Figure 3). Eleven genes were specifically expressed in the stomatal lineage (meristemoids, GMC and/or guard cells). The majority of these genes (7) were expressed only in guard cells, with expression peaking in mature GCs. The expression pattern of ERF51 (At3g57600) and an LRR (At1g03440), however, more closely mirrored that of FAMA. Expression of these genes began in GMCs nearing their transition to guard cells, and peaked in young guard cells, though reporter expression of these two genes persisted in maturing guard cells longer than the FAMA reporter ( Figure 3 A-C). Our reporter results indicate that induction of FAMA expression could enrich for genes expressed in the late stomatal lineage and potential FAMA target genes.
Another test of the connection between these genes and FAMA is monitoring their expression when FAMA is inactivated. We were concerned, however, that the poor viability of the fama mutant might confound results if we compared our overexpressors to these plants. We instead profiled plants expressing Est:FAMA-EAR, an inducible dominant negative construct that recapitulates the fama mutant phenotype (Ohashi-Ito and Bergmann, 2006), grown under the same conditions as the 48h samples. Many (but not all) of the genes we chose to follow more extensively were oppositely regulated in the Est:FAMA and Est:FAMA-EAR plants (Supplemental Table 3).

ChIP reveals potential direct FAMA targets
Differential expression following FAMA induction and expression in late stomatal lineage cells are characteristic of direct targets; however, to test this directly, we performed ChIP using FAMA and the 5' regions of candidate targets. To specifically pull down FAMA/DNA complexes, we created a MYC-tagged version of the estrogen-inducible FAMA construct used for the transcriptional profile and demonstrated that it was functional in inducing the formation of single-celled stomata (Supplemental Figure 2). FAMA overexpression was then induced by 10 μM β-estradiol treatment for 4h or 48h and the efficiency of induction was assayed by phenotypic analysis and a Western blot against the MYC epitope (data not shown). For each of 10 promoter regions tested, induced and noninduced samples were IP-ed with antibodies against MYC and histone H3 (as a positive control).
The abundance of the PCR product obtained from the MYC-IP with and without induction was then compared. ChIP was performed on regulatory regions of three candidate genes, CDKB1;1, FAMA and ICE1/SCRM and seven genes chosen from the microarray and reporter results, with an exon of the PIP2;1 gene (AT3G53420) used as a negative control for IP specificity. As we expected direct targets to be in the 4h induction, we focused primarily on this timepoint, however, we did also retest several of our candidates at 48h to reflect the stage at which they were most highly expressed. For these experiments, we considered a result positive only if there was no PCR amplification of the target in the non-induced sample.
We initially used ChIP analysis to test whether FAMA could directly regulate CDKB1;1 expression. CDKB1;1 is a cyclin-dependent kinase that promotes completion of stomatal lineage cell divisions (Boudolf, et al., 2004). FAMA-promoted terminal cell differentiation leads to a steep downregulation of CDKB1;1 expression level and in fama mutants, CDKB1;1 reporter expression persists ( Figure   1D, E and Ohashi-Ito and Bergmann, 2006). Here we show that this effect on expression may be direct as FAMA-MYC is associated with the ~400 base pair (bp) 5' regulatory region of CDKB1;1 ( Figure 4), but not with the negative controls.
One of the most extensively studied transcriptional regulatory paradigms was established around the animal muscle bHLH, MyoD. MyoD modulates its own expression and that of its partner bHLHs (reviewed in Berkes and Tapscott, 2005). We used ChIP analysis to test whether either of these regulatory features was also found in the FAMA and ICE1/SCRM module. Both the FAMA and ICE1/SCRM promoters contain G-box elements. We previously showed that FAMA was not absolutely required for its own expression (Ohashi-Ito and Bergmann, 2006); however these data did not rule out other autoregulatory relationships such as negative and positive feedback to fine-tune expression. ChIP assays showed that FAMA was indeed associated with its own promoter ( Figure 4); however, we did not see association with the 1000 bp 5' regulatory region of ICE1/SCRM ( Figure   4). We next used ChIP to detect in vivo interactions of FAMA with promoter regions of new candidate targets. We continued to address the relationship among bHLH and HLH factors by testing At5g08130 and At2g40435, the latter of which, while highly upregulated by FAMA, is fairly broadly expressed at low levels in the leaf epidermis ( Figure  At1g03440, though after 48h we did detect some association (data not shown). Our final choices for ChIP were expressed late in GC development, the upregulated LRR-RLK At5g53890, and downregulated SCRL17 (At3g46600). In both cases, FAMA was associated with the 5' regulatory regions (Figure 4), suggesting that it may not just set off an initial wave of gene expression, but plays a continued role in GC development.

Functional analysis of differentially expressed genes
After identifying potential target genes based on expression pattern and FAMA-ChIP results, we sought to assign functions to these genes, and to assess their contribution to a FAMA-led differentiation program. The fama phenotype of GMCs that are division-competent, but differentiation-defective, could arise from misregulation of many independent downstream programs, each of modest effect, or may be coordinated in a steep hierarchy with a few major regulators (also likely transcription factors) each controlling major subprograms. These two scenarios predict different phenotypic consequences of mutations in FAMA downstream genes.
Among the previously characterized genes we also identified in our study, a few qualitatively mirror the effects of misregulating FAMA. For example, scrm1-2 infrequently produces epidermal tumors that consist of redividing GMCs (Kanaoka et al., 2008) and elimination of CDKB1;1 and CDKB1;2 produces some single celled stomata that resemble those formed by FAMA overexpression (Xie, et al., 2010). Many of the genes highly upregulated at 48h when ectopic guard cells are present are likely to contribute more to guard cell function than development. These genes include signaling pathway components MPK9 and MPK12 (Jammes, et al., 2009) and ion channels KAT1 (Nakamura, et al., 1995) and GORK (Hosy, et al., 2003) that modulate stomatal pore aperture.
To test potential functions of the not-previously characterized genes identified in the microarray experiment, we analyzed stomatal lineage phenotypes and overall growth in 99 T-DNA insertion lines (representing 82 genes) and 17 CMV35S promoter-driven overexpression lines corresponding to a subset of those genes of interest from the array analysis (Supplemental Table 4 and 5).
Individual seedlings were screened for cell division and cytokinesis abnormalities such as arrested GMCs, malformed guard cells, stomatal clusters or altered stomatal index (number of stomata/total epidermal cells). Cotyledons and leaves of twenty 5-10 day old seedlings were analyzed per line and genotypes confirmed as noted in Supplemental Table 5. Under the screening criteria we applied, deregulated expression of no single previously uncharacterized gene produced a reproducible stomatal development phenotype (data not shown).
One explanation for lack of phenotypes is that these genes may share functions with close relatives; to address this possibility, we made higher order mutants with paralogs of the two genes expressed in GMCs, ERF51 and the LRR At1g03440. ERF51 is closely related to three other genes, ERF49 (At1g75490), ERF50 (At5g18450) and ERF52 (At2g40220). Based on microarray data, these homologues are found in guard cells, but are also expressed more broadly (e-FP browser, Toufighi, et al., 2005). LRR At1g03440' s expression pattern mirrors FAMA, but we noticed that its closest paralog, At4g03010, was also upregulated in our array (4.1 fold at 48hr) and is also stomatal lineage expressed ( Figure 3). Multiple mutant combinations of T-DNA insertions lines were constructed by crossing and genotypes confirmed by PCR, however we did not observe any stomatal development phenotype in the double LRR mutant, nor did we find reproducible, highly penetrant, stomatal phenotypes in quadruple mutant plants bearing T-DNA insertions in ERF49-52 (data not shown).

Discussion
FAMA encodes a bHLH transcription factor that acts as a "master controller" of the GMC to guard cell transition, regulating both cell divisions and cell fate. In this role, FAMA shares properties with animal bHLHs that guide muscle and neural development, with its close paralogs MUTE and SPCH involved in stomatal development and with members of other transcription factor families that regulate formative divisions ( for example, Cui, et al., 2007;Welch, et al., 2007;Willemsen, et al., 2008;Breuninger, et al., 2008;Levesque, et al., 2006). In each of these situations, the same question arises: how do these transcription factors initiate the cascade of events that culminate in precise cell identities and division behaviors?
In the experiments reported here, we used inducible expression of FAMA to identify potential targets of this transcription factor, and took advantage of FAMA's ability to redirect many cells into guard cell identity to identify genes that are enriched in this specialized cell type. We substantiated these findings by demonstrating that the expression of 11 previously uncharacterized genes was highly enriched in the stomatal lineage. We further showed that FAMA can be associated with the 5' regulatory regions of eight genes, consistent with these genes being its direct targets. Our microarray study effectively identified previously characterized stomatal-expressed proteins, many of which lead to distinctive phenotypes when misregulated (e.g. Kanaoka, et al., 2008;Nakamura, et al., 1995;Pilot, et al., 2001;Jammes, et al., 2009). We were unable, however, to identify novel stomatal phenotypes associated with any single T-DNA insertion line, or with multiple mutants of ERf51 and LRR relatives.
Recent large-scale studies of transcription factors involved in light signaling (Lee, et al., 2007), brassinosteroid response (Sun, et al., 2010) or root development (Levesque, et al., 2006) indicate that the regulatory networks downstream of major regulators are likely to be broad. Our results add to this growing observation. As an example, in the case of the cell fate regulator SHORT ROOT (SHR), Levesque et al. (2006) defined eight potential direct targets of this GRAS domain transcription factor, including one known target, SCARECROW (SCR). Later studies of some of these factors pointed to important regulatory relationships between SHR and the new targets (Welch, et al., 2007), but it is interesting to note that with the exception of SCR, loss of function mutations in these targets do not exhibit morphological phenotypes on their own (Welch, et al., 2007;Levesque, et al., 2006;Strader, et al., 2004). This theme of lack of observable phenotypes is repeated in other microarray studies on developmental regulators (Nawy, et al., 2005). Based on these previous studies with transcription factor families and our own data with the genes regulated by FAMA (and the ERF subfamily in particular), it appears that multiple overlapping functions of close relatives are the norm.
In our experiments it is also possible that the lack of observable morphological phenotypes may be because the genes we identified affect guard cell functions such as stomatal opening and closing, ion transport or alter cell wall composition. Given our primary interest in developmental regulators and the fact that functional redundancy has also been documented in these responses (e.g. channels KAT1/KAT2 (Pilot, et al., 2001) and MPK12/ MPK9 (Jammes, et al., 2009)), we did not pursue extensive physiological analyses in this study. We also cannot rule out that there are genes whose altered expression would lead to stomatal development defects in our data set and we simply did not discover them because we sampled only the currently available T-DNA insertion resources.

Refinement of FAMA activity based on its relationship to targets
As member of the H-E-R containing group of bHLHs, FAMA is predicted to bind to a G-box sequence motif; these motifs are sufficiently short, however, that they occur by chance every 4kb, and thus are predicted to be present in the regulatory regions of many genes. We previously demonstrated that the HER motif was required for FAMA activity (Ohashi-Ito and Bergmann, 2006); however, our attempts to confirm in vitro FAMA binding to G-boxsites (or to identify others by SELEX) by electromobility shift assays using purified FAMA protein have not been successful (Ohashi-Ito, unpublished), consistent with findings that FAMA poorly homodimerizes, but efficiently heterodimerizes with several bHLH partners (Ohashi-Ito and Bergmann, 2006;Kanaoka, et al., 2008). The best opportunities to find FAMA targets, therefore, were through identifying the genomic regions associated with FAMA by ChIP. In this study, we selected a small number of differentially regulated genes and found that FAMA indeed was found associated with their 5' regulatory regions. Eight of the 10 genes possess G-(or the related E-) boxes within the FAMA associated regions, but we could demonstrate association with two--CDKB1;1 and At3g49600--that do not. We also found that FAMA could interact directly both with genes that were up regulated and genes that were down regulated by its induction, indicating that FAMA could serve as a repressor as well as an activator in vivo. Such dual behavior is not uncommon; for example, WUSCHEL is primarily a repressor, but in combination with a specific partner, is part of an activating complex (Ikeda, et al., 2009).
ChIP data revealed that FAMA is able to bind to its own promoter. FAMA also binds directly to the regulatory regions of some bHLH and HLH-like proteins it induces (though not of putative partner their own expression; for example other bHLHs (eg. MyoD, Weintraub, et al., 1994) induce the expression of dimerization partners that, in addition to initiating cascades of downstream gene expression, later feedback regulate themselves. In the case of FAMA, it is particularly interesting that the most highly upregulated early gene (At2g40435) is not a bHLH, but resembles the antagonistic HLH proteins that dimerize with, and inactivate, bHLHs (Weintraub, et al., 1994).
Previously we showed that a FAMApro:GUS reporter was still expressed in the GMC tumors of a fama mutant, indicating that FAMA is not absolutely necessary to promote its own expression (Ohashi-Ito and Bergmann, 2006), however, we have noticed that the normal window of FAMA expression is very tight and even modest overexpression of FAMA can drive cells into differentiation without division, forming malfunctional 1-celled stomata (Ohashi-Ito and Bergmann, 2006). Based on these observations, an attractive model is that FAMA ensures its temporally restricted activity by both inducing antagonists and binding to its own promoter as part of negative feedback loops. FAMA's potential "active" partner ICE1/SCRM is also upregulated in response to FAMA induction, though it is not a direct FAMA target. ICE1/SCRM normally is expressed earlier than FAMA and persists in many stomatal lineage cells (Kanaoka, et al., 2008) and this higher expression level may be an indirect effect of the induction of guard cell identity by FAMA.

Implications for combinatorial control of the GMC to guard cell transition
The transition from GMC to guard cell and the subsequent maturation of the guard cells into mature, highly responsive and functional valves is an important control point and is regulated by a number of different genes in addition to FAMA. The nature of the products encoded by these genes suggests the involvement of signaling cascades (YODA and downstream MAPKs, Lampard, et al., 2009), cell cycle regulators CDKB1;1 and 1;2 (Xie, et al., 2010;Boudolf, et al., 2004), as well as other transcription factors of the bHLH (Ohashi-Ito and Bergmann, 2006;Kanaoka, et al., 2008) and MYB classes (Lai, et al., 2005;Xie, et al., 2010). We have an ongoing interest in establishing functional relationships among these different elements. Some of these factors are likely to act upstream of FAMA. For example, signaling through MAPK cascades promotes GC formation and expression of a dominant negative version of YODA phenocopies the fama mutant (Lampard, et al., 2009). We previously showed a direct connection between MAPK signaling and FAMA's paralogue SPCH (Lampard, et al., 2008), but the relationship between MAPKs and FAMA differs in two important ways. First, FAMA is not a target of MPK3 and MPK6 (Lampard, et al., 2008) Table 1).
Recent studies on FLP and MYB indicate that these proteins directly regulate the transcription of many cell cycle genes, including CDKB1;1, as they regulate the GMC to guard cell transitions (Xie, et al., 2010). Here we have demonstrated that FAMA is also associated with the CDKB1;1 regulatory region, and a vast literature supports the existence of MYB/bHLH partnerships for plant transcriptional control (Payne, et al., 2000;Zimmermann, et al., 2004;Goff, et al., 1992). Despite the common expression pattern and biological function, however, no physical interaction could be  Table 1). Identifying the precise binding sites of FAMA relative to the MYBs on the well-defined CDKB1;1 promoter may point to a new way that MYBs and bHLHs can promote specific transcriptional responses in plants.
A number of trends are revealed by the analysis of the regulators downstream of FAMA. First, there are no genes whose misregulation accounts for a large part of the FAMA phenotype, suggesting that the network downstream of FAMA is quite broad. Consistent with this, our initial ChIP analysis showed that FAMA was associated with the promoters of genes representing other transcription factors, cell cycle regulators and genes we might expect to be downstream effectors of the differentiation process, a suite of levels similar to that seen for SHR direct targets (Welch, et al., 2007;Levesque, et al., 2006). In Figure 5, we outline a working model of FAMA activity based on experimental data presented in this paper.
In the future, genome scale approaches using ChIP-seq could significantly extend the network immediately downstream of FAMA. In addition, analysis of the promoter occupancy by FAMA's potential partner ICE1/SCRM could demonstrate whether there is an obligate relationship between these two transcription factors, or whether each also participates in its own activities. The latter interpretation is suggested by the fact that ICE1/SCRM is involved in roles outside of the stomatal lineage (Chinnusamy et al., 2003) and that the penetrance of the guard cell phenotype in ice1-2/scrm null mutant plants is quite low (only 2-4 defective GMCs/leaf, DCB, unpublished). Monitoring the promoter occupancy by the new bHLH and HLH genes induced by FAMA, relative to occupancy by FAMA itself would determine whether feedforward and feedback regulation manifested by bHLH partner switching is a part of this system.

Construction of lines and choice of induction times.
We tested FAMA in several inducible expression systems available for use in Arabidopsis including ethanol (Deveaux, et al., 2003), estrogen (Zuo, et al., 2000) and glucocorticoid (Lloyd, et al., 1994;Craft, et al., 2005). After preliminary trials of all three systems, we found that the estrogen system (Estpro:FAMA) consistently gave a strong, reliable induction of the FAMA transcript in germinating seedlings without inducing severe stress responses. Two timepoints were selected: 4h and 48h based on the rationale that such timepoints should help identify either direct targets of FAMA (4 h) and/or GC expressed genes (48 h).

Microarray analysis
RNA was obtained from seedlings treated as follows: sterilized seeds of transgenic Estpro:FAMA in the Col ecotype were sown on half strength MS + 10g/l sucrose plates, then kept at 4° C for 3 days.

Seedlings were grown on vertically oriented plates in a 22° C incubator under 24 h light for 5 days.
Approximately 100 mg of seedlings/sample were then transferred with forceps to plates containing MS + Sucrose + 10 μ M β-estradiol (or new MS + sucrose plates for controls). After 4h or 48h on new plates, samples were collected (three independent replicates for each sample), frozen in liquid nitrogen and stored at -80° C. From this point on, the 3 replicates from each timepoint and genotypes were processed in parallel. RNA was isolated with the RNeasy Plant Mini Kit (QIAGEN), and samples were processed and labeled with standard protocols as provided by the manufacturer (Enzo), and hybridized to Affymetrix ATH1 array chips using protocols as described in Bergmann (2004). The array images were analyzed with the Affymetrix Microarray Suite 5.0 with the target intensity set to 500. The expression levels of the genes were analyzed using GeneSpring 4.2 software (Silicon Genetics). Genes that consistently increase or decrease in their expression levels relative to the expression level in the corresponding wildtype samples were identified using the Significance Analysis of Microarrays (SAM) software with delta value >10 and FDR <.01%. Kmeans clustering using Jexpress 9.1 (Molmine AS) was used to generate the clusters in Figure 2 with the following parameters: 9 clusters built from a maximum of 200 iterations (random seed: 647575666) using a Pearson correlation as distance matrix. Data were pre-filtered, as recommended in Orlando et al., (2009) to remove genes which were not expressed above the cutoff or that were expressed uniformly across all the measurements.   Table 5). Stomatal phenotypes were scored initially in cotyledons and first leaves of live 5-10 day old plate-grown seedlings using bright field microscopy. Confirmation of T-DNA insertions was performed by PCR using primers designed using the iSECT tools (salk.signal.edu) and standard genotyping protocols (Lukowitz et al., 2000). Double and higher order mutants were generated by crosses and PCR genotyping in subsequent generations. Supplemental table 1 summarizes the different T-DNA lines screened in this study.

Reporter and overexpression constructs and observations
Translational and transcriptional reporter lines were made by PCR amplifying a 2.5kb promoter region using D-TOPO cloning compatible primers from a genomic DNA extract (Supplemental Table 4). In most cases, a NotI restriction site was introduced on the reverse primer to facilitate subsequent cloning steps. PCR fragments were recombined into a pENTR vector using pENTR TM Directional TOPO (R) cloning Kit (Invitrogen) and further recombined into a pMDC107 (GFP) or

pMDC163 (GUS) destination vectors (Curtis and Grossniklaus, 2003). For overexpression or fusion
protein constructs, cDNA of interest were obtained by PCR amplification from total cDNA (extracted from 7day old seedlings) or directly from ABRC clones and transferred, if needed, to a gateway compatible entry clone. When necessary, 2.5 kb promoter region of the cDNA of interest was inserted upstream of the initial ATG by utilizing a Not1 site in the vector. Plant transformation vectors pMDC107, pMDC163, pMDC43, pMDC83, pMDC32 and p35GY (Kubo, et al., 2005) were used as noted in Supplemental Table 2. Every construct was verified by DNA sequencing and introduced into Agrobacterium tumefaciens strain GV3101. Arabidopsis plants were then transformed by floral dipping (Clough and Bent, 1998). GUS staining and fluorescent protein observations were performed using standard protocols as described in (Ohashi-Ito and Bergmann, 2006) and visualized on a Leica DM5000 microscope or a Leica SP5confocal microscope. For fluorescent protein fusions, cell walls were visualized by bathing seedlings in 0.02 mg/ml propidium iodide (Molecular probes P3566) for 1-5 minutes. Images were converted from Leica file format with Image J (NIH) and false colored and prepared for publication with Photoshop (Adobe) by increasing brightness and enhancing contrast. All digital enhancements were applied equally across all regions of a given image.

Accession Numbers
Microarray datasets are MIAME compliant and have been submitted to the GEO database under the GEO Series accession number GSE21786 (http://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE21786). All named genes are associated with an AGI code at first mention in text.        Genes that were upregulated at 4h and 48h compared to controls at similar timepoints. FAMA expression levels are indicated by green trace. Cluster II: Genes that were strongly upregulated at 48h, but not significantly changed at 4h relative to respective controls. Cluster III: Genes whose expression level is higher in the 48h control than in the 4h control (not of interest). Cluster IV:

Supplemental
Genes that were strongly downregulated at 4h and 48h relative to respective controls. Cluster V: Genes that were upregulated at 4h but downregulated at 48h relative to respective controls. Cluster VI: Genes whose expression level is lower in the 48h control than in the 4h control (not of interest).
Cluster VII: Genes that were downregulated at 4h but whose expression level is mainly unchanged or slightly increased at 48h relative to respective controls. Cluster VIII: Genes that were strongly downregulated at both 4h and 48h relative to respective controls. Cluster IX: Genes whose expression level exhibits a mixed response, with downregulation at 48h as a common feature (not of interest).   To form typical two-celled stomata, genes promoting guard cell identity must be upregulated while cell division processes must be coordinately downregulated during the transition from the GMCs.
FAMA expression is strong in late GMCs and early GCs. With this temporal and spatially restricted timing, FAMA may both directly activate differentiation genes while repressing cell cycle control genes such as CDKB1;1 to limit GC division. In this model, FAMA acts as both transcriptional  T  r  a  n  s  c  r  i  p  t  i  o  n  f  a  c  t  o  r  s   A  T  1  G  2  2  4 A  T  5  G  5  3  8  9  0  L  e  u  c  i  n  e  -r  i  c  h  r  e  p  e  a  t  r  e  c  e  p  t  o  r  k  i  n  a  s  e  (  L  R  R  -R  K  )  A  T  5  G  5  7  0  5  0  A  B  I  2  P  r  o  t  e  i  n  p  h  o  s  p  h  a  t  a  s  e  2  C  A  T  1  G  2  6  6  0  0  C  L  E  9  H  o  m  o  l  o  g  o  u  s  t  o  t  h  e  C  l  a  v  a  t  a  3  g  e  n  e  A  T  3  G  2  4  9  8  2  L  e  u  c  i  n  e  -r  i  c  h  r  e  p  e  a  t  f  a  m  i  l  y  p  r  o  t  e  i  n  A  T  4  G  2  9  2  4  0  L  e  u  c  i  n  e  -r  i  c  h  r  e  p  e  a  t  f  a  m  i  l  y  p  r  o  t  e  i  n  A  T  1  G  3  2  3  2 A  T  1  G  5  3  8  4  0  P  M  E  1  P  e  c  t  i  n  m  e  t  h  y  l  e  s  t  e  r  a  s  e  A  T  1  G  0  5  3  1  0  P  e  c  t  i  n  e  s  t  e  r  a  s  e  A  T  2  G  1  5  4  7  0  G  l  y  c  o  s  i  d  e  h  y  d  r  o  l  a  s  e  f  a  m  i  l  y  2  8  p  r  o  t  e  i  n  /  p  o  l  y  g  a  l  a  c  t  u  r  o  n  a  s  e  (  p  e  c  t  i  n  a  s  e  )  f  a  m  i  l  y  p  r  o  t  e  i  n  A  T  5  G  2  0  7  4  0  I  n  v  e  r  t  a  s  e  /  p  e  c  t  i  n  m  e  t  h  y  l  e  s  t  e  r  a  s  e  i  n  h  i  b  i  t  o  r  A  T  5  G  6  5  7  3  0  X  y  l  o  g  l  u  c  a  n  :  x  y  l  o  g  l  u  c  o  s  y  l  t  r  a  n  s  f  e  r  a  s  e  A  T  5  G  6  2  3  6  0  I  n  v