Research ResourcePlant biology

Combinatorial interaction network of transcriptomic and phenotypic responses to nitrogen and hormones in the Arabidopsis thaliana root

See allHide authors and affiliations

Science Signaling  25 Oct 2016:
Vol. 9, Issue 451, pp. rs13
DOI: 10.1126/scisignal.aaf2768

A nitrogen-hormone interaction network in plants

Both intrinsic factors, such as hormones, and extrinsic factors, such as water and nutrient availability, shape plant development. Ristova et al. evaluated the short-term transcriptional responses and the long-term developmental responses of Arabidopsis thaliana roots to the nitrogen-containing nutrients nitrate and ammonium and the hormones auxin, cytokinin, and abscisic acid. The authors identified genes that were stimulated or inhibited in response to each single factor and to all possible combinations of these three hormones and two nitrogen sources. By combining these transcriptomic data with quantification of changes in root architecture after each treatment, the authors built a multivariate network model of the interaction between nitrogen and hormones that may predict changes in the architecture response of the Arabidopsis root. These data will be useful for future studies on the molecular mechanisms that mediate these interactions and may help identify combinations of nutrients and hormones that improve plant growth under specific environmental conditions.


Plants form the basis of the food webs that sustain animal life. Exogenous factors, such as nutrients and sunlight, and endogenous factors, such as hormones, cooperate to control both the growth and the development of plants. We assessed how Arabidopsis thaliana integrated nutrient and hormone signaling pathways to control root growth and development by investigating the effects of combinatorial treatment with the nutrients nitrate and ammonium; the hormones auxin, cytokinin, and abscisic acid; and all binary combinations of these factors. We monitored and integrated short-term genome-wide changes in gene expression over hours and long-term effects on root development and architecture over several days. Our analysis revealed trends in nutrient and hormonal signal crosstalk and feedback, including responses that exhibited logic gate behavior, which means that they were triggered only when specific combinations of signals were present. From the data, we developed a multivariate network model comprising the signaling molecules, the early gene expression modulation, and the subsequent changes in root phenotypes. This multivariate network model pinpoints several genes that play key roles in the control of root development and may help understand how eukaryotes manage multifactorial signaling inputs.


Living organisms are under the influence of fluctuating and combined external constraints. Sessile organisms, such as plants, cannot flee from adverse conditions and therefore have evolved regulatory mechanisms that help them to adapt to such fluctuations during the course of development. The interplay between external factors and internal signaling molecules is likely at the core of plant adaptations, but the effects of combinations of such molecules remain largely unknown.

Upon changes in nutrient availability, primary regulatory functions very often involve a rapid fine-tuning of growth processes through the modification of the hormonal status of the plant [reviewed in (13)]. Furthermore, these growth processes can then feed back through hormones to regulate the expression of metabolic genes, which are often at the forefront of nutrient sensing. This cycle is likely to be supported by multilayered dedicated signaling pathways (2). At a genetic level, genes involved in hormone sensing and signaling control some aspects of plant adaptation to fluctuating nutrient conditions. For example, this is true for the auxin-resistance gene AXR4 (4), the abscisic acid (ABA)–related transcription factors ABI4 and ABI5 (5), the auxin-related transcription factor ARF8 (6), the auxin receptor AFB3 (7), the ethylene signaling components ETR1 and EIN2 (8), the cytokinin (CK) signaling–related proteins ARR and AHK (9), and the proteins involved in CK biosynthesis: IPT3, IPT5, and IPT7 (10). Conversely, nutrient sensing mutants are affected in their response to hormones. For example, this is true for the hexokinase (sugar sensor) HXK1 (11) and the nitrate transceptor (nitrate sensor) NRT1.1 (12). This line of evidence suggests that there are dedicated signaling pathways mediating a feed-forward cycle connecting nutrient and hormone signaling (2, 3).

Genome-wide approaches (13) have refined the definition of plant response to external factors as well as to hormonal treatments (13, 14). Hormonal treatments are usually made under standard nutrient conditions, and very few studies have reported the effects of combinatorial treatments (15, 16). Furthermore, the effects of combinations of hormonal treatments together with nitrogen availability on genome-wide expression reprogramming remain largely unknown. Indeed, hormonal treatments tend to be usually performed in full nutrient media, and it is rare to perform hormonal treatments on nutrient-starved versus well-fed plants. Thus, how nutrient availability affects the responses to hormones and how the presence of hormones affects the responses to fluctuating nutrient conditions are largely unexplored.

The plant root is a complex organ, and its development is very sensitive to hormone and nutrient availability (4, 1723). The root system is also an important target for selective breeding and biotechnological approaches (24). Many questions about how combinations of nutrients and hormones affect root architecture remain. For example, does the auxin-mediated arrest of primary root (PR) growth that is followed by lateral root (LR) outgrowth (17, 25) occur under all nutrient conditions? Do CKs affect responses to auxin or ABA? If so, how do the roots respond when treated with a combination of ABA, CK, and auxin simultaneously? These are just a few questions that can only be answered by studying complete combinatorial treatment matrices of nutrients and hormones (2628).

Thus, here, we studied the effect of all combinations of five nutrient and hormone signaling molecules [nitrate (NO3), ammonium (NH4+), auxin, CK, and ABA] known to affect root development and influence one another. We recorded short-term (hours) genome-wide variations in the transcriptomic profiles of treated roots and examined the effect of these combinatorial treatments on root development following days of treatment. Finally, we used modeling methods to interconnect both data sets in a comprehensive multivariate model that embodies signaling molecules, gene regulation, and root parameters. This approach provides an important resource data set and pinpoints the role of several loci in controlling root development and potentially mediating signal integration.


Combinatorial treatment for discovery of signal interactions

To study the effects and interactions of five nitrogen nutrient and hormone signaling molecules, we developed a matrix of all 32 possible combinatorial treatments (table S1), including NO3 (1 or 0.5 mM), NH4+ (1 or 0.5 mM), the auxin indole-3–acetic acid (IAA; 500 nM), the CK kinetin (CK; 500 nM), and ABA (1 μM). When NO3 and NH4+ were used in combination, the total nitrogen treatment was maintained at 1 mM to record the effect of these specific nutritional molecules and not the effect of total nitrogen content. Mock treatments included KCl and dimethyl sulfoxide as controls for salt and solvent provision, respectively (table S1). The signal molecule concentrations have been chosen to be comparable, yet lower than those used in a previous study (14) to record clear marker gene activation and thus be able to record signal interaction. Arabidopsis thaliana plants were grown on ammonium succinate for 12 to 14 days until the first true leaves formed, transferred to fresh N-free media for 24 hours, and then transferred to fresh media containing the combinations of signals being tested. Roots were harvested after 4 hours of treatment for gene expression analysis or left for 4 days on agar plates containing the signaling molecules (table S1) for imaging and subsequent analysis of root architecture.

Root response to combinatorial treatments

The roots of wild-type Arabidopsis Col-0 seedlings subjected to the 32 treatment conditions (12 to 20 plants from two to three independent experiments) were scored for developmental characteristics. We measured 10 root traits (Fig. 1A). Results are presented as box plots (fig. S1) and as a cluster of mean values (Fig. 1B). We performed analysis of variance (ANOVA) on the data to identify the significant effects of treatment and signal interaction in the control of each trait (Fig. 1C). As expected, auxin had the most marked effect on root development, reducing PR growth and increasing LR length and density (1719). Auxin interaction with CK and ABA was very pronounced when molecules were provided in combinations (Fig. 1C and fig. S1). For instance, CK repressed PR length but only in the absence of auxin. Inversely, CK counteracted the inductive effect of auxin on LR number and density. ABA repressed PR and LR length (fig. S1). It is noteworthy that NO3 and NH4+ provision influenced root development. NO3 effect was the most pronounced on LR density and length. NH4+ affected LR density through LR number (29). Interestingly, despite the fact that NO3 and NH4+ were known to affect nitrogen assimilation pathway by interacting with one another (3032), here, we detected no significant interaction concerning their influence on root development (Fig. 1C; see cluster line named NH4*NO3). By contrast, NO3 and NH4+ each clearly interacted with hormonal signaling. A significant NO3 and IAA interaction was recorded for the control of PR and LR length, but the most marked interaction was recorded for the inductive effect of NH4+ on IAA control of LR density. Indeed, auxin promoted LR density (from 2 to 6 LRs per centimeter of PR), but this induction was enhanced when NH4+ was present in the media (from 2 to ~8 LRs per centimeter of PR) (Fig. 1, B and C, and fig. S1). We also detected more complex phenomena. For example, we observed a very strong signal interaction between CK, ABA, NO3, and NH4+ in the control of the PR response. This was illustrated when the neo-grown PR (P2, root grown after transfer) was measured. For instance, ABA repressed root growth in the absence of nitrogen and CK (NO3- or NH4+-free conditions; fig. S1, compare bars 1 and 2). This repression was lost in the presence of nitrogen (fig. S1, compare bars 9 and 10, bars 17 and 18, and bars 25 and 26). This indicates that nitrogen (NO3 or NH4+) modified the effect of ABA on PR growth. Inversely, ABA repressed PR growth in the presence of NO3 and CK (fig. S1, compare bars 19 and 20, bars 27 and 28) but not in the absence of NO3 or in the presence of NH4+ alone (fig. S1, compare bars 3 and 4, bars 11 and 12).

Fig. 1 Hormones and nitrogen interact to shape root system architecture.

(A) Scheme of Arabidopsis root traits that were scored. P, total PR length; P1, PR length before treatment; P2, PR growth after treatment; LR.length, mean length of all LRs. (B) Clustering of root traits. Root trait values have been mean-centered normalized to display their responsiveness to the treatment matrix. Absolute quantifications are available in fig. S1. LR.nb, number of visible (>0.5 mm) LRs; LR.nb.P1, number of visible (>0.5 mm) LRs in P1; LR.dens, number of LRs per centimeter of P; LR.P.ratio, ratio of LR length to PR length; T.LR.L, total LR length; T.R.L, total root length (primary plus lateral roots). (C) ANOVA results. Heat map display of F values, which represent the strength of regulation by a given factor or combination of factors. The dendogram values reflect the distance used here (Pearson correlation). Gray color indicates no significant effect (P > 0.05) of the factor on a particular root trait.

These combinatorial effects suggest the presence of complex signal integration mechanisms that enable roots to integrate input from multiple signals and adapt their development accordingly. These root architecture adjustments show that NO3 or NH4+ can influence hormone action at multiple levels. Much more can be gleaned from this data set, but for the purposes of clarity and concision, we present only a subset of signal interactions to illustrate the prominent and fine regulations (Fig. 1C and fig. S1).

Transcriptional responses to combinatorial treatments

Genome-wide investigation of gene expression in response to the 32 combinatorial treatments of nitrogen and hormones was recorded using an A. thaliana genome array (Affymetrix ATH1) (32 conditions × 2 independent experiments). These data sets were modeled through ANOVA as previously described (26, 33, 34) (results are provided in table S2). We considered a gene to be regulated by any factor or combination of factors if P < 0.01 for a given factor or a combination of factors (for details, see Materials and Methods). This modeling retrieved 9507 regulated probe sets (referred to hereafter as “probes” for simplification). Clustering analysis was first performed on NO3 (fig. S9), NH4+ (fig. S10), IAA (fig. S11), CK (fig. S12), and ABA (fig. S13) responsive genes separately. To clarify the most predominant modes of genome regulation and signal management, (i) we grouped and sorted gene expression models according to the number of genes that they contain (Fig. 2A and fig. S2) (26), and (ii) we clustered the ANOVA results to monitor the signal convergence at the genome-wide level (Fig. 2B; see Materials and Methods). In brief, genes were used as sentinels of signal convergence. For example, if a set of genes was controlled by ABA and IAA, it marks a significant signal convergence. At the genome-wide level, this represents a survey of signal convergence in the control of gene expression (Fig. 2A and fig. S2) (26).

Fig. 2 Genome-wide features of hormone-hormone and nitrogen-hormone interactions.

(A) Dominant models of gene regulation classified by the number of probes affected by a given factor or combination of factors. (B) Dendogram illustrating the convergence of the different signaling molecules and their combinations on the control of the gene expression. (C) Heat map display of the degree of overlap (z score) between lists of gene controlled by simple signals and binary combinations of signals. (D) Examples of an AND logic gate behavior in the control of mRNA accumulation. The expression of At2g32280 increased only when CK and IAA were both present in the media, and the expression of At4g18560 increased only when ABA and NO3 were both present in the media. The histogram (means ± SEM of each bar, n = 16) was built from the 64 data points (25 = 32 treatments × 2 biological replicates).

The most influential factors were ABA, which affected 3499 probes, and IAA, which affected 365 probes (Fig. 2A, models 1 and 2, respectively). The next most influential signal was NO3*CK, which controlled 328 probes (model 3). This indicates that a large set of genes was regulated differently when NO3 and CK were present in combination and not regulated by either as a simple (single) factor alone. To understand this important phenomenon detected by ANOVA, we clustered genes according to their expression (fig. S3). CK tended to invert the gene expression in response to NO3. A gene (and relatives in the cluster) that tented to be repressed by provision of NO3 (At1g15480; fig. S3) was induced by NO3 when CK was present, and the other way around (see At5g14320; fig. S3). Although the magnitude of the interaction was not very marked, this interaction was significant (passed ANOVA threshold), and the importance was revealed by the number of genes responding to this mode of regulation (figs. S2 and S3). This analysis identified a new level of interaction between NO3 and CK signals in addition to those previously described (2, 9, 10).

The next model of regulation, including a large number of genes, corresponded to genes regulated by ABA and IAA or ABA and CK without interaction (Fig. 2A, models 4 and 5). This meant that a large number of genes were controlled at the same time by these pairs of signals and that their effect was largely additive. It is noteworthy that, again, the effect of CK largely overlapped with the effect of ABA and that this combination controlled more genes than CK alone. These observations, together with previous statements, show that CK signaling is highly dependent on other signals. Finally, the next three simple models (models 6, 7, and 8) corresponded to NH4+ alone, NO3 alone, and CK alone, which controlled 198, 183, and 175 probes, respectively. This logic can be followed for the top 100 models (fig. S2).

To probe genome-wide signal convergence, we built a signal cluster using ANOVA-binarized matrices as a template (Fig. 2B) (26). If two signals clustered together, this indicated that they controlled similar sets of genes. At a genome-wide level, (i) NO3, NH4+, and CK simple signals and (ii) NO3*CK and NO3*NH4+ composite signals controlled a large set of common genes. This NO3-NH4-CK cluster was related to the IAA signal likely through CK convergence. To support these observations, we performed a randomization test using the GeneSect algorithm (35) (available at that monitors the significance of gene list overlap. We observed that the overlap between gene lists was more significant than expected for most of the signals we studied (P < 0.001). However, z scores themselves can highlight the stronger signal convergence (Fig. 2C). For example, the set of genes controlled by NO3 strongly overlapped with genes controlled by CK, ABA*NO3, and NH4+. A strong convergence was also noted between CK and IAA and between IAA and ABA (Fig. 2C).

We then explored whether the biological functions related to the gene lists displaying simple signal control (such as ABA or IAA) and composite signal control (such as NO3*CK) were specific to a particular signal combination. To do so, we subjected the different gene lists corresponding to the dominant mode of gene regulation (defined in Fig. 2A) to GeneCloud analysis, which reports semantic terms that are highly represented in the particular gene list (36). Particular functions seemed to be specifically targeted by different combinations of signals (table S3). For instance, ABA as a simple signal (model 1) controlled genes related to the following terms: “chloroplast” [GeneCloud Fisher test, false discovery rate (FDR)–corrected P = 9.94 × 10−19], “aba” (P = 4.16 × 10−5), or “abscisic” (P = 3.03 × 10−4). However, genes targeted by ABA*IAA (model 4) did not display significant enrichment in these terms but rather displayed strong enrichment for the terms “hair” (P = 1.63 × 10−6), “extensin-like” (P = 3.18 × 10−6), and “shade” (P = 2.53 × 10−4). Similarly, ABA*CK-targeted genes (model 5) were enriched in the terms “cysteine-histidine-rich” (P = 1.71 × 10−4), “vacuole” (P = 6.34 × 10−4), and “iron” (P = 4.44 × 10−3). These observations show that simple signals or combinations of signals control genes likely related to specific functions. Likewise, NO3*CK-controlled genes (model 3) were enriched in the terms “kip1-like” (P = 4.47 × 10−3), “nadh-ubiquinone-plastoquinone” (P = 6.14 × 10−3), “deaminase” (P = 1.43 × 10−2), and “spermine” (P = 4.01 × 10−2). On the other hand, genes regulated by NO3 alone (model 7) were enriched in the terms “trehalose” (P = 2.16 × 10−3), “nitrogen” (P = 8.92 × 10−3), and “plastidic” (P = 2.81 × 10−2), and genes regulated by CK alone (model 8) were enriched in the terms “arr” (P = 4.52 × 10−3), “calcium-dependent” (P = 5.76 × 10−3), and “differentiation” (P = 1.86 × 10−2). The fact that these terms were not enriched in the genes controlled by NO3*CK showed again that signal combination was likely targeting specific regulatory modules. We detected a strong signal convergence between IAA and CK, between IAA and ABA, between NO3 and NH4+, and between CK and ABA (Fig. 2C). In evaluating the terms enriched for these pairs of signals (table S3), we noted that NO3 and NH4+ converged to control genes related to the redox state of the cell or auxin-related genes, that IAA and CK converged on a high number (11) of peroxidases, and that IAA and ABA converged to control root hair–related genes (among other terms, see table S3).

A very strong signal convergence has been observed above through genome-wide analysis but can be manifested through single gene analysis. To also exemplify the power of combinatorial approaches at the single gene level, we retrieved genes whose expression patterns highlight an AND logic gate behavior in signal interactions (Fig. 2D). Finally, a set of genes previously known to be responsive to the studied signals (14, 32, 37, 38) and their interactions were validated by quantitative real-time polymerase chain reaction (qRT-PCR) on the two biological replicates used for ATH1 chips plus one biological replicate that was not used for genome-wide assay (fig. S4).

Identification of individual signal markers

According to the results presented above, signal interaction is a general rule rather than an exception. For instance, of the 1545 genes controlled by the nutrient-related signals NO3, NH4+, or both together (union of gene lists of NO3, NH4, and NO3*NH4+), 982 genes (64%) were also regulated by a hormonal signal (fig. S5).

Thus, conversely, this data set was also a great occasion to isolate gene markers of the simple signals as those particularly responsive to only one particular stimulus. We used pattern matching to identify genes that are the best markers for each of the five simple signals as well as markers of nitrogen addition (NO3-OR-NH4+ logic gate behavior) (fig. S6). Whether we treated plants with NO3 or NH4+ alone or with NO3 and NH4+ simultaneously, we maintained the total nitrogen content of the media at 1 mM (table S1), enabling us to identify genes that were regulated by the presence of nitrogen. The accumulation of the BT5 (At4g37610) transcript was one such marker of the availability of nitrogen. Other marker genes that were regulated in response to a simple signal included the ARR5 (At3g48100), the expression of which was specific to CK treatment, as well as IAA13 (At2g33310) for auxin, NIR1 (At2g15620) for NO3, GDH2 (At5g07440) for NH4+, and ABI1 (At4g26080) for ABA, among others (fig. S6 and table S4).

A multidimensional model linking signals to root developmental parameters through gene expression modulation

To unify our observations of the effects of the five nutrient and hormonal signals and their interactions on gene expression and root development, we developed a model linking these two levels of information (Fig. 3A). This integrated network model included the rapid transcriptional response assessed 4 hours after treatment, the long-term phenotypic adaptation of roots over 4 days after treatment, and the effect of combinatorial treatments. This approach makes the assumption that early gene expression modulation may reflect the root developmental plasticity later observed. This idea is supported by a large body of literature showing that important transcriptional reprogramming events precede adaptive root development (3942) and that key regulatory genes, such as for example ACR4 (39, 43), PUCHI (44), IAA28 and GATA23 (41), and AFB3 and NAC4 (7, 45), are included among these regulated genes.

Fig. 3 Multivariate model linking signaling molecules to root developmental response through early gene expression modulation.

A multidimensional model was built linking root traits to early gene regulation. Nodes are of three types: phenotypic traits (yellow hexagons), genes (gray diamonds), and signals [ABA (green circles), IAA (red circles), CK (dark blue circles), NH4 (purple circles), and NO3 (light blue circles)]. Black squares depict interactions between signals. Edges were built from two kinds of information: edges between phenotypic traits and genes are based on Pearson correlations [cor > 0.7 (red edge); cor < −0.7 (blue edge); randomization test, P < 10−5], with the thickness of the edge proportional to the correlation value. Edges between signals and genes were based on ANOVA results (influence of signals on gene expression, P < 0.01). (A) Hypothetical scheme showing a gene (At) that is transcriptionally controlled by ABA, IAA, and an interaction between CK and IAA. (B) The complete model linking phenotypic traits, 280 genes, and the five signals. The complete network is provided as a Cytoscape readable file (Nitrogen+Hormones.cys). (C) Subset of the multidimensional network presented in (B) extracting the potential influence of the genes already known to control developmental parameters and their connections to the five studied signals. (D) Venn diagram and heat map display of the degree of overlap (z score) between sets of genes related to LR development as identified by key previous studies (39, 42) and our study.

This integrated model (Fig. 3B) reiterated that IAA and ABA were the two major forces shaping root architecture under our experimental conditions. IAA influenced PR and LR growth through its effect on LR number and P2 length, whereas ABA mainly influenced parameters related to LR length (Fig. 3B). This model connected signals to traits through gene expression. A variant of the model directly linking traits to the signals (fig. S7) that revealed gene modules (46) related to a particular group of signals and particular traits.

To evaluate the predictive power of the network model (Fig. 3B), we performed an extensive literature search through the TAIR (The Arabidopsis Information Resource) website ( to see which genes in the network were already experimentally associated with root development. A large number of genes in this model had previously been linked to root development. This included genes encoding the kinase PINOID (47); the transcription factors GATA23 (41), PUCHI (44), IAA18 (48), and LBD29 (4951); and many others (Fig. 3C). We compared the predictions of this model to two landmark transcriptomic data sets reflecting LR developmental programs (39, 42). This comparison of gene lists regulated in response to treatments that initiate LRs, together with the genes belonging to the model presented in Fig. 3D, showed a high degree of overlap. The amount of convergence was evaluated by randomization testing (35). We found that the genes we identified, whose expression may predict root phenotypic traits, belonged to common modules regulated upon LR activation (fig. S8). The overlap between these three gene lists was greater than would be expected by chance (Fig. 3D). Furthermore, we investigated semantic terms overrepresented in this gene list. Genes with terms related to “root” (GeneCloud Fisher test, P = 6.97 × 10−19), “auxin” (P = 4.93 × 10−5), and “inhibitor-lipid-transfer” (P = 2.27 × 10−4), along with others, were present more than expected by chance. The InterPro domain ipr006459 (P = 3.68 × 10−5) was found in five genes (CASPL1C2 At1g03700, CASPL4D1 At2g39530, CASPL2B2 At2g35760, CASPL1D1 At4g15610, and CASPL1E1 At4g15630) from this list, which was 23 times greater than the representation of this domain in the entire genome. This protein domain is found in Casparian strip membrane domain protein (CASP)–like proteins, which have four membrane-spanning domains (52). CASPs recruit the lignin polymerization machinery necessary for formation of the Casparian strip, a specialized region of the cell wall that prevents the passive flow of solutes in the endodermis and thus modifies its cell wall properties (53). Knowing that endodermis cells contiguous with the initiating LRs are adapting themselves to facilitate the formation of the LRs (54), it is interesting to hypothesize that expression of these five genes correlates with subsequent LR development, perhaps through cell wall modification.

The multivariate network model to reveal new genes involved in shaping root architecture

Our network model (Fig. 3B) highlighted genes already known to shape root architecture, but we also sought to determine whether it could reveal new ones. We scored root development in the presence of nitrogen (1 mM NO3) in 10 knockdown or knockout mutant lines (fig. S14) that affect genes distributed throughout the multivariate network (Fig. 4). The respective mutants displayed altered root development phenotypes when compared to wild-type plants (Fig. 4A). The genes affected in these mutants encode diverse proteins and display clear transcriptional regulation in response to gravistimulation, which induces LR formation and emergence (Fig. 4B) (40). Together, these observations demonstrate that our model can predict genes that affect root development. Further studies will of course be required to determine the molecular mechanisms by which these gene products affect root development.

Fig. 4 Identification of potential regulators of LR development.

(A) Heat map of mean-centered normalized values showing the root development phenotypes when the indicated mutant plants were grown on standard media containing 1 mM KNO3 [ratio = mutant/wild type (WT)] (gray, not significant from the WT REML mixed model). LR.No, LR number. (B) Transcriptional responses concomitant to LR development after gravistimulation. Data are means ± SEM, n = 4 (40).

To demonstrate that our approach could also identify genes involved in mediating responses to the factors we tested, we focused our analysis on the gork mutant phenotype. The GORK protein is an outward-rectifying potassium (K+) channel that is highly abundant in guard cells (55). GORK mediates K+ fluxes that modulate cell turgor and thus control stomata closing (56). Although GORK is a well-characterized protein, it was not previously linked to root development. ABA promoted the accumulation of gork transcripts by 10-fold in the roots of seedlings, reaching Affymetrix MAS5 signals above 400, which was clearly above background expression. This regulation of gork by ABA treatment was consistent with previous reports of ABA-stimulated gork expression in other tissues (57, 58). Two independent gork mutant lines (gork-1 and salk_082258) displayed LR elongation and meristem size phenotypes (Fig. 5A and fig. S15). Furthermore, the modulation of LR length by ABA was also affected in the gork mutants in different nitrogen-containing media (Fig. 5, A to C, and fig. S15). This shows that the modulation of gene expression within hours may be related to subsequent phenotypic adaptations. We also noted that the meristem of the PR and LR of the salk_082258 mutant display a cell swelling phenotype (fig. S15). These preliminary insights provide a basis for future investigation of the role of GORK in root development. Altogether, our model (Fig. 3B and data file S1) may be a useful resource for the identification of genes involved in the hormonal and nutritional control of root architecture.

Fig. 5 Mutations in the gene encoding the K+ channel GORK affect LR development and the response of roots to ABA.

(A) ANOVA for six root traits. Plants (WT and gork mutants) were grown on 0.5 mM ammonium succinate media for 10 days and transferred to the same media [control or supplemented with abscisic acid (1 μM ABA)]. After 4 days plates were scanned and root traits quantified. (B) Two independent gork mutant lines (salk_082258 in the Col-0 background and gork-1 in Ws background) display root development phenotypes, both under control conditions (white boxes) and on media containing ABA (gray boxes), n > 10. (C) Representative examples of root phenotypes. Genotypes and the presence or absence of ABA are indicated above the images.


It is technically simpler to study the effects of signaling molecules by applying them one at a time to reveal their role in a particular phenomenon. However, in nature, biological systems are not exposed to isolated stimuli. Instead, they are exposed to combinations of diverse types of signals. It is therefore likely that molecular networks evolved to be responsive to combinations of molecules rather than individual factors (27). Here, our goal was to integrate the genomic and developmental responses to nutritional (nitrogen) and hormonal cues in the model plant Arabidopsis. We believe that this resource provides a first step toward a better understanding of such signal interactions.

Signal interactions in the control of root development

At the developmental level, each of the five individual signals (NO3, NH4+, IAA, CK, and ABA) has a different effect on LR traits, recapitulating previously reported effects of these molecules (17, 59). We found that combinations of these signals strongly influence LR traits as well (Fig. 1 and fig. S1). Whereas IAA and CK interact mainly in the control of PR growth and influence LR density, ABA and IAA have a combined effect only on LR length (Fig. 1C). In our experimental system, ABA*CK did not affect root traits (Fig. 1C). This absence of interaction was not reflected in the transcriptomic response. Indeed, CK and ABA tend to control a largely overlapping set of genes (Fig. 2C) that are potentially involved in other physiological functions. To get better insight into what functions may be coregulated by ABA and CK signals, we performed a GeneCloud analysis (36) on the intersection of the ABA- and CK-responsive genes. This revealed the semantic terms “urea” (Fisher test, P = 2.95 × 10−4), “suberin” (P = 2.95 × 10−4), and “thalianol” (P = 2.91 × 10−5). This demonstrates that, even whether ABA and CK do not seem to interact in the control of root architecture, they may converge to the control of other physiological processes.

Signal interactions in the control of genome-wide transcription

We have provided a first step toward a comprehensive analysis of combinatorial interactions at a genome-wide scale in a multicellular biological system. Our results demonstrate that the signaling molecules we studied (hormones, NO3, and NH4+) largely influence one another’s effects because a large proportion of genes are influenced by a signal interaction (meaning that the gene expression cannot be deduced from the addition of the individual signal effects). For instance, if we consider the whole set of 9507 regulated probes at a threshold of P < 0.01, 4006 probes were also influenced by a composite signal (fig. S16). In other words, ~42% of the response is managed by an intricate web of regulation where each studied signaling molecule influences the signaling pathways of the others. Signal interactions between hormones have already been revealed at several levels of integration (60, 61). Here, we report the genome-wide magnitude of such interconnections between hormones and between hormones and nutrient signals and propose that this is a key feature of the underlying regulatory networks.

Nutrition and hormone interconnections

Because nitrogen-containing compounds are fundamental components of plant anabolism, these nutrients are absolutely essential for sustaining plant growth and development. Thus, a reduction in nitrogen provision may affect plant growth through simply the law of mass action, which is very likely during long periods of starvation. However, experimental investigations also demonstrated that rapidly responsive and dedicated signaling pathways may also exist for detecting available nutrients and coupling this information to the hormonal signaling pathways that direct development (2, 3). For instance, Walch-Liu et al. reported that shoot growth is very quickly modified upon nutrient treatment, even before the nitrogen content of the shoot has been modified (62). The current data set provides many entry points into the potential networks coupling nutritional sensing and hormonal control.

First, at the transcriptomic level, it is interesting to note that the CK signal tends to control many common genes together with the NO3 and NH4+ signals (Fig. 2, B and C). Second, NO3 and CK tend to deeply influence one another’s effect on several hundred genes (Fig. 2A and fig. S3). Third, the hormone-responsive genes are deeply influenced by NO3 and NH4+ signals. For instance, if we only consider the genes controlled by individual signals and the first order of composite signals, such as IAA*CK, 33% of the IAA-regulated genes are also regulated by NO3 or NH4+, and 55% of the CK-regulated genes are also regulated by NO3 or NH4+. Fourth, many genes that correlate with root development parameters in our data set (Fig. 3B) are influenced by NO3 or NH4+ more than expected by chance [randomization test, z score = 12 (for NO3) and 7 (for NH4+)]. Potential influences of the NO3 and NH4+ signals are also manifested in the LR transcriptional response during development (fig. S8B) (40). Fifth, even at the level of single genes, the interaction between signals is obvious. For example, the gene At4g18650 is under the control of an ABA-NO3 AND logic gate pathway (Fig. 3D). We believe that the evolution of such a precise transcriptional mechanism entangling ABA and NO3 signals is an example of the importance of hormonal and nutritional interactions. Last, genes identified by previous studies as reporters of hormone activity, such as NAC4 for ABA or ARR8 for CK (14), were found in our study to be regulated by nitrogen nutritional signals as well, demonstrating once more that hormonal and nutritional signaling pathways are very closely connected (fig. S4).

In conclusion, we provide a developmental and transcriptional landscape of nutritional and hormonal interactions. This study provides insights into a largely unexplored world of nitrogen and hormone signal interactions and is a first step toward elucidating how combined signaling cues interact to shape a sessile organism whose plasticity is crucial for survival and adaptation to fluctuating environmental conditions. These phenomena are critical for plant adaptation and are thus crucial to maintain food networks in the face of a changing climate.


Plant culture and treatment

A. thaliana plants, ecotype Columbia (Col-0), were grown in sterile hydroponics as adapted from (6). Sterilized seeds were sown on Nitex 03-250/47 mesh (Sefar America) supported by a plastic platform to allow roots to grow in hydroponics inside a sterile Phytatray (Sigma-Aldrich). Hydroponic media consisted of 1× Murashige and Skoog basal medium containing no nitrogen or sucrose (custom-ordered from Gibco BRL) supplemented with the following: 3 mM sucrose, 0.5 mM ammonium succinate, and MES buffer (pH 5.7) (0.5 g/liter). Plants were grown for 14 to 16 days in day/night cycles (16/8 hours, 150 μmol photons m−2·s−1) (Percival Scientific Inc.) at 22°C. Twenty-four hours before the treatments, plants were transferred to equivalent fresh N-free media. Treatments (table S1) were then applied for 4 hours. Root development analysis was performed under the same conditions in petri dishes on agar-solidified media. Four days after treatment, plates were scanned and root development parameters were scored with Optimas image analysis software as previously described (10, 63).

Mutant strains

Mutant strains salk_098602, salk_082258, salk_095938, salk_136802, salk_069429, salk_023374, salk_117623, salk_074376, salk_019552, and salk_125384 (in the Col-0 background) were obtained from the Arabidopsis Biological Resource Center stock center. The gork-1 mutant was provided by J.B. Thibaut.

Gene expression analysis

Total RNA extraction was performed using TRIzol reagent according to the manufacturer’s recommendations (Invitrogen). For qRT-PCR, 1 to 2 μg of total RNA was digested by deoxyribonuclease I (Sigma-Aldrich). RNA was then reverse-transcribed to generate single-stranded complementary DNA (cDNA) using ThermoScript Reverse Transcriptase (Invitrogen) according to the manufacturer’s protocol. Gene expression was measured by qRT-PCR (LightCycler; Roche Diagnostics) using gene-specific primers and LightCycler FastStart DNA Master SYBR Green (Roche Diagnostics). The abundance of transcripts was normalized to ACT2, ACT8, and Clathrin as previously described (31). For microarray analysis, cDNA was synthesized from 2 μg of total RNA. Labeled cDNA (8 μg) was hybridized to an Arabidopsis ATH1 gene chip (Affymetrix) for 16 hours at 45°C. Washing, staining, and scanning were performed as recommended by Affymetrix. Image analysis and normalization to a target median intensity of 150 were performed with the Affymetrix MAS5 set at default values. We analyzed the reproducibility of replicates using the correlation coefficient and visual inspection of scatterplots of pairs of replicates. The raw data (.CEL files) have been deposited in the Gene Expression Omnibus (GEO) database (accession no. GSE71737).

Modeling of gene expression patterns

All data manipulations were performed on R ( The ANOVA analysis was carried out using the R aov() function on log 2 MAS5-normalized data (26, 33, 34). A probe signal has been modeled as follows: Yi = α1.ABA + α2.NO3 + α3.CK + α4.IAA + α5.NH4 + α6.ABA*NO3 + α7.ABA*CK + α8.NO3*CK + α9.ABA*IAA + α10.NO3*IAA + α11.CK*IAA + α12.ABA*NH4 + α13.NO3*NH4 + α14.CK*NH4 + α15.IAA*NH4 + α16.ABA*NO3*CK + α17.ABA*NO3*IAA + α18.ABA*CK*IAA + α19.NO3*CK*IAA + α20.ABA*NO3*NH4 + α21.ABA*CK*NH4 + α22.NO3*CK*NH4+ + α23.ABA*IAA*NH4 + α24.NO3*IAA*NH4 + α25.CK*IAA*NH4 + α26.ABA*NO3*CK*IAA + α27.ABA*NO3*CK*NH4 + α28.ABA*NO3*IAA*NH4 + α29.ABA*CK*IAA*NH4 + α30.NO3CK*IAA*NH4 + α31.ABA*NO3*CK*IAA*NH4 + ε, where α1 to α31 represent the coefficient quantifying the effect of each of the factors (ABA, NO3, CK, IAA, and NH4) and their interactions, and ε represents the nonexplained variance. We determined the FDR to be <5%.

Clustering algorithm, Sungear analysis, and interpretations

Hierarchy between signals was evaluated by average linkage hierarchical clustering. First, Euclidean distances were calculated using the dist() function in the R software. Second, clusters were generated by the hclust() function. Last, plots were generated using the plot() function with default values. Dendrogram interpretations were carried out as previously described (26). In brief, if a given gene behaves similarly in response to two factors (for example, IAA and CK), it will decrease the distance (increase the linkage) between those two factors. Generating dendograms allows one to visually capture the relative relationship of the signals in the control of the regulation of the gene set considered (Fig. 2A). Note that branch length is set to a constant value and is not related to the data [plot() function with default values]. Only the height of the node reflects the distance between the branches and the associated leaves of the tree. Heat map clusters were generated with the MeV software using Pearson correlation as distance (


Fig. S1. Hormones and nitrogen interact to shape root system architecture (absolute quantification).

Fig. S2. Revealing genome-wide features of nitrogen-hormone interactions, the top 100 models.

Fig. S3. Study of the NO3*CK-only responsive genes.

Fig. S4. qRT-PCR validation of sentinel genes.

Fig. S5. Sungear figures (generalized Venn diagram) measuring the convergence of nutritional and hormonal signaling pathways.

Fig. S6. Identification of marker genes of simple signals.

Fig. S7. Multivariate model with signal-to-trait connections.

Fig. S8. Clustering analysis of genes correlated with root traits.

Fig. S9. Clustering of all the genes controlled by NO3.

Fig. S10. Clustering of all the genes controlled by NH4+.

Fig. S11. Clustering of all the genes controlled by IAA.

Fig. S12. Clustering of all the genes controlled by CK.

Fig. S13. Clustering of all the genes controlled by ABA.

Fig. S14. Characterization of mutants at the transcriptional level.

Fig. S15. Mutation in a K+ channel affects root development and ABA responsiveness on nitrate-containing media.

Fig. S16. Venn diagram illustrating the influence of signal combinations on the whole root transcriptome.

Table S1. Matrix of treatments applied in this study.

Table S2. Genome-wide ANOVA results.

Table S3. GeneCloud analysis of the 10 dominant modes of regulation and gene clusters.

Table S4. Best marker genes for the simple signals.

Data File S1. Nitrogen+Hormones.cys Cytoscape file for Fig. 3B.


Acknowledgments: We thank J.B. Thibaut for providing the gork-1 seeds. Funding: This work, especially the modeling procedures, was supported by the French Agence Nationale de la Recherche (ANR) (NitroNet: ANR 11 PDOC 020 01) and the CNRS (PEPS Bio math Info 2012–2013: SuperRegNet). Genome-wide approaches were supported by European-FP7-International Outgoing Fellowships (Marie Curie) (AtSYSTM-BIOL; PIOF-GA-2008-220157) to G.K. C.C. work was supported by the Labex NUMEV (SuperRegNet). Work on networks was supported by NIH R01-GM032877 to G.M.C. Results on root architecture responses were supported by NSF MCB-0929338 to G.M.C. Bioinformatics analysis was supported by the VirtualPlant platform ( developed under NSF DBI-0445666 to G.M.C. The work of D.R. at New York University was supported by the International Fulbright Science & Technology Award for Outstanding Foreign Students, and her work in the Busch laboratory was supported by an EU FP7 COFUND PLANT FELLOWS grant. A.M. was supported by the French ANR (ANR-JCJC-NUTSE). The collaboration by G.K., S.R., and G.M.C. was supported by the creation of an international-associated laboratory (LIA-CoopNet) from the CNRS. Author contributions: G.K. and G.M.C. designed the project. D.R., A.M., M.P., S.R., G.J.K., D.S., and G.K. performed the experiments. D.R., A.M., M.P., S.R., C.C., and G.K. analyzed the data. B.L., S.R., K.D.B., W.B., G.M.C., and G.K. contributed to the study design during the project course. D.R. and G.K. wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Transcriptomic data are available through the GEO platform (accession no. GSE71737). Arabidopsis mutants can be retrieved from G.K. or D.R.
View Abstract

Stay Connected to Science Signaling

Navigate This Article