Transcriptome Responses to Combinations of Stresses in Arabidopsis

Biotic and abiotic stresses limit agricultural yields, and plants are often simultaneously exposed to multiple stresses. Combinations of stresses such as heat and drought or cold and high light intensity have profound effects on crop performance and yields. Thus, delineation of the regulatory networks and metabolic pathways responding to single and multiple concurrent stresses is required for breeding and engineering crop stress tolerance. Many studies have described transcriptome changes in response to single stresses. However, exposure of plants to a combination of stress factors may require agonistic or antagonistic responses or responses potentially unrelated to responses to the corresponding single stresses. To analyze such responses, we initially compared transcriptome changes in 10 Arabidopsis ( Arabidopsis thaliana ) ecotypes using cold, heat, high-light, salt, and ﬂ agellin treatments as single stress factors as well as their double combinations. This revealed that some 61% of the transcriptome changes in response to double stresses were not predictable from the responses to single stress treatments. It also showed that plants prioritized between potentially antagonistic responses for only 5% to 10% of the responding transcripts. This indicates that plants have evolved to cope with combinations of stresses and, therefore, may be bred to endure them. In addition, using a subset of this data from the Columbia and Landsberg erecta ecotypes, we have delineated coexpression network modules responding to single and combined stresses. foliar spray application of a bacterial elicitor FLG peptide g22) performed in environmentally controlled chambers (RISØ DTU National Laboratory). A pilot study using Col and L er ecotypes performed to identify sublethal doses of combined stress treatments. This identi ﬁ ed an optimal period of 3 h before the onset of visible phenotypic responses such as wilting. To ensure independence between biological replicates, the stress treatments and plant growth were done in three independent batches. Each stress treatment lasted 3 and was done on 3-week-old plants. The high-salinity treatments were performed by soil irriga- tion with 100

Plants are often simultaneously exposed to various biotic and abiotic stresses in their natural or agronomic habitats (Ahuja et al., 2010). Roughly 300 cellular stress genes are conserved in all organisms to defend or repair vital macromolecules against environmental factors (Kültz, 2005). However, stress response genes also evolve rapidly as organisms adapt to changing environments. Thus, antifreeze proteins evolved separately in different phyla (Cheng, 1998), and roughly half of the osmoresponsive genes in the model plant Arabidopsis (Arabidopsis thaliana) are plant specific (Rabbani et al., 2003). Because biotic and abiotic stresses reduce harvest yields, considerable research has aimed to understand the responses of model plants and crops to single stresses (for review, see Hirayama and Shinozaki, 2010;Chew and Halliday, 2011). This work has identified sets of canonical response genes induced by heat, cold, osmotic, or high-light stresses (Kreps et al., 2002;Seki et al., 2002;Rizhsky et al., 2004;Oono et al., 2006;Kleine et al., 2007;Hannah et al., 2010;González-Pérez et al., 2011) and in response to pathogen infection and exposure to pathogen-associated molecular patterns (Navarro et al., 2004;Nielsen et al., 2007). It has also revealed that plant responses to different stresses are coordinated by complex and often interconnected signaling pathways regulating numerous metabolic networks (Nakashima et al., 2009). Nonetheless, apart from a notable study on the effects of simultaneous drought and heat stress (Rizhsky et al., 2004), the effects of stress combinations have been little studied (Mittler, 2006;Atkinson and Urwin, 2012). Therefore, further work is needed if we wish to understand the full complement of stress responses by comparing data on single stresses with data on multiple stress responses. Such data will be relevant to agronomy (Oerke et al., 1994) and provide tools to answer basic questions about signaling "cross talk" in systems biology (Mundy et al., 2006).
Whole-genome expression profiling with microarrays is a useful tool to monitor changes in transcript levels and thereby gene expression in response to stresses and other factors (Seki et al., 2009). Most such studies have used Arabidopsis as a model because of its amenability to subsequent forward and reverse genetic analyses (Somerville and Koornneef, 2002). Robust algorithms have been developed for high-throughput microarray data to decipher global biological processes and to generate testable biological hypotheses (Harr and Schlötterer, 2006;Chawla et al., 2011). For example, lists of transcripts differentially responding to different stresses can be generated and ranked by biological criteria (Dudoit et al., 2002;Nielsen et al., 2007). Network-based algorithms can successfully deal with some of the complexities of biological and other physical systems (Barabási and Oltvai, 2004). Different methods and algorithms have been developed (Chawla et al., 2011) to construct major types of biological networks, including gene-metabolite, protein-protein interaction, transcriptional regulatory, gene regulatory, and coexpression networks (Yuan et al., 2008). Coexpression or coregulation of genes may implicate them in similar biological processes, such that individual modules of genes can be attributed to specific processes. A primary assumption is that, in such networks, strongly coregulated or coexpressed groups of genes participate in similar biological processes, such as signaling or metabolic pathways (Williams and Bowles, 2004). For example, a study by Weston et al. (2008) showed that a coexpression network-based analysis could delineate population-level, adaptive physiological responses of plants to abiotic stress. In addition, by meta-analysis of microarray data and other publicly available information (Mentzen and Wurtele, 2008), modular coexpressionbased analysis can dissect regulon organization in the Arabidopsis genome by identifying functional modules that share a similar expression profile across multiple spatial, temporal, environmental, and genetic conditions.
We conducted a large-scale microarray experiment to analyze plant responses to multiple, concurrent stresses and to identify the levels and functions of stress regulatory networks. To this end, 10 ecotypes of the model Arabidopsis were subjected to five individual stress treatments and six combinations of these stress treatments under the same growth and experimental conditions. Here, we present and analyze this homogeneous data set of ecotype responses to single and combined stresses. Importantly, our analysis shows that, when two stresses were combined, an average of 61% of the transcripts responded in modes that could not be predicted from individual single stress treatments. In addition, only a minor fraction (6%) of the transcripts exhibited antagonistic responses to stress combinations under which the plants apparently must prioritize between the responses. Given the novelty of the responses we uncovered, we explored the modular organization of transcription networks using weighted gene coexpression network analysis (WGCNA; Zhang and Horvath, 2005). This permitted us to identify stress-responsive modules and potentially key regulatory genes to further understand plant responses to multiple stresses.

Transcript Profiling of Stress Treatments
To investigate the effects of five single stress treatments (cold, high light, salt, heat, and flagellin [FLG]) and six combined stress treatments (cold and high light, salt and heat, salt and high light, heat and high light, heat and FLG, and cold and FLG) on global transcript levels in Arabidopsis, labeled RNA was hybridized in triplicate to Arabidopsis NimbleGen ATH6 microarrays (Supplemental Table S1). A total of 210 arrays covering stress experiments from 10 different ecotypes (Columbia [Col], Landsberg erecta [Ler], Columbia 24, Cape Verde Islands, Kashmir, Antwerpen, Shakdara, Kyoto 2, Eringsboda, and Kondara) were hybridized, and three arrays were identified as outliers (data not shown) and removed to achieve a total of 207 arrays. The array contained probes for 30,380 transcripts for which significant changes in levels were determined by comparing single or double stress treatments with ecotype-matched controls and used in further analyses.

Responses to Single Stresses and Intraspecific Variation
We initially compared the Col and Ler transcript responses of the single stress treatments with a benchmark set of responses to similar stress treatments. Despite the fact that the benchmark sets are composed from previously published studies using various ecotypes and experimental setups, there was a good overall overlap (Supplemental Table S2), including key stressresponsive transcripts (Supplemental Table S3). These results indicate that the individual stress treatments we applied had effects similar to those described previously in other analyses and that our individual stress treatments were appropriate.
We then investigated how the ecotypes responded to single stress treatments to identify ecotype differences (Supplemental Fig. S1). In general, the responses of the 10 ecotypes were highly correlated, except for Col, which behaved as an outlier for responses to heat, salt, or high light. Such intraspecific variation in responses to environmental stimuli have been well documented in a number of plant species, including Arabidopsis ecotypes (Koornneef et al., 2004). As Col and Ler are widely used for Arabidopsis research, we focused on these two ecotypes to obtain consensus results. Using the Col and Ler transcript sets, we compared transcript overlap and congruency between each of the responses. This revealed low overlap between abiotic and biotic transcript stress responses. Among the single abiotic stresses, cold and high light were most similar (33% transcript overlap) with very congruent responses such that 87% of the transcripts responded in the same direction (increased or reduced levels; Supplemental Fig. S2). In contrast, considerable dissimilarity was observed between responses to salt and cold or salt and heat single stresses, while responses to salt and high light were more congruent.

Responses to Combined Stresses
In principle, when plants are exposed to combined stresses, their responses to the single stresses must be modulated to produce a combined response. The responses of a given transcript to two single stresses may be neutral, agonistic, antagonistic, or unrelated, and the response to the combined stresses may be a combination of such responses. However, as described below, the response of transcripts to combined stresses is not easily predictable. To describe these responses, we clustered significantly responding transcripts (Supplemental Table S4) from the two single stresses and the combined stress treatments to predefined expression profiles (see "Materials and Methods") and defined five transcript response behaviors or modes: "combinatorial," "canceled," "prioritized," "independent," and "similar" (Fig. 1). Transcripts in combinatorial mode have similar responses to the two single stresses and different responses to the combined stresses. Transcripts in canceled mode respond differently to the single stresses and are returned to control levels in response to the combined stress treatment. Transcripts in prioritized mode respond differently to the single stresses and remain at one of these levels in response to the combined stresses. Importantly, these three modes represent transcript responses or regulatory modes between two stress factors that cannot be predicted from single stress experiments. In contrast, transcripts in independent mode, whose regulation patterns do not respond to the addition of a second stress, as well as transcripts in similar mode for the two single and the combined stresses might be more readily identified from single stress experiments.
Comparison of the responses of single versus combined stresses revealed that, on average, 61% of the transcripts responded in a mode (combinatorial, canceled, or prioritized) in which the two single stress responses interact such that the response to the combined stresses cannot be predicted from the single stress experiments alone ( Fig. 2A). The extents of these interacting modes were dependent on the particular stresses applied and ranged from 49.3% of the transcripts in the heat and FLG experiment to 73.8% in the salt and heat double stress experiment (Fig. 2B). The majority of the interacting transcript responses were canceled or combinatorial, with averages of 29.3% and 24.7%, respectively. However, on average, only 6.8% of the transcripts responded in the prioritized mode, in which the plant must decide between antagonistic responses. This indicates that responses to multiple stresses involve relatively few, oppositely responding transcripts. The experiment with the smallest fraction of transcripts responding in prioritized mode is the cold and high light experiment (3.0%; Fig. 2B). This is in good agreement with the responses to the two single stress treatments, which share a large overlap in transcripts and thus exhibit very congruent responses (Supplemental Fig. S2). In contrast, the combined salt and heat treatments show the highest level of prioritized transcripts (12.1%).
The independent and similar response modes, in which stress responses apparently do not interact upon combined stress exposure, constituted on average 39% of the transcript responses, and the number of transcripts regulated in the independent mode generally was the larger fraction of these (28%). The response to combined salt and heat treatment was the most oppositely directed in that it had the highest level of unpredictable responses (combinatorial, canceled, or prioritized), with the fewest number of transcripts in similar mode (1.5%), whereas the combined cold and high light treatments had the largest number of transcripts in similar mode (18.6%). This corresponds well to the congruency and similarity of the single stress treatments and the levels of prioritized transcripts found above for these double stress combinations.
We then investigated whether transcripts of the particular response modes could be associated with biological functions via their corresponding, significant GO terms (Supplemental Fig. S3; Supplemental Table S5). Several patterns appear clear, as three of the response modes are associated with sets of GO terms specific to them. The independent mode is associated with chloroplastic or photosynthetic terms including thylakoid membrane organization and responses to cyclopentenones and high light intensity. This indicates that the effects of high-light treatment may contain a specific set of transcripts that are not influenced by the other stress treatments and may solely be associated with light treatment. The canceled mode is primarily associated with terms related to secondary metabolism (anthocyanin, indoleacetic acid, phenylpropanoid, etc.) and growth regulation (ethylene and auxin responses). This may indicate that different stresses promote different secondary metabolite pathways and differentially affect growth in response to auxin and ethylene. The combinatorial mode is primarily associated with defense terms (systemic acquired resistance, programmed cell death, salicylate biosynthesis, etc.), which Plant Physiol. Vol. 161, 2013 1785 may reflect the intersecting pathways regulating defense against varied pathogens as well as alterations in defense-related programs, such as endoplasmic stress, in response to abiotic stresses. Not surprisingly, other, more general GO terms, such as response to hormones and water deficit, are shared among five and four of the response modes. We note that the transcript response mode assignments were stable using different transcript significance thresholds (Supplemental Fig. S4) and verified the profiles and assignments for six transcripts by quantitative PCR (qPCR; Supplemental Fig. S5).
To further investigate the dominance in a pair of stress treatments, we quantified the extent of the influence regulation that each stress imposed on the other in a stress combination (Fig. 3). Interestingly, this indicated that the single stress responses were not always dominated by another stress or vice versa. For example, transcripts that significantly responded to the single high-light treatment were regulated to a lesser extent when combined with cold or salt than when high light was combined with heat. However, transcripts responding to single heat treatment were regulated less in this combination (high light and heat) than when combined with a different stress than high light. In contrast, the transcript response to FLG treatment alone in both cases was regulated to a lesser extent than the abiotic transcript responses (heat and cold), whereas the cold and salt single stress transcript responses were regulated to a greater extent in their stress combination experiments.

WGCNA
To group these stress-responding transcripts into stress regulatory modules, we performed a WGCNA on the 2,236 most differentially regulated transcripts with a thresholding power of 7, which was the lowest power for a good fit of the scale-free topology index (see "Materials and Methods"; Supplemental Fig. S6). This identified nine significant coexpression modules (Table I; Supplemental Table S6). These describe Figure 1. Clustering of transcripts to predefined expression profiles generating the transcriptional response modes. For each stress combination, transcript sets were created by the union of the 500 most significant transcripts for each single stress and for the combination. These transcripts were clustered to 20 predefined expression profiles, each categorizing a potential expression change that may occur when multiple stresses are applied. Each transcript was assigned to the profile with the highest Pearson correlation coefficient. Boxes at left represent the 20 predefined expression profiles described by the transcript expression pattern for stress 1 (S1), stress 2 (S2), and the combination of stress 1 and stress 2 (S1+S2). The dotted line represents transcript expression with no change compared with the control. In the boxes at right, each column represents the union of the two single stresses and the combined stress for a given double stress experiment. In each box, the value and color represents the percentage of transcripts that correlate with the particular predefined expression profile (green is higher). The transcriptional response modes are composed of a given set of the predefined expression profiles as indicated on the right: combinatorial, similar levels in the two individual stresses but a different response to combined stresses; canceled, transcript response to either or both individual stresses returned to control levels; prioritized, opposing responses to the individual stresses and one stress response prioritized in response to combined stresses; independent, response to only one single stress and a similar response to combined stresses; similar, similar responses to both individual stresses and to combined stresses. HL, High light. transcript sets that have similar response profiles throughout the sample series of Col and Ler ecotypes. These modules were associated with stress treatments by singular value decomposition (Langfelder and Horvath, 2007), which significantly associated modules 2, 6, 7, and 8 to abiotic, cold and high light, cold, and biotic stress treatments, respectively (Fig. 4). In addition, GO enrichment of transcripts in each module (Table I) was performed. Only 104 transcripts could not be placed in any of the modules and were assigned to module 10. Furthermore, the network node degree distribution of the WGCNA showed scale-free behavior and could identify highly connected transcripts inside the modules (Supplemental Fig. S7; Supplemental Table S6). Some of these are known to be central stress regulators, in particular a number of transcription factors (TFs). Four of the most significant, functionally associated modules are described below.

Module 8: Biotic Stress Response Module
Module 8 showed significant association with biotic stress (FLG) and with combinations of biotic and temperature stresses. It did not show any association with single cold, heat, salt, or high-light stress. This module included 72 annotated TFs (Guo et al., 2005). Some of these TFs with higher connectivity within the module were WRKY6 (At1g62300; Robatzek and Somssich, 2002)  Module 7 showed significant association with cold, high light, combined cold and high light, and combined cold and FLG. It had very little association with combined salt and high light and none with single FLG treatment. The latter indicates that combined biotic stress (FLG) and cold generates a unique stress response pattern different from the single stress treatments. Based on their connectivity within the module, two TFs apparently important in cold temperature responses are PRR7 (At5g02810; Salomé et al., 2010) and HMGB2 (At4g23800; Kwak et al., 2007). PRR7 has been implicated as a morning loop component in temperature compensation, while HMG2 has been shown to be induced by cold treatment. Other highly connected transcripts in this module with functional associations to cold stress include COR47 (At1g20440), PGM (At5g51820), RSA4 (At5g01410), LTI30 (At3g50970), HVA22D (At4g24960), and ERD7 (At2g17840; Oono et al., 2006).

Module 6: Cold and High Light Stress Response Module
Module 6 showed significant association with two independent abiotic treatments, cold and high light. The gene enrichment results clearly reflect these associations ( Table I). The presence of some previously reported cold response regulators such as CBF1 (At4g25490), CBF2 (At4g25470), and DREB1A (At4g25480) clearly implicates this module in responses to cold. The module includes 17 other TFs, and some of these that are highly connected are IAA19 (At3g15540), At2g46670, APRR9 (At2g46790), APRR5 (At5g24470), ATHB-2 (At4g16780), CCA1 (At2g46830), HFR1 (At1g02340), PIL1 (At2g46970), as well as a MYB and homeodomain-like protein (At3g10113), a basic helix-loop-helix protein (At3g21330), and three zinc finger proteins (At1g73870, At5g48250, and At5g44260). Most of these TFs have clear functional connections to temperature-and light-dependent developmental programs. For example, APRR5 and APRR9, CCA1, and PIL1 are involved in temperature compensation in circadian rhythms (Salomé et al., 2010), HFR1 regulates a phytochrome A-dependent photomorphogenesis pathway (Yang et al., 2009), ATHB-2 regulates photomorphogenesis and shade avoidance (Steindler et al., 1999) in part by modulating auxin-responsive growth mediated by IAA19 (Tatematsu et al., 2004), and At5g48250 appears to be a target of Flowering Locus C-independent effects of the autonomous floral pathway. The module also has some association with combined cold and FLG but no association with single FLG treatment. This again highlights the fact that the interaction of a biotic stress factor (FLG) with cold was unique.

Module 2: Abiotic Stress Response Module
Module 2 exhibited significant association with both single and combined abiotic stresses such as single Response to abiotic stimulus, cellular response to red or far-red light, circadian rhythm, response to radiation, shade avoidance, response to cold, response to hormone stimulus 7 69 Response to cold, response to blue light, cold acclimation, auxin homeostasis, response to far-red light, cellular response to carbohydrate stimulus, sugar-mediated signaling pathway, response to nonionic osmotic stress, response to abscisic acid stimulus, hyperosmotic salinity response, detection of gravity 8 907 Response to biotic stimulus, response to abiotic stimulus, multiorganism process, response to bacterium, response to heat, response to wounding, response to fungus, response to oxidative stress, response to light intensity, innate immune response, response to jasmonic acid stimulus, response to cold, indole glucosinolate metabolic process, flavonol metabolic process, host programmed cell death, response to hormone stimulus, salicylic acid metabolic process, response to ethylene stimulus, posttranslational protein modification, response to ozone, lignin metabolic process 9 166 Response to stress, electron transport or energy pathways, cell organization and biogenesis 10 104 These transcripts were not placed in any of the modules heat, high light, and salt, combined salt and high light or salt and heat, and slight association with cold and high-light. The module did not exhibit any significant association with FLG treatments. Functional enrichment analysis (Table I)   internode growth, which are regulated by auxin and gibberellin. Little is known of the functions of other TFs in this module, including DREBA-4 (At5g52020), SRS6 (At3g54430), emb2746 (At5g63420), HMG1/2-like (At1g76110), DOF (At5g65590), ZN-finger (At2g37430), SMAD/FHA (At2g21530), and two BLH genes (At2g16400 and At5g5091).

DISCUSSION
This study was designed to investigate combinations of stresses that mimic harsh environmental conditions that occur in the field. Therefore, we examined high light intensity that may occur due to diminished ozone layer protection, in combination with low or high temperature, as well as saline irrigation combined with high temperature. Furthermore, we were interested in the interaction of abiotic stresses such as low or high temperature combined with biotic stresses such as the immune response to pathogens induced by the FLG elicitor. Here, we provide a benchmarked set of transcript profiles responding to such single and combined stresses. We initially determined that the transcript responses we detected to five single stresses were similar to those of a benchmark set of responses to similar stresses described by others. This meant that the stresses we applied were comparable to other studies and permitted us to then analyze the transcript profiles responding to combinations of these stresses. We note, however, two potential limitations to the stress applications that others and we have used. First, due to logistics and costs, we harvested rosettes after stress applications at a single time point based on previous studies. The assessment, therefore, will not detect the temporal dynamics of single stress responses or of the interactions between combinations of stresses. Nonetheless, we chose the sampling time point before the onset of visual stress symptoms to attempt to detect responses caused specifically by the environmental insult and not systemic responses that occur "downstream" or as indirect consequences of the specific stress combinations. Second, the assessment does not detect the relative intensities of single stresses. Despite these caveats, we expect that our profiles, as the largest robust such data set at present, can be productively mined by other researchers.
The second aim of this study was to analyze the response profiles to identify the behaviors or modes of regulation of sets or modules of transcripts in response to combined stresses. To these ends, we used two types of analyses to attempt to capture the range of responses displayed in the data. Our analysis of transcript behaviors or modes was based on clustering the top 500 most significantly responding transcripts for each single stress and for the stress combinations. Five transcript response modes (combinatorial, canceled, prioritized, independent, and similar) were identified that describe potential transcript regulatory complexity and assign predictabilities to the responses of individual transcripts. Importantly, the combinatorial, canceled, and prioritized response modes on average constitute 61% of the total transcripts, and these modes cannot be predicted from the corresponding single stress experiments. That the majority of the transcript responses are not predictable when two stresses are combined points to the limitations of attempts to delineate common stress responses or points of cross talk between signaling pathways during multiple stresses by simply identifying overlapping sets of genes that are regulated by both stresses (Kreps et al., 2002;Mentzen and Wurtele, 2008;Carrera et al., 2009). Additionally unpredictable are transcripts that respond only to the combined stresses and not to either individual stress. For example, using the 500 most significant transcripts from each stress experiment, we found that 55.8% to 79.6% of the transcripts regulated in the double stress experiment were not among the most significant transcripts in the corresponding single stress experiments. Therefore, these potentially novel stress-related transcripts will be absent in analyses using only single stress response experiments.
Overall, the most abundant transcript response modes are canceled, independent, and combinatorial, which together on average include 85% of the total transcripts. The combinatorial response (27.1%) of transcripts with similar responses to two single stresses and different responses to the combined stresses is indicative of the level of interaction between responses to different stresses. For example, combined heat and FLG treatment has the lowest level of combinatorial transcripts (7.1%) and the highest level of independent transcripts (39.7%) along with a relatively high level of prioritized transcripts (9.1%). The low overlap in transcripts between the abiotic and biotic single stress experiments (approximately 5%; Supplemental Fig. S2) could indicate that the early transcriptional responses to an innate immune elicitor versus an extreme physical change in the environment target different transcript sets.
The independent mode (28.2%) contains transcripts that are regulated in either of the single stresses and whose regulation is maintained without interference from the other stress. These transcripts define the proportion of the responses to combined stresses that are not shared by, and do not interfere with, the responses to the single stresses. For example, responses to combined cold and high light have the lowest level of independent (18.5%) and prioritized (3.0%) transcripts. This is in line with the highest (33%) and most congruent (0.87%) transcript overlap between responses to cold or high light singly and the highest similar response mode of combined cold and high light (18.6%). In addition, 60% of the transcripts respond to combined cold and high light in canceled or combinatorial mode, indicative of strong regulatory interactions between these two stresses. This may mirror the ecology of the temperate Arabidopsis ecotypes used here, in that combined cold and highlight stress is common in temperate regions (Ivanov et al., 2012).
The most common mode is canceled (30.6%), in which transcripts responding oppositely to single stresses return to control levels in response to the combined stresses. For example, responses to salt or heat stress are significantly dissimilar (0.31%) and the transcript response modes to combined salt and heat have the highest level of prioritized transcripts (12.1%) of all the combined stress experiments. In addition, combined salt and heat have the lowest level of similar transcripts (1%) compared with an average of 11.2% for all combined stress experiments and have very high fold changes compared with the single salt or heat stresses (Fig. 3). The latter indicates that, under the conditions we applied, responses to heat stress largely dominate responses to salt stress. Despite the caveats noted above concerning stress durations and intensities, this result suggests that adaptation to combined salt and heat stresses is more difficult than adaptation to the other combined stresses assessed here. Nonetheless, such a conclusion may only apply to temperate plants with similar ecologies to Arabidopsis, given the extent of adaptive diversity in related plants such as the halophyte Thellungiella salsuginea (Wu et al., 2012).
Our description of the transcript sets or modules responding to combined stresses employed weighted gene coexpression analysis (Langfelder and Horvath, 2008) to hierarchically cluster modules with similar transcript profiles. The resulting weighted transcriptional networks exhibited scale-free, modular topology, and the clustering resulted in nine significant modules. Eigengene significance-based analysis showed that among the nine modules, four have significant associations with different biotic and abiotic stresses. Module 2 appears to be associated with abiotic stress, while module 6 exhibits cold-and high light-associated response signatures. Modules 7 and 8 have significant association with biotic stress (FLG) and with combinations of biotic and temperature stresses. Transcripts in modules 7 and 8 may be useful for addressing agronomic problems, such as reduced crop productivity due to new pathogen invasions coupled with temperature stress in predicted scenarios of global climate change. Using network connectivity and gene significance measures, supported by existing information from Arabidopsis TF databases, we have identified a number of TFs as targets for future translational experiments to engineer increased stress resistance in crops.

Plant Stress Treatments
Arabidopsis (Arabidopsis thaliana) plants of ecotypes Col, Ler, Columbia 24, Cape Verde Islands, Kashmir, Antwerpen, Shakdara, Kyoto 2, Eringsboda, and Kondara were subjected to the following stress treatments: salt, cold, heat, high light, salt + heat, salt + high light, cold + high light, heat + high light, as well as FLG (flg22 peptide), cold + FLG, and heat + FLG. In addition, FLGtreated plants were grown with two control conditions, control and control + Silwet (control for the effect of Silwet detergent used for FLG application). Stress treatments were selected from previous studies (Kreps et al., 2002;Seki et al., 2002;Kilian et al., 2007) and microarray experiments compiled at www.weigelworld.org/resources/microarray/AtGenExpress. Combinations of high light (800 mmol photons m 2 2 s 2 1 ), cold (10°C), heat (38°C), high salinity (100 mM NaCl), and foliar spray application of a bacterial elicitor (20 mM FLG peptide flg22) were performed in environmentally controlled chambers (RISØ DTU National Laboratory). A pilot study using Col and Ler ecotypes was performed to identify sublethal doses of combined stress treatments. This identified an optimal period of 3 h before the onset of visible phenotypic responses such as wilting. To ensure independence between biological replicates, the stress treatments and plant growth were done in three independent batches. Each stress treatment lasted 3 h and was done on 3week-old plants. The high-salinity treatments were performed by soil irrigation with 100 mM NaCl solution. In order to saturate the soil, irrigation with the saline solution started at the end of the light period the night before collection and was refilled at the onset of the combined treatment. For the cold + high-light treatment, heat from three sodium lamps was displaced by circulating fans and a plexiglass shield, and ambient plant temperature was maintained by ice trays and monitored at 10°C with an infrared thermometer (ThermaTwin TN408LC). To reduce the effects of circadian rhythmicity, treatments were performed 5 h after chamber dawn. After stress treatments, leaves were collected and frozen in liquid nitrogen.

RNA Isolation and Microarray Hybridizations
Total RNA samples were isolated (RiboPure kit; Ambion), and reverse transcription of mRNA was performed from total RNA using a doublestranded cDNA synthesis kit (Superscript; Invitrogen). The copy DNA (cDNA) obtained was subsequently labeled with Cy3, and the product was precipitated using NimbleGen kits according to the NimbleGen Gene Expression protocol for microarrays. Four micrograms of the labeled products was loaded onto microarrays, hybridized overnight, and washed in the NimbleGen Wash Buffer Kit following the NimbleGen protocol. Scanning was performed on a Roche 2-micron scanner, and the images were analyzed with the NimbleScan software. The microarray used was the Arabidopsis NimbleGen 12-plex chips using the ATH6 build (Gene Expression Omnibus no. GPL16226) in a Latin Square design with four independent probes per transcript. A total of 210 arrays were hybridized.

Microarray Data Preprocessing
Data were imported into R (R Core Team, 2012) using the oligo (Carvalho and Irizarry, 2010) and pdInfoBuilder (Falcon et al., 2012) packages using the AgilentAT6 build. If more than one scan was available for an array, the best scan was selected using singular value decomposition as the array with the lowest residuals. The data were normalized using quantiles , and three outliers were removed by comparing the arrays using Pearson correlation coefficient and singular value decomposition plots, giving a total of 207 arrays for the analysis. Expression indexes were calculated using robust multiarray average . All statistical comparisons between experiments were performed by Student's t test using the normalized log 2 transcript expression indexes. All treatments were compared with the control experiments except treatments including FLG, which were compared with control and Silwet samples. Transcript annotation was acquired from The Arabidopsis Information Resource 10 (Lamesch et al., 2012) using the biomaRt data-mining tool (Guberman et al., 2011). The microarray data are available at the Gene Expression Omnibus with the record GSE41935.

Benchmarking
Because the benchmarking gene sets are derived from various experiments, ecotypes, and sources, we used all transcripts from the Col and Ler single stress experiments with P # 0.01 as input for the benchmarks. Additionally, we used the top 500 most significant transcripts from each treatment from comparisons using all ecotypes as a single group. The benchmarking gene sets were derived from the following sources (Ashburner et al., 2000;Huala et al., 2001;Navarro et al., 2004;Thimm et al., 2004;Kilian et al., 2007;Kleine et al., 2007;Kant et al., 2008;Papdi et al., 2008;Shameer et al., 2009;González-Pérez et al., 2011;Less et al., 2011;Avin-Wittenberg et al., 2012;Kilian et al., 2012) and are listed in Supplemental Table S2.

Transcriptional Response Modes
For each stress combination, transcript sets were created by the union of the top 500 most significantly responding transcripts for each single stress and for Plant Physiol. Vol. 161, 2013 1791 the combination of the two stresses. Thereafter, the transcript sets were clustered using Pearson correlation coefficient to 20 predefined expression profiles, each categorizing a potential expression pattern that may occur upon multiple stress application. The predefined expression profile with the highest Pearson correlation coefficient was selected for each transcript. The transcriptional response modes (combinatorial, canceled, prioritized, independent, and similar) were a multiset of several predefined expression profiles. Association of transcriptional response modes with GO terms was performed using GO-slim from The Arabidopsis Information Resource as of December 12, 2012 (Berardini et al., 2004) using Fisher's exact test with a significance threshold of 10 25 after Bonferroni correction. To simplify the network for GO terms with identical connectivity, we only kept the most specific term.
qPCR Verification of Microarray Transcripts qPCR was performed using the Brilliant II SYBR Green one-step kit (Agilent Technologies) with 10 pmol of each primer and 12.5 ng of total cDNA in 10 mL, and the reactions were run on a CFX 96 Thermocycler (Bio-Rad). For this verification, biological triplicate Col samples were used and relative log 2 expression was determined using ACT2 (AT3G18780), which was determined to be highly expressed with minimal variation across the different treatments in the microarray data. All primer efficiencies were within 100% 6 2%, and expression levels were calculated assuming 100% efficiency. Primers used and an agarose gel of PCR products matching the expected product sizes are shown in Supplemental Table S7.

WGCNA
A weighted gene coexpression network was constructed using the R package WGCNA (Langfelder and Horvath, 2008) using a united list of significant transcripts (P # 0.01) generated from Student's t tests between control and treated plants in the Col and Ler ecotypes. A total of 2,236 transcripts that responded to at least two of the stress treatments were used to construct the weighted network from the normalized expression data by transforming the pairwise gene correlation matrix to a weighted matrix with a scaling factor beta (b = 7) and using the assumption that biological networks are scale free. Here, weight represents the connection strength between gene pairs (Zhang and Horvath, 2005). To minimize the effects of intrinsic noise in highthroughput transcriptomic data, we transformed the adjacency into a topological overlap matrix (Yip and Horvath, 2007). The tree was created using hierarchical clustering, and the dynamic tree-cut algorithm was used to identify modules with similar expression patterns. A minimum module size of 30 and a height cut of 0.25 corresponding to a correlation of 0.8 were used to merge similar transcripts. The module eigengenes were used to define measures of module membership (at the significance level P # 0.001), intramodular connectivity, and gene significance (Langfelder and Horvath, 2007). Intramodular connectivity of transcripts was used to identify hubs in the modules and was measured by computing the whole-network connectivity (kTotal), the within-module connectivity (kWithin), and the outgoing connectivity (kOut = kTotal 2 kWithin). BiNGO (Maere et al., 2005), an opensource Java tool, was used to determine which GO terms were significantly overrepresented in our module transcript lists (P values were Bonferroni corrected).

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S1. Comparison of ecotype responses to the single abiotic stresses (salt, cold, heat, and high light).
Supplemental Figure S2. Overlap congruence between responses to individual stress treatments based on comparisons between the top 500 transcripts of each single treatment.
Supplemental Figure S3. Network of transcriptional response modes and their associated GO terms.
Supplemental Figure S4. Comparison using different thresholds for generating the transcriptional response modes shows consistent profiles.
Supplemental Figure S5. Microarray and qPCR transcript profiles for six transcripts verifying the transcriptional data and profile assignments.
Supplemental Figure S6. Hierarchical clustering, soft threshold, and clustering dendrogram of transcripts of the WGCNA network.
Supplemental Figure S7. Scale-free behavior of network node degree distribution of the WGCNA network.
Supplemental Table S1. Overview of the experimental setup with regard to ecotypes, stress treatments, and biological replicates.
Supplemental Table S2. Benchmark of single stress experiment data versus previous single stress experiments showing overlap, including the reference gene sets used for benchmarking.
Supplemental Table S3. Selection of the key single stress benchmark genes identified in benchmarking.
Supplemental Table S4. Top 500 regulated transcripts for each of the stress treatments from the combined analysis of Col and Ler ecotypes.
Supplemental Table S5. P values of the GO terms and transcriptional response modes for overrepresentation used to build the network.
Supplemental Table S6. WGCNA module membership and connectivity for each of the included transcripts.
Supplemental Table S7. Primers used for qPCR and agarose gel showing PCR products with expected sizes.