Identi ﬁ cation of Cytokinin-Responsive Genes Using Microarray Meta-Analysis and RNA-Seq in Arabidopsis 1[C][W][OA]

Cytokinins are N 6 -substituted adenine derivatives that play diverse roles in plant growth and development. We sought to de ﬁ ne a robust set of genes regulated by cytokinin as well as to query the response of genes not represented on microarrays. To this end, we performed a meta-analysis of microarray data from a variety of cytokinin-treated samples and used RNA-seq to examine cytokinin-regulated gene expression in Arabidopsis ( Arabidopsis thaliana ). Microarray meta-analysis using 13 microarray experiments combined with empirically de ﬁ ned ﬁ ltering criteria identi ﬁ ed a set of 226 genes differentially regulated by cytokinin, a subset of which has previously been validated by other methods. RNA-seq validated about 73% of the up-regulated genes identi ﬁ ed by this meta-analysis. In silico promoter analysis indicated an overrepresentation of type-B Arabidopsis response regulator binding elements, consistent with the role of type-B Arabidopsis response regulators as primary mediators of cytokinin-responsive gene expression. RNA-seq analysis identi ﬁ ed 73 cytokinin-regulated genes that were not represented on the ATH1 microarray. Representative genes were veri ﬁ ed using quantitative reverse transcription-polymerase chain reaction and NanoString analysis. Analysis of the genes identi ﬁ ed reveals a substantial effect of cytokinin on genes encoding proteins involved in secondary metabolism, particularly those acting in ﬂ avonoid and phenylpropanoid biosynthesis, as well as in the regulation of redox state of the cell, particularly a set of glutaredoxin genes. Novel splicing events were found in members of some gene families that are known to play a role in cytokinin signaling or metabolism. The genes identi ﬁ ed in this analysis represent a robust set of cytokinin-responsive genes that are useful in the analysis of cytokinin function in plants.

The phytohormone cytokinin was discovered based on its ability to stimulate cell division in cultured plant cells (Skoog and Miller, 1957). In intact plants, cytokinins play diverse roles in nearly every aspect of plant cell growth and development, including meristem function, leaf senescence, sink/source relationships, vascular development, and biotic and abiotic responses (Taniguchi et al., 2007;Séguéla et al., 2008;Argueso et al., 2009;Nishiyama et al., 2011;Argueso et al., 2012). Cytokinins regulate these processes through interactions with other plant hormones, which depend on cell type, developmental stage, and exogenous inputs (Werner and Schmülling, 2009;Wang et al., 2011b;Vanstraelen and Benková, 2012). The availability of mutants affecting distinct steps in cytokinin signaling has facilitated a general analysis of how cytokinins regulate gene expression (Kieber, 2002).
Cytokinin signal transduction involves a multistep phosphorelay system that begins with autophosphorylation of a membrane-bound receptor (ARABIDOPSIS HISTIDINE KINASE2 [AHK2], AHK3, and CYTOKININ RESPONSE1 [CRE1]/AHK4) upon cytokinin binding and subsequent transfer of the phosphate group to His phosphotransfer proteins (AHPs) and then finally to the response regulators (ARRs). There are two types of response regulators in Arabidopsis (Arabidopsis thaliana). Type-B ARRs are activated by phosphorylation and mediate cytokinin-regulated gene expression, including the up-regulation of the type-A ARRs, which act as negative feedback elements in cytokinin signaling (Brandstatter and Kieber, 1998;D'Agostino et al., 2000;Sakai et al., 2000;Imamura et al., 2003;Tajima et al., 2004;To et al., 2004;Taniguchi et al., 2007;Argyros et al., 2008). A subclade of APETALA2/Ethylene Response Factor (AP2/ERF) transcription factors, the Cytokinin Response Factors (CRFs), some of which are up-regulated in response to cytokinin, also play a role in the cytokinin-regulated gene expression (Rashotte et al., 2006).
Microarrays have been widely used in the genomewide analysis of hormone responses in Arabidopsis, including the response to altered cytokinin function (Rashotte et al., 2003;Kiba et al., 2004;Brenner et al., 2005;Lee et al., 2007;Taniguchi et al., 2007;Yokoyama et al., 2007;Argyros et al., 2008;Dello Ioio et al., 2008;Goda et al., 2008;Müller and Sheen, 2008;Argueso et al., 2010;Bishopp et al., 2011;Köllmer et al., 2011;. These data sets have identified many genes regulated by cytokinin, and several have been validated by other methods (Supplemental Table  S1). Nevertheless, there is only modest overlap in the sets of genes identified in these studies, which likely results from multiple factors, including lab/growth-specific effects, treatment-specific effects, and the identification of false positives/negatives. A complementary approach to query transcriptome changes is RNA-seq (Mortazavi et al., 2008;Chen et al., 2010;Filichkin et al., 2010;Lu et al., 2010). Despite biases related to transcript size and sequence composition, RNA-seq data represent mostly random samplings from expressed sequences and thus, in addition to determining transcript abundance, can provide new information about splicing patterns, alternative splicing, and previously unannotated genes, even in well-annotated, well-characterized genomes such as the model plant Arabidopsis. An important advantage of RNA-seq data over microarrays is that prior knowledge of the expressed transcriptome is not required. This is particularly important in the case of Arabidopsis and cytokinin signaling, because the majority of hormone-related array studies were performed using the older model ATH1 Affymetrix array, which lacks specific probe sets for more than 11,000 annotated genes. However, the statistical methods needed to analyze RNA-seq data are less well developed compared with microarray data analysis.
Meta-analysis of microarrays has been used extensively in animal systems to define robust sets of regulated genes (Nazri and Lio, 2012;Plank et al., 2012). Thus, a meta-analysis of the microarray data associated with the cytokinin responses, taking advantage of the variety of conditions employed for these analyses, would be a powerful tool to define a robust set of cytokinin-regulated genes. A recent study presented a preliminary meta-analysis of cytokinintreated microarrays in Arabidopsis, though this study was limited to comparable growth conditions, treatment, and time points, and only a small subset of the genes identified were presented .
The aim of this work was to identify a robust set of genes regulated by cytokinin, defining genes that are induced in multiple conditions using meta-analysis of the publicly available, as well as previously unpublished, microarray data sets. Filters were applied using empirically defined criteria with appropriate statistical limits. Most of the genes from this set were verified by an independent RNA-seq experiment using very NASCARRAY S-173 stringent statistical criteria. Furthermore, the RNA-seq analysis identified a set of cytokinin-regulated genes that were not present on microarrays and also provided evidence for alternative splicing events among cytokinin signaling/metabolism genes. The gene sets identified were analyzed for enrichment of Gene Ontology (GO) terms as well as using MapMan to identify processes and pathways affected by cytokinin.

Hierarchical Clustering of Expression Profiles of Experiments Used for Microarray Meta-Analysis
As the basis for our meta-analysis, we employed 13 experiments querying cytokinin-regulated gene expression changes, all of which were performed using the ATH1 Affymetrix array platform. These experiments  Table I. The expression data were gene-wise normalized, and hierarchical clustering of gene expression patterns was carried out using GeneSpring software with Euclidean distance correlation for the two series of experiments, entities and variables. The color scale at the bottom represents the relative signal intensities. included data from seven published studies as well as six previously unpublished data sets from our laboratories ( Table I). The raw data (.CEL) files from these experiments were obtained from either ArrayExpress (http://www.ebi. ac.uk/arrayexpress/) or NASCArrays databases (http:// affymetrix.arabidopsis.info/narrays/experimentbrowse.pl) or were generated in our laboratories. These microarray experiments included distinct treatment regimes, growth conditions, and/or developmental stages and were generated in multiple laboratories, resulting in potential differences in global gene expression and cytokinin responses. We incorporated all of these experiments in our analysis to identify the robustly cytokinin-regulated genes (i.e. those that are differentially regulated in multiple conditions).
To gain insight into the effect of the varying conditions/treatments on transcript expression profiles, we performed a hierarchical clustering analysis of all the transcripts in control and cytokinin-treated samples across different experiments as described in Table I. Using GeneSpring software, the expression data were gene-wise normalized (Cheng et al., 2009), and a hierarchical clustering of signal intensities of gene expression patterns was carried out using Euclidean distance correlation for two entities (the genes and the experimental conditions). This clustering revealed the relatedness of the various data sets (Fig. 1). Duplicate or triplicate microarrays clustered as nearest neighbors in all cases. The microarrays that examined root tissues clustered together, as did the experiments using etiolated seedlings. However, experiments involving shoots and whole seedlings did not clade separately from each other, perhaps reflecting the preponderance of RNA from the shoot when RNA is prepared from a whole seedling. The method of application (spraying versus immersion) did have an effect on gene expression, as these experiments did tend to cluster together. Overall, this analysis revealed substantial differences in overall gene expression profiles based on the experimental designs and tissues used for these microarrays, which suggests that genes regulated by cytokinin across these diverse conditions will reflect robustly regulated genes.

Identification of a Set of Cytokinin-Responsive Genes Using Microarray Meta-Analysis
A Fisher's combined probability test (Brown, 1975;Kippola and Santorico, 2010) was used to analyze the microarray data and identified 4,485 differentially regulated genes (P # 0.01; Supplemental Table S2). To identify genes that are robustly regulated by cytokinin, a variety of empirically defined filters were then applied to this data set. For the analysis of microarray data, 1.5-to 2-fold changes are often used as minimum cutoffs when identifying regulated genes. We explored various cutoffs values using a set of genes that were previously demonstrated to be regulated by cytokinin using additional methods beyond microarray analysis (e.g. quantitative reverse transcription [RT]-PCR and northern blotting; Supplemental Table S1), including genes involved in cytokinin signal transduction (type-A ARRs and AHK4), metabolism (CYTOKININ OXIDASEs CKX3, CKX4, and CKX5) and transcription factors. The latter were particularly useful because transcription factors are generally expressed at low levels (Holland, 2002;Vaquerizas et al., 2009), and thus their expression is difficult to accurately assess using microarrays. The relative expression of these genes in control and cytokinintreated samples as determined by each microarray was plotted with different midpoint values ( Fig. 2A). Most of these transcription factors fell below a threshold of 1.8-fold induction in this analysis, indicating that a cutoff of greater than or equal to 1.8-fold excluded many known cytokinin-regulated genes ( Fig. 2A). By contrast, a midpoint of 1.5-fold induction captured most of the known cytokinin-regulated transcription factors ( Fig. 2A).
We next examined the relationship between the number of array experiments in which a gene was upregulated greater than or equal to 1.5-fold and its identification based on the Fisher's combined probability test (P # 0.01; Fig. 2B). There is a sharp decline in the percentage of genes that passed the threshold of P # 0.01 for the Fisher's test for genes that appear in less than six out of 13 experiments (Fig. 2B). In addition, we examined a subset of genes confirmed to be regulated by cytokinin, including several transcription factors, and almost all (15 out of 16) were found to be induced by greater than or equal to 1.5-fold in at least six microarray experiments (Fig. 2C). We thus queried for all genes that displayed a greater than or equal to 1.5-fold change in response to cytokinin in at least six experiments (Table I) and that had a P value of 0.01 or less in the Fisher's test. This analysis identified 226 genes (Tables II and III), including 39 out of 55 genes previously validated by other methods in addition to microarray (Supplemental Table  S1). We refer to this list of 226 genes as the "golden list" of cytokinin-regulated genes, containing 158 upregulated genes (Table II) and 68 down-regulated genes (Table III).  Table S1) was used for this analysis. The relative expression of each gene in cytokinin versus control samples in each microarray experiment (Table I) was plotted with a midpoint value on the color scale of 1.3-, 1.5-, 1.8-, or 2-fold change as indicated. B, Examination of the number of genes that are up-regulated at least 1.5-fold in at least the number of microarray experiments noted on the x axis (blue line, left y axis), plotted along with the percentage of those genes thus identified that test as cytokinin regulated using a Fisher's combined probability test with a P value of 0.01 or less (red line, right y axis). C, Number of microarrays in which a set of 15 genes previously demonstrated to be induced by cytokinin (Supplemental Table S1) show at least a 1.5-fold induction in the microarray experiments shown in Table I. The dashed line indicates the six or more microarrays chosen as the cutoff for inclusion in the list of robustly regulated genes.
The up-regulated genes on the golden list include nine out of the 10 type-A response regulators. The sole exception, ARR17, did not have a significant meta-P value, which may reflect its low level of expression and hence poor detection on microarray experiments. Multiple cytokinin metabolism genes are also part of the golden list, including CKX4 and CKX5, which encode cytokinin oxidases/dehydrogenases (Werner et al., 2001), and a UDP-glucosyl transferase gene (UGT76C2) encoding a cytokinin O-glycosyltransferase (Hou et al., 2004). Ten cytochrome P450 genes are present in the golden list as up-regulated genes, including CYP735A2, which encodes a cytokinin transhydroxylase (Takei et al., 2004), all of which are up-regulated in response to cytokinin. The cytokininregulated CRF genes (Rashotte et al., 2006) are also part of this list, as are multiple additional AP2/ERF transcription factors. Four thioredoxin genes, which likely play a role in biotic and/or abiotic stress, are present on the golden list as up-regulated in response to cytokinin, as was a glutaredoxin gene. Other interesting genes induced by cytokinin on the golden list include two annexins, three dirigent-like proteins, an abscisic acid-dependent aquaporin (NOD26-LIKE MAJOR INTRINSIC PROTEIN1), five acyl transferases, and several genes involved in auxin signaling, including two SMALL AUXIN UP-REGULATED (SAUR)-like genes, an auxin efflux carrier, and AUXIN RESISTANT3 (AXR3). Furthermore, there were multiple genes involved in cell wall function, including expansins, pectin lyase-like, and xyloglucan endotransglucosylase/ hydrolase32. Among the 68 down-regulated genes on the golden list (Table III), there is an uncharacterized UDP-glucosyl transferase gene, UGT76B1, that is downregulated in most of the microarray experiments, suggesting a potential role in cytokinin homeostasis. Several genes encoding transporters are down-regulated, as are three peroxidase genes. The golden list also includes Table II. Up-regulated genes identified by microarray meta-analysis Genes in boldface indicate candidates validated by methods other than microarrays.

Validation of Genes Identified by the Meta-Analysis
We validated the regulation by cytokinin and characterized the response kinetics of a representative group of genes encoding transcription factors from the golden list using the NanoString nCounter system. The NanoString system allows for assessment of transcript levels with similar sensitivity to what is found with quantitative RT-PCR but does not require enzymatic reactions (e.g. reverse transcription) or amplification of targets and thus eliminates the potential biases introduced by these steps (Geiss et al., 2008;Heyl et al., 2008;Malkov et al., 2009). We chose genes encoding transcription factors for this first analysis, as they tend to be expressed at relatively low levels and thus represent a fairly stringent test of the golden list. Noncytokinin-regulated genes were included as controls, and the cytokinin-treated seedlings were compared with the dimethyl sulfoxide (DMSO)-treated seedlings at matched treatment times. All nine of these up-regulated genes from the golden list retested as being induced by cytokinin, though the kinetics of the response to cytokinin varied (Fig. 3A). AXR3 (AT1G04250) is only induced at later time points (.120 min). ZINC FINGER PROTEIN5 (ZFP5; AT1G10480) and ZFP6 (AT1G67030) showed a rapid but transient response to cytokinin, reaching almost basal level of expression in 120 min. We also examined seven genes encoding transcription factors from the golden list that were down-regulated in response to cytokinin (Fig. 3B). The transcript levels for all seven were lower than their corresponding controls.
In addition to the NanoString analysis, we validated additional genes from the golden list using quantitative RT-PCR. For this, we examined the expression of 11 up-regulated genes ( Fig. 3C) and three downregulated genes (Fig. 3D) in response to treatment with cytokinin. With the sole exception of the GIANT KILLER (GIK; AT2G35270) gene, this quantitative RT-PCR analysis confirmed the regulation of these genes in response to cytokinin. Overall, these analyses validated the quality of the golden list and support the notion that these represent a set of genes that is robustly regulated by cytokinin.

A Type-B Response Regulator Binding Motif Is Enriched in the Promoters of the Up-Regulated Genes
Previous studies have shown that type-B ARRs mediate the majority of cytokinin-regulated gene expression (Argyros et al., 2008). We thus examined the putative regulatory regions of the genes comprising the golden list to determine if they were enriched for type-B binding sequences. The core DNA sequence responsible for in vitro binding of ARR1, ARR2, ARR10, and ARR11 has been identified as AGAT (T/C; Sakai et al., 2000;Hosoda et al., 2002;Imamura et al., 2003). Further analysis revealed that the ARR1 consensus binding site includes an additional two bases on either side of the core sequence (AAGAT [T/C] TT; Taniguchi et al., 2007). While the AGAT (T/C) core motif is likely functional in planta (Ross et al., 2004), its frequent occurrence precludes its use as a diagnostic for the cytokinin responsiveness of a given promoter. This core sequence is not significantly enriched in the putative regulatory regions (defined for this purpose Table III. Down-regulated genes identified by microarray meta-analysis Genes in boldface indicate candidates validated by methods other than microarrays.

Microarrays/13
Down-Regulated Genes 11 AT2G35270 (GIK); AT3G13760 (C1 domain); AT5G43520 (C1 domain) 10 AT1G72360 (ERF73); AT4G33420 (Peroxidase); AT5G15830 (bZIP3) 9 AT2G28160 (bHLH029); AT1G77330 (ACC oxidase); AT1G78000 (SULTR1;2); AT5G43180 (Unknown); AT5G50760 (SAUR-like) 8 AT1G23390 ( Figure 3. Validation of golden list candidates by nCounter gene expression analysis and quantitative RT-PCR. A, NanoString nCounter gene expression analysis of up-regulated representative genes from the golden list. B, NanoString nCounter gene expression analysis of down-regulated representative genes from the golden list. For A and B, 10-d-old light-grown seedlings were treated with 5 mM BA or a DMSO vehicle control for 60 or 120 min as indicated, and the normalized counts for indicated transcripts were determined. Error bars indicate the mean 6 SE from three biological replicates. C, Relative transcript level of the up-regulated representative genes from the golden list as determined by quantitative RT-PCR. D, Relative transcript level of the down-regulated representative genes from the golden list as determined by quantitative RT-PCR. Ten-day-old light-grown seedlings treated with 5 mM BA or a DMSO vehicle control for 120 min were used for quantitative RT-PCR. The level of transcript relative to the DMSO control is shown as a mean of three biological replicates 6 SE of the mean. as 1 kb upstream of the translational start site) of the genes comprising the golden list (Fig. 4). By contrast, the extended version of the ARR1 binding site is substantially overrepresented in the putative regulatory regions of the cytokinin-induced genes present in the golden list ( Fig. 4) relative to the putative regulatory regions of the cytokinin nonresponsive genes (i.e. genes that displayed fold change of ,1 in response to cytokinin in all 13 microarray experiments), consistent with prior studies . As expected, there was no enrichment of the TATA box element (TTATTT) in the putative regulatory regions of the cytokinin-responsive genes and also no enrichment of the auxin response element (Aux-RE; TGTCTC) that we used as a negative control (Fig. 4). Another putative cytokinin response element (TATTAG) has been identified in cucumber (Cucumis sativus; Kuroda et al., 2000;Fusada et al., 2005), but this element is not overrepresented in the putative regulatory regions of the genes in the golden list. The enrichment of the type-B binding sites in the up-regulated genes suggests a direct regulation of these genes by type-B ARRs in response to cytokinin, consistent with the role of the type-B ARRs in cytokinin-regulated transcription (Mason et al., , 2005Argyros et al., 2008). Interestingly, there is no enrichment of the type-B element motifs in the putative regulatory regions of the cytokinin-down-regulated genes on this list, which may reflect an indirect effect of the type-B ARRs on these genes or the presence of an alternative type-B binding site not represented by the extended ARR1 binding site.

Differential Expression Analysis of RNA-Seq Data
To further explore cytokinin-regulated gene expression, we used RNA-seq to characterize the response of the transcriptome to cytokinin. We analyzed three biological replicates of 10-d-old seedlings treated for 2 h with 5 mM cytokinin or a DMSO vehicle control. The sequencing reads were aligned to the Arabidopsis genome using TopHAT, and differential expression analysis was performed with DESeq with a false discovery rate of 5% as described in "Materials and Methods." This analysis revealed 573 genes (P # 0.05) that were differentially regulated by cytokinin (423 upregulated and 150 down-regulated), including many genes previously demonstrated to be up-or downregulated under similar conditions, as well as new genes not previously identified as cytokinin regulated (Supplemental Table S4).
As expected, the type-A response regulators were among the most highly up-regulated genes, with log 2 fold changes (as calculated by DESeq) ranging from 1.68 to 5.5. Nine out of 10 type-A ARRs showed differential regulation in this experiment, with the sole exception again being ARR17. Overall, the fold induction of the type-A ARRs in RNA-seq was similar to that found in a microarray experiment that employed similar experimental conditions (Rashotte et al., 2003;Lee et al., 2007). Other known cytokinin-responsive genes appeared among the top 50 differentially upregulated genes sorted by ascending P values, including cytokinin metabolism (CKX4 and CKX5) and signaling genes (AHK4), annexins (ANNAT3 and  (Tables II and III) using the Motif Mapper of the TOUCAN 2 workbench. The number of motifs present in the promoters of the regulated genes was divided by the number of motifs present in an equal number of promoters from random genes whose level was found to be unchanged by cytokinin in all microarray experiments used in the meta-analysis. Aux-RE was chosen as a negative control.   Table S4).
Most prior studies exploring the response of the Arabidopsis transcriptome to cytokinin have used the ATH1 microarray, which was developed using the TIGRv5 annotations of the Arabidopsis genome (Redman et al., 2004). Many additional genes (approximately 11,000) have been annotated in the genome since the release of TIGRv5, and these are not represented on the ATH1 microarray. RNA-seq allows an exploration of the effects of cytokinin on these recently annotated genes. The RNA-seq analysis of cytokinin-treated seedlings revealed 73 new genes that were differentially regulated in the presence of cytokinin (60 up-regulated and 13 down-regulated; Table IV). The up-regulated genes encode proteins such as eight transcription factors (including five bHLH family members), three defensin-like proteins, two F-box proteins, a cytochrome P450, and CKX7 (Table IV). The down-regulated genes include one bHLH family member and a peroxidase gene. We used quantitative RT-PCR to independently examine the expression of the top (by fold change) 18 up-regulated and three down-regulated genes among the list of genes uniquely identified by RNA-seq (Fig. 5). With the sole exception of an uncharacterized peroxidase (AT1G77100), all of the other genes tested were confirmed to be regulated in response to cytokinin.
We identified several noncoding genes as being differentially expressed in response to treatment with benzyladenine (BA; Fig. 6A), including two microRNA genes that were up-regulated, MIR163 and MIR414. The top candidate for a target of miRNA163 identified by the AthaMap prediction tool (Bülow et al., 2012) is an S-adenosyl-L-Met-dependent methyltransferases (AT1G15125). Consistent with the up-regulation of its cognate microRNA, AT1G15125 was down-regulated in response to cytokinin (Fig. 6B). No potential targets of MIR414 were predicted using this prediction tool. The other nonprotein-coding genes (AT3G27884, AT1G31935, AT5G03668, and AT5G01215) have not been characterized. To query physiological processes potentially targeted by the transcript changes that occur in response to cytokinin, we analyzed the GO term enrichment of the cytokinin-regulated genes identified by RNA-seq using AmiGO (Carbon et al., 2009) and DAVID (Supplemental Fig. S1; Supplemental Table S5; Huang da et al., 2009). This analysis revealed not only the expected enrichment of genes involved in cytokinin signaling and metabolism, but also enrichment of genes involved in secondary metabolism (anthocyanin, flavonoid, and suberin biosynthetic pathways), nutrient assimilation (nitrate transport and iron binding), and reductase activities. An analysis of the golden list resulted in an enrichment of similar GO term categories and is also consistent with the GO term enrichment from a previous meta-analysis of cytokinin-regulated genes .
The cytokinin-regulated genes present in the GOenriched group cytokinin signaling and metabolism included almost all the type-A ARRs, the CRFs, and AHP1 (AT3G21510), as well as a cluster of cytokinin oxidases (CKX3, CKX4, CKX5, and CKX7) and a group of five N-glucosyl transferases, only one of which, AT5G05860, has been shown to act on cytokinins. The largest group of cytokinin-regulated genes (approximately Disease resistance genes were also found to have a high GO term enrichment in the cytokinin-regulated gene list. The cytokinin-regulated genes from this group included defense-associated dirigent-like proteins, defense-associated oxygenase (DOWNY MILDEW RESISTANT6), methyltransferase (AT1G77530), and chalcone isomerase like (AT5G05270). This is consistent with recent studies linking cytokinin to defense responses (Choi et al., 2010;Argueso et al., 2012). Other enriched groups of genes include GDSL family of Ser esterases/lipases, carboxypeptidases, thioredoxins, phototropic-responsive Nonphototrophic hypocotyl3 (NPH3) family proteins, ubiquitin-conjugating enzyme family members, several cytochrome P450s, and peroxidases. A number of kinase genes were identified that could link cytokinin to other signaling pathways in the plant, including many Ser/Thr kinases such as WITH NO LYSINE9, CBL-interacting kinases, Nek6, BR1-LIKE3, WAG1, CLAVATA1, several Leu-rich receptor-like protein kinases, a MAPKKK, and a senescence-induced kinase (AT2G19190).
There was an enrichment of genes related to auxin function in the set of cytokinin-regulated genes. These included a number of auxin responsive genes, including GH3 family and SAUR-like auxin responsive genes, and multiple ATP-binding cassette transporters, some of which may play a role in auxin transport. The induction of these genes involved in auxin function by cytokinin is consistent with the previously described cross talk between these two plant hormones (Moubayidin et al., 2009).
As a complementary approach to the GO term enrichment analysis, we explored the putative functions of the genes identified by microarray meta-analysis and RNA-seq experiment using MapMan (Thimm et al., 2004), which allows the visualization of cytokinininduced changes in different metabolic processes ( Fig.  7; Supplemental Fig. S3). Each gene in MapMan is initially organized in bins rather than as pathways, which allows genes to be assigned into a pathway even when their function is only tentatively known. Figure 7 shows the overall mapping of the cytokinin-regulated genes identified by microarray meta-analysis and RNA-seq, and Supplemental Figure S3 highlights relevant transcriptional changes related to overall metabolism, cellular responses, and regulation. As with the GO analysis, the cytokinin-regulated genes identified by the metaanalysis and RNA-seq results in the enrichment of similar functional categories and pathways in MapMan. As expected, there were a large number of genes assigned to the hormone metabolism bin ( Fig. 7; Table V). A large Figure 5. Validation of the cytokinin regulation of a subset of genes identified by RNA-seq that is not present on the ATH1 microarrays. A, Relative transcript level of the up-regulated indicated genes. B, Relative transcript level of the down-regulated indicated genes. Ten-day-old lightgrown seedlings were treated with 5 mM BA or a DMSO vehicle control for 120 min and were determined by quantitative RT-PCR. The level of transcript relative to the DMSO control is shown as a mean of three biological replicates 6 SE of the mean.
number of genes were linked to biotic and abiotic stress, consistent with a role of cytokinin in the response to pathogens (Choi et al., 2010;Argueso et al., 2012) and abiotic stress (Argueso et al., 2009). An additional bin identified from both the golden list and RNA-seq is comprised of genes involved in secondary metabolism (Fig. 7), and in most cases, these genes are up-regulated in response to cytokinin. This includes genes involved in the biosynthesis of terpenes, flavonoids, and phenylpropanoids. This suggests that cytokinin may upregulate the synthesis of these compounds. Consistent with this, cytokinin increases the synthesis of anthocyanin, a flavonoid compound synthesized via the phenylpropanoid pathway (Deikman and Hammer, 1995). Another striking category in this analysis was a group of genes regulating the redox state of the cell; more specifically, cytokinin caused the up-regulation of multiple genes encoding glutaredoxins. Glutaredoxins are oxidoreductases that use glutathione to catalyze the reduction of disulphide bridges and protein-GSH adducts (S-glutathionylated proteins). Nine of the 30 glutaredoxin genes in the Arabidopsis genome were identified as being up-regulated by cytokinin. The functions of glutaredoxins in plants are only beginning to be elucidated (Rouhier et al., 2008), but include roles in the protection against photooxidative stress (Laporte et al., 2012), ironsulfur cluster assembly, flower development (Wang et al., 2009), and pathogen response (Wang et al., 2009). The cytokinin regulation of multiple glutaredoxin genes suggests a link to these processes.

New Gene Models and Alternative Splicing in Genes Involved in Cytokinin Signaling and Metabolism
In addition to providing information regarding overall transcript levels, RNA-seq read alignments can reveal alternative splicing, alternative transcription start and termination sites, and the relative prevalence of transcript variants. According to The Arabidopsis Information Resource 10 (TAIR10) gene annotations, nearly 20% of annotated genes produce multiple transcripts via one or more of these mechanisms. For most of these, little is known regarding relative expression levels of individual transcript isoforms or which isoforms are the most abundant. Moreover, it is possible that many of these genes may produce novel, previously undiscovered variants. We investigated these questions using Integrated Genome Browser (IGB), a visualization tool that allows in-depth inspection of individual loci (Nicol et al., 2009).
Using IGB, we inspected read alignments and junction files for 62 genes, focusing on gene families whose members are known to play a role in cytokinin signaling or metabolism (Supplemental Table S6). These included cytokinin receptors (AHKs), phosphotransfer proteins (AHPs), type-A and type-B ARRs, CRFs, cytokinin oxidases (CKXs), isopentyltransferases (IPTs), and the Lonely Guy (LOG) gene family.
Of the genes examined, 14 (23%) were already annotated in the TAIR10 genome annotations as producing more than one transcript variant. Of the 14 genes annotated as producing multiple variants, 11 had sufficient coverage to allow identification of the most abundant isoform (Table V). In the majority of cases, more than 90% of spliced reads overlapping differentially transcribed regions supported expression of one isoform, recapitulating a trend observed previously in ESTs (English et al., 2010). We also asked whether the cytokinin treatment triggered changes in the relative amounts of alternative isoforms; we found no significant difference in isoform ratios in the treated versus control samples in any of the alternatively spliced genes, suggesting that alternatively spliced isoforms are produced at constitutive levels independent of cytokinin signaling.
The RNA-seq data also revealed previously unannotated introns and allowed improvement of gene structure annotations. Among the 62 genes examined, 45 new splicing events affecting 28 genes were detected (Supplemental Table S7). However, among these 45 new splicing events, 25 were supported by only one spliced read, suggesting these were either low abundance events or artifacts of library construction or alignment. Of the novel events supported by more than one read, 13 events (in 11 genes) were supported by at least three reads from at least two libraries, suggesting that these events are less likely to be artifacts of library preparation or alignment (Table VI). We selected eight of these well-supported events in eight genes for validation by PCR (Table VII).
RT-PCR analysis confirmed new splicing events in five of the eight genes (Supplemental Fig. S4). For four genes (ARR4, LOG6, LOG8, and IPT2), the validated events included examples of alternative splicing, and new gene models were created that represent the new alternatively spliced isoforms. In the case of ARR17, the validated event added a new exon and intron to the 59 end of the gene and was compatible with the existing model. Rather than create a new gene model   Transcription factors, miRNA, Type A ARRs, histone modifiers AT2G46310 (CRF5) 21 , AT1G74890 (ARR15) 22 , AT1G59940 (ARR3) 23 , AT1G19050 (ARR7) 23 , AT3G48100 (ARR5) 23 , AT3G57040 (ARR9) 23 , AT5G62920 (ARR6) 23 , AT2G41310 (ARR8) 23 , AT2G28160 (bHLH029) 24 , AT2G22770 (NAI1) 25 , AT4G26150 (CGA1) 26,27 , AT1G10480 (ZFP5) 28,29,30 , AT1G16530 (ASL9) 31 , AT1G04240 (SHY2) 13 , AT1G72360 (ERF73) 32 , AT2G40970 (MYBC1) 33 , AT1G71030 (MYBL2) 34 , AT5G65210 (TGA1) 35  for ARR17, we updated the boundaries of its existing gene model AT3G46380.1. In LOG6, the RNA-seq data contradicted the 39 boundary of the final intron. Because the expression level for LOG6 was relatively low, we also examined spliced reads from another RNA-seq data set (Gulledge et al., 2012). Altogether, the RNAseq data contradicted the annotated gene model (AT5G03270.1), consistent with a previous report investigating the LOG6 gene structure (Tokunaga et al., 2012). To capture LOG6 alternative splicing, we created five new gene models, including one new model (AT5G03270. 2) encoding what appears to be a full-length protein compared with the other variants (Supplemental Fig. S4). In most cases, the conceptual translations for the RT-PCR-verified, novel alternatively spliced isoforms were shorter than the annotated gene models' translation. Alternatively spliced isoforms for IPT2 (AT2G27760.2), LOG6 (AT5G03270.3, AT5G03270.4, AT5G03270.5, and AT5G03270.6) and LOG8 (AT5G11950.3) introduced frame shifts and a premature termination codon. The alternatively spliced isoform for ARR4 (AT1G10470.2) contained a new 59 intron that removed the annotated start codon, resulting in a conceptual translation that removed 14 amino acids on the N terminus of the gene's conceptual translation. According to the RNAseq data, most novel isoforms were lower in abundance than their annotated counterparts, with the exceptions of IPT2 and LOG6. In IPT2, both variants were equally abundant. For LOG6 (involved in cytokinin activation), the annotated form was not detected, and the alternatively spliced isoforms appeared equally abundant. Taken together, these results highlight how RNA-seq experiments, in addition to providing new information about differential expression, allow revision and improvement of gene models. However, tools and mechanisms to permit sharing of new models are needed to maximize the impact and relevance of RNAseq data sets. Pending development of new data sharing tools, these new gene models will be made available via the publicly accessible IGBQuickLoad site (www. IGBQuickLoad.org).

DISCUSSION
Various microarray studies have shown that many processes influenced by cytokinin are reflected in changes of the abundance of transcripts encoding proteins of diverse function under specific conditions (Rashotte et al., 2003; Kiba et al., 2004;Brenner et al., 2005;Lee et al., 2007;Argyros et al., 2008;Goda et al., 2008). Our meta-analysis shows that there are many genes identified as cytokinin regulated by only one or a few studies (Fig. 2B), which likely reflects a combination of false positives/negatives and/or context-dependent regulation. The meta-analysis involving microarrays from multiple labs analyzing plants grown under differing conditions and treated with cytokinin in different ways defines a set of robustly cytokinin-regulated genes, a so-called "golden list" that serves as a valuable resource for subsequent studies of cytokinin function and as a foundation for future network analyses in Arabidopsis.
To construct the golden list of cytokinin-regulated genes, we considered only those genes whose transcripts were differentially regulated by at least 1.5-fold in at least 40% of the microarray experiments (i.e. six out of 13). This percentage is higher than used in previous microarray meta-analysis , where only 25% of the microarrays were used without any further filters. While the prior approach is useful to identify as many genes as possible (though only the top 25 were reported), we sought to identify genes that can be assigned as cytokinin regulated with high confidence. Our choice of filters used for the meta-analysis in this work is supported by the observation that with a 1.5-fold change (up-regulation), 73% of the identified genes in six out of 13 microarrays are shown to be up-regulated by an independent experiment (RNA-seq; Fig. 8B). Several published studies further support the quality of the golden list. Taniguchi et al. (2007) identified 23 genes that are direct targets of ARR1, and 20 of these are present on the golden list. The top 25 genes identified as up-regulated by cytokinin that are presented in the  meta-analysis are all found on the golden list. The golden list also contains 39 genes previously validated to be cytokinin responsive by methods other than microarrays (Supplemental Table S1), accounting for 71% of the cytokinin-regulated genes described in the literature. The filtering criteria determined to be optimal for the set of microarrays used here may not be optimal for the meta-analysis of other sets of microarrays, but it does provide a starting point and a reasonable rubric for such studies.
The generation of a list of regulated genes from transcriptome analysis is always a tradeoff between the level of stringency and inclusiveness, trying to achieve an optimal balance between the number of false positives and negatives. The existence of multiple experiments querying cytokinin-regulated gene expression provides an opportunity to identify cytokinin-regulated genes with a higher degree of confidence, independent of experimental variation due to differing growth and treatment conditions. Our meta-analysis focuses on genes regulated by cytokinin across a range of growth conditions, treatment regimens, and, in some cases, distinct tissue types. It is clear that this analysis will exclude genes induced in limited conditions, such as only in etiolated seedlings, or only at specific time points or under specific growth conditions. The hierarchical clustering of the microarray data shows distinct patterns of gene expression in etiolated seedling and roots relative to shoots and whole seedlings. A prior transcript profiling study has shown that roots and shoots display largely similar responses to cytokinin , consistent with the notion that many genes will be induced by cytokinin in multiple tissues and supporting the inclusion of microarray data from both root and shoot tissues in our meta-analysis.
There was an overall 54% overlap in the genes identified as cytokinin regulated by the meta-analysis and the RNA-seq. However, when only up-regulated genes were considered, this overlap rose to 73% and dropped to 14% for the down-regulated genes (Supplemental Fig.  S2). This low percentage of overlap with down-regulated genes may reflect the limitations of RNA-seq and/or of microarrays to identify down-regulated transcripts, the difficulty in robustly identifying transcript levels for genes whose expression level is decreased further by application of a stimulus, or an inherent variability in these genes (Fig. 3). In any case, the lower level of reproducibility for down-regulated genes compared with up-regulated genes has also been found when examining Arabidopsis responses to the phytohormone ethylene (Hall et al., 2012), and thus the difference we observe here is not unusual.
In addition to known early-responsive cytokinin signaling and metabolism genes (type-A ARRs, CKXs, etc.), there were many groups of genes that were overrepresented in both the gene list obtained by microarray meta-analysis and by RNA-seq. A particularly interesting group of genes are the transcription factors, as they likely play a role in mediating the diverse responses to cytokinin. For example, prior research has implicated the cytokinin-regulated CRF and CGA1 transcription factor families in regulating the stimulatory effect of cytokinin on chloroplast division (Rashotte et al., 2003;Okazaki et al., 2009;Chiang et al., 2012). Genetic analysis of the novel cytokinin-regulated transcription factors uncovered here may lend similar insight into how cytokinin regulates its diversity of downstream responses. Several other regulated genes can be associated with demonstrated functions of cytokinin, such as cross talk with other hormones, including auxin (e.g. SAUR like, IAA3 and AXR3), or its functions in cell wall (e.g. AT5G64620; CELL WALL/VACUOLAR INHIBITOR OF FRUCTOSIDASE2 [C/VIF2], AT3G15720; and pectinlyase). Another set of interesting genes identified encodes TNF-receptor associated factor (TRAF)-like proteins (AT5G26290, AT3G20370, AT5G26260, AT5G52330, and AT1G582700), identified by both microarray metaanalysis and RNA-seq. The TRAF-like genes are a large gene family in Arabidopsis encoding potential E3 ligase proteins (Huang et al., 2004;Gingerich et al., 2007;McBerry et al., 2012), only a very few of which have been functionally characterized. These genes may play a role in the regulation of protein turnover in response to cytokinin.
RNA-seq allowed the examination of regulation by cytokinin of nonprotein-coding genes. We found both the up-regulation of miRNA163 and the down-regulation of its potential target, S-adenosyl-L-Met-dependent methyltransferase (AT1G15125), in response to cytokinin. This family of enzymes has diverse roles (Loenen, 2006) in biotic and abiotic stress and in secondary metabolism. RNA-seq also allows the refinement of existing gene models, the definition of novel gene models, and quantification of alternate splicing variants. We used the data from the RNA-seq, coupled with the IGB tool, to analyze the gene models for various genes involved in cytokinin function. This analysis has substantially refined these gene models, which is an important foundation for subsequent analysis of these genes.
The defining of a list of genes that are robustly regulated by cytokinin will help elucidate the mechanisms by which this phytohormone achieves its pleiotropic effects. It provides markers useful for subsequent research and will enable the delineation of a transcriptional network linking cytokinin to other signaling networks in the coordination of growth and development, as well as to biotic and abiotic stress responses.

Plant Materials, Growth Conditions, and Cytokinin Treatment
Arabidopsis (Arabidopsis thaliana; accession Columbia-0 [Col-0]) was used in all experiments. Seeds were sown onto six separate plates containing Murashige and Skoog (MS) media (Murashige and Skoog, 1962) and grown under constant light (150 mmol m -2 s -1 ) at 22°C. After 10 d of growth, seedlings from each plate were collected and added to flasks containing liquid MS media supplemented with 5 mM BA dissolved in DMSO or an equal aliquot of DMSO only. Seedlings from each plate were added to separate flasks. The six flasks were then incubated on a rotating shaker in constant light for 120 min. Samples were collected by and flash frozen in liquid nitrogen. Seedlings from the same plate and the same flask served as biological replicates and were pooled for RNA isolation, complementary DNA (cDNA) library synthesis, and sequencing.

Microarray Meta-Analysis
Microarrays used in this study are described in Table I. Raw data files (.CEL) files were downloaded from NASCArrays (http://affymetrix.arabidopsis.info/narrays/ experimentbrowse.pl) or ArrayExpress (http://www.ebi.ac.uk/arrayexpress/). The raw files were investigated using GeneSpring version 12.1-GX (Agilent Technologies) with default values. Because some microarrays were represented with only a single replicate, no P-value filter was applied. The fold change and P value for each gene were further processed using Microsoft Excel. The P values obtained from this analysis were analyzed using Fisher's method. This test combines P value probabilities from each test into one test statistic (X 2 ) using the formula where p i is the P value for the i th hypothesis test. When all the null hypotheses are true and p i (or its corresponding test statistic) is independent, X 2 has a x 2 distribution with 2k degrees of freedom, where k is the number of tests being combined. This fact can be used to determine the P value for X 2 . In this experiment, k = 13 (13 microarrays used). For the data sets for which no P value was available (one replicate, experiments number 1 and 2), the P value was assumed to be 1 for all the genes. This test provides a meta-P value; a cutoff of meta-P # 0.01 was used to filter the differential regulated genes. The gene list was further filtered using a fold change cutoff of 1.5 or greater and further filtered for genes displaying a fold change of 1.5 or greater in six or more microarray experiments.

Hierarchical Clustering of the Microarray Data
GeneSpring was used to gene-wise normalize the expression data, and hierarchical clustering of signal intensities with respect to gene expression patterns was carried out with Euclidean distance correlation for the two series of experiments, all genes (referred to in GeneSpring as entities) and variables (referred to as conditions). The graphic result was downloaded from the analysis.

In Silico Promoter Analysis
The genomic sequences 1,000 bp upstream of the translational start sites of the genes identified as cytokinin regulated from the meta-analysis (Tables II  and III) were downloaded from the TAIR Web site (http://www.arabidopsis. org/tools/bulk/sequences/index.jsp). This set of sequences was tested for conserved promoter elements, including the core type-B binding element (AGAT [T/C]; Sakai et al., 2000;Hosoda et al., 2002;Imamura et al., 2003), the extended type-B element (AAGAT[T/C]TT; Taniguchi et al., 1998), a cucumber (Cucumis sativus) putative cytokinin response element (TATTAG; Kuroda et al., 2000;Fusada et al., 2005), a TATA box (TTATTT), and Aux-RE (TGTCTC) as a negative control using the Motif Mapper of the TOUCAN 2 workbench (Aerts et al., 2005). Random lists derived from the same number of promoters were generated from genes that were unchanged by cytokinin treatment (i.e. no differential regulation in any microarray). The motifs were counted, and the frequency of their occurrence in the promoters of cytokininregulated genes was divided by their frequency in promoters of an equal number of nonregulated genes.

RNA Isolation
Total RNA was isolated using RNAeasy Mini RNA kit (Qiagen) according to the manufacturer's instructions. Concentration, integrity, and extent of contamination by ribosomal RNA were monitored using the ND-1000 spectrophotometer (Thermo Fisher Scientific) and the Bioanalyzer 2100 (Agilent Technologies).

cDNA Libraries
Barcoded cDNA libraries were prepared for multiplex sequencing on the Illumina HiSeq2000 platform using the TruSeq sample preparation kit according to manufacturer's instructions. Six libraries were prepared from each of three controls and three treatment groups. Starting amount for each sample was 5 mg of RNA. Libraries were named BA1, BA2, and BA3 and Control 1, Control 2, and Control 3 according to the treatment. Libraries BA1, BA2, and BA3 received barcoded TruSeq adapters 9, 10, and 11 and libraries Control 1, Control 2, and Control 3 received adapter barcodes 1, 3, and 8. To obtain sufficient material for sequencing, the libraries were amplified by PCR for 15 cycles following the recommendations of the TruSeq sample preparation protocol. Final elution of each library was in 25 mL of total volume. Library concentrations and quality were assessed using the Qubit 2.0 Fluorometer (Life Technologies) and the Bioanalyzer 2100. BA1 and Control 1, BA2 and Control 2, and BA3 and Control 3 were combined for clustering and sequencing in three separate flow cell lanes. Sequencing was performed on the HiSeq2000 instrument, 50 cycles per flow cell, at the Lineberger Cancer Research High-Throughput Sequencing Center at University of North Carolina, Chapel Hill. Initial sequencing yields for libraries BA1 and Control 1 were low, so a second round of sequencing was performed on these two libraries, thus providing technical replicates for BA1 and Control 1. To distinguish technical replicates, we designate data from these different runs as BA1_1, BA1_2, Control 1_1, and Control 1_2.

Sequence Processing
Sequence data in FastQ format were downloaded from the sequencing facility Web site and processed as follows using a combination of publicly available tools and custom programs. Quality was assessed using FastQC; all samples, including the first low-yielding runs of Control 1 and BA1 libraries, had high overall quality. Sequences were aligned onto the latest Arabidopsis Col-0 genome assembly (released June 2009) using TopHat version 2.0.5 (Langmead et al., 2009) and BowTie version 2.0.5 (Langmead et al., 2009) with maximum and minimum intron size parameters set to 2,000 and 30 bases, reflecting the unusual compactness of Arabidopsis genes compared with mammalian genomes.
The accepted_hits.bam file produced by TopHat was used in subsequent analysis steps. This binary alignment file (Li et al., 2009) contained both spliced and unspliced read alignments. SAMtools version 0.1.18 (Li et al., 2009) and the "NH" (number of hits) Sequence Alignment Map flag were used to separate alignments according to the number of times a read mapped onto the reference genome. Reads mapping to a unique location in the genome were used for differential expression and splicing analysis. The wiggles program (from the TopHat distribution) was used to created coverage graphs in bedGraph format for each binary alignment file. Browser extensible format and bedGraph files were sorted and indexed using the tabix utility (Li, 2011) and deployed on an IGB QuickLoad site to facilitate visualization in IGB (Nicol et al., 2009). The number of singlemapping reads that overlap each annotated gene in this Arabidopsis TAIR10 gene model annotations release was counted and supplied as inputs to DESeq for statistical analysis of differential gene expression. Differential expression analysis of RNA-seq read alignment counts was performed using Bioconductor 2.11 (Gentleman et al., 2004) and Bioconductor library DESeq version 1.2.1 (Anders and Huber, 2010). Probe set annotations used to identify genes not represented on the ATH1 microarray (Gene Expression Omnibus platform accession GPL138) were from the TAIR probe set annotation file available from ftp://ftp.arabidopsis.org/ Microarrays/Affymetrix/affy_ATH1_array_elements-2010-12-20.txt. The false discovery rate was controlled at 5% using the method of Benjamini and Hochberg (1995). For each gene, probe sets that uniquely match each gene (unique.ps) and probe sets that also cross hybridized with related genes (promiscuous.ps) were identified using annotations from the TAIR file.
For splicing analysis, a custom program (FindJunctions) was developed to predict exon-exon junctions from spliced alignments produced by TopHat. FindJunctions was run with options -u (unique) and -n 5 (five flanking bases) so that junction features were predicted from uniquely mapping reads with gapped alignments in which at least five bases of the read aligned to both sides of the gap representing a possible intron. Each junction feature thus represents a distinct splicing event and defines the boundaries of a putative intron inferred from the data. Strand of origin for junctions was inferred from the corresponding genomic sequence using canonical splice site (GT-intron-AG) consensus sequences, where possible. The program (developed in Java) is freely available as part of the open-source GenoViz project (http://sourceforge. net/projects/genoviz/). Novel splicing variants for individual genes were detected by visual analysis of junction features displayed in IGB alongside TAIR10 structural annotations. Numbers of spliced reads supporting annotated or novel introns were tabulated and compared to identify the most abundantly expressed isoform. For genes in which alternative splicing was detected, single-mapping read alignments were used to extend, where possible, the 59 and 39 ends of genes. IGB was used to identify the longest possible reading frame in cases where new splicing events disrupted the annotated open reading frame.

RNA Extraction Quantitative RT-PCR and NanoString nCounter Gene Expression Analysis
Seedlings were grown on MS media (Murashige and Skoog, 1962) for 10 d in continuous light, treated with either DMSO or BA for 120 min, and subjected to total RNA extraction using an RNeasy Plus kit (Qiagen). cDNA was prepared from the total RNA using SuperScript III reverse transcriptase (Invitrogen) as described by the manufacturer. Quantitative RT-PCR was performed using SYBR Premix Ex Taq polymerase (TaKaRa) in a DNA Engine OPTICON 2 (MJ Research). Primers used are listed in Supplemental Tables S8  and S9. Three biological samples each were analyzed with three technical replicates per sample. The relative expression for candidate genes (normalized to b-tubulin as a reference gene and to the wild type grown on DMSO as a control sample) and SEs were determined using REST 2009 software (Qiagen). Error bars represent SE of the mean of the three biological replicates, and Student's t test was performed as a test for significance.

Validation of Splicing Variants by RT-PCR
Wild-type Arabidopsis (accession Col-0) seedlings were sown onto plates containing MS media (Murashige and Skoog, 1962) under constant light of 150 mmol m -2 s -1 at 22°C. After 10 d of growth, seedlings from each of five plates were divided between two flasks containing liquid MS and either 20 mM BA dissolved via DMSO or the comparable quantity of DMSO. Seedlings were treated for 120 min, with 90 min of the treatment on a shaker set to 150 rpm. Seedlings from the same plate and treatment were pooled for each sample and collected with RNALater (Invitrogen). RNA was extracted using a Qiagen RNA mini kit following manufacturer's instructions. cDNA for PCR was synthesized by SMARTMMLV reverse transcriptase using 1 mg of RNA per 20 mL reaction.
Splicing events were tested based using the primers shown in Supplemental Table S10 using a DNAEngine Peltier Thermal Cycler set to either 30 or 35 cycles of 94°C denaturation, 58°C or 60°C annealing, and 72°C elongation using 0.2 mL or more of cDNA in a 25-mL reaction. For most primers, a series of PCR reactions over a gradient of 50°C to 70°C annealing was done to optimize annealing temperatures and check for nonspecific products.

GO and Pathways Analysis
GO term enrichment was performed using the GO Web site (http://amigo. geneontology.org/cgi-bin/amigo/term_enrichment; Carbon et al., 2009) and the DAVID tool (http://david.abcc.ncifcrf.gov/home.jsp; Huang da et al., 2009 ). The graphic result was downloaded from the Web site.

MapMan Analysis
MapMan (Thimm et al., 2004) was downloaded following author's instructions (http://mapman.gabipd.org/web/guest/mapman), and maps were produced to visualize processes affected by the cytokinin treatment. Averaged fold changes of 13 experiments used for microarray meta-analysis was used for the genes from the golden list, and log 2 -transformed fold change for RNA-seq (Supplemental Table S4) served as input to MapMan. Pathway images were saved as .JPG files. Genes were assigned to individual reactions using BINCodes (Thimm et al., 2004).
Sequence data are available from GenBank under BioProject accession no. PRJNA194426.

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S1. GO term enrichment of cytokinin-regulated genes.
Supplemental Figure S2. Comparison of genes from RNA-seq and the golden list.
Supplemental Figure S4. New gene models.
Supplemental Table S2. Genes identified as cytokinin responsive in the microarray meta-analysis by Fisher's combined probability test.
Supplemental Table S3. Golden list genes and functions.
Supplemental Table S5. Genes in functional groups significantly enriched from GO term analysis (DAVID).
Supplemental Table S7. Novel splicing events and supporting reads.
Supplemental Table S8. Primers used for RT-PCR validation of cytokininregulated genes identified in the golden list.
Supplemental Table S9. Primers used for RNA-seq candidates validation.
Supplemental Table S10. Primers used to test novel splicing events.